gauss_plume_cartesian_integral_func Function

    function gauss_plume_cartesian_integral_func(x_s,y_s,z_s,cos_val,sin_val,x_r,y_r,ay,by,az,bz,sig_y_0,sig_z_0,H1,H2,delta)

        implicit none
        real x_s,y_s,z_s,x_r,y_r
        real ay,by,az,bz,sig_y_0,sig_z_0,H1,H2,delta
        real gauss_plume_cartesian_integral_func
        real sig_y,sig_z,x,y
        real cos_val,sin_val
        real pi,sig_limit
        parameter (pi=3.141592,sig_limit=3.)

        !r=sqrt((x_s-x_r)**2+(y_s-y_r)**2)
        !if (abs(u_s).lt.001) u_s=0.001
        !th=atan(v_s/u_s)
        !if (u_s.lt.0) th=th+pi
        !cos_val=cos(th)
        !sin_val=sin(th)
        x=(x_r-x_s)*cos_val+(y_r-y_s)*sin_val
        y=-(x_r-x_s)*sin_val+(y_r-y_s)*cos_val

        gauss_plume_cartesian_integral_func=0.
        sig_y=sig_y_0+ay*exp(by*log(x))+x*abs(delta)
        if (x.ge.0.and.abs(y).lt.sig_y*sig_limit) then
            sig_z=sig_z_0+az*exp(bz*log(x))

            !gauss_plume_cartesian_integral_func=1./(2.*pi*sig_y)*exp((-y**2)/2./sig_y**2) &
            !    *sqrt(pi/2.)*(erf((z_s-H1)/sqrt(2.)/sig_z)-erf((z_s-H2)/sqrt(2.)/sig_z)+erf((z_s+H2)/sqrt(2.)/sig_z)-erf((z_s+H1)/sqrt(2.)/sig_z))/(H2-H1)
            gauss_plume_cartesian_integral_func=1./(2.*pi*sig_y)*exp((-y*y)/2./(sig_y*sig_y)) &
                *sqrt(pi/2.)*(erf((z_s-H1)/sqrt(2.)/sig_z)-erf((z_s-H2)/sqrt(2.)/sig_z)+erf((z_s+H2)/sqrt(2.)/sig_z)-erf((z_s+H1)/sqrt(2.)/sig_z))/(H2-H1)
        endif

    end function gauss_plume_cartesian_integral_func