gauss_plume_second_order_rotated_reflected_integral_func Function

public function gauss_plume_second_order_rotated_reflected_integral_func(r, ay, by, az, bz, sig_y_0, sig_z_0, z_s, z_pbl, H1, H2)

Arguments

Type IntentOptional Attributes Name
real :: r
real :: ay
real :: by
real :: az
real :: bz
real :: sig_y_0
real :: sig_z_0
real :: z_s
real :: z_pbl
real :: H1
real :: H2

Return Value real


Called by

proc~~gauss_plume_second_order_rotated_reflected_integral_func~~CalledByGraph proc~gauss_plume_second_order_rotated_reflected_integral_func gauss_plume_second_order_rotated_reflected_integral_func proc~uemep_subgrid_dispersion uEMEP_subgrid_dispersion proc~uemep_subgrid_dispersion->proc~gauss_plume_second_order_rotated_reflected_integral_func program~uemep uEMEP program~uemep->proc~uemep_subgrid_dispersion

Source Code

    function gauss_plume_second_order_rotated_reflected_integral_func(r,ay,by,az,bz,sig_y_0,sig_z_0,z_s,z_pbl,H1,H2)

        implicit none
        real r,ay,by,az,bz,sig_y_0,sig_z_0,z_s,z_pbl,H1,H2
        real gauss_plume_second_order_rotated_reflected_integral_func
        real sig_th,sig_z,B,c_y_int,c_z_int
        real order_1,order_2
        real :: correction=2.
        integer k,n_loop
        real z_loop(6)
        real pi
        parameter (pi=3.141592)

        !Corrected for the B**2 falut in the taylor expansion and for the fact that the integral was only half a circle. 20.08.2019
        !Still need to implement reflections
        r=max(0.001,r)
        order_1=1.
        order_2=1.

        sig_th=(sig_y_0+ay*(exp(by*log(r))))/r
        sig_z=sig_z_0+az*(exp(bz*log(r)))

        B=-(sig_th**2)*(bz*(sig_z-sig_z_0)/r/sig_th+by*(r*sig_th-sig_y_0)/sig_z)

        if (z_s.gt.z_pbl.or.z_s+sig_z.lt.z_pbl/3.) then
            n_loop=2
            z_loop(1)=z_s;z_loop(2)=-z_s
        else
            n_loop=5
            z_loop(1)=z_s;z_loop(2)=-z_s;z_loop(3)=2.*z_pbl-z_s;z_loop(4)=2.*z_pbl+z_s;z_loop(5)=-2.*z_pbl+z_s
        endif

        if (sig_z.gt.0.9*z_pbl) then
            c_z_int=1./z_pbl
        else
            c_z_int=0.
            do k=1,n_loop
                c_z_int=c_z_int+sqrt(pi/2.)*sig_z*(erf((H2-z_loop(k))/sqrt(2.)/sig_z)-erf((H1-z_loop(k))/sqrt(2.)/sig_z))/(H2-H1)
            enddo
        endif

        if (B.gt.-1.) then
            c_y_int=1./(2.*pi*sqrt(2.*pi)*r*sig_z*sqrt(1.+B))*erf(pi/(2.*sqrt(2.))/sig_th*sqrt(1.+B)*correction)
        else
            c_y_int=correction/(4.*pi*sig_th*r*sig_z)*(1-order_1*pi**2*(1.+B)*correction**2/(24*sig_th**2)+order_2*pi**4*((1.+B)**2*correction**4/(640.*sig_th**4)))
        endif
        !write(*,*) c_y_int,c_z_int

        gauss_plume_second_order_rotated_reflected_integral_func=c_y_int*c_z_int

    end function gauss_plume_second_order_rotated_reflected_integral_func