gauss_plume_cartesian_func Function

    function gauss_plume_cartesian_func(x_s,y_s,z_s,cos_val,sin_val,x_r,y_r,z_r,ay,by,az,bz,sig_y_0,sig_z_0,delta)

        implicit none
        real x_s,y_s,z_s,x_r,y_r,z_r
        real ay,by,az,bz,sig_y_0,sig_z_0,delta
        real gauss_plume_cartesian_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_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))

            !write(*,*)sig_y_0,sig_y,sig_z_0,sig_z
            !gauss_plume_cartesian_func=1./(2.*pi*sig_y*sig_z)*exp((-y**2)/2./sig_y**2) &
            !    *(exp((-(z_r-z_s)**2)/2./sig_z**2)+exp((-(z_r+z_s)**2)/2./sig_z**2))
            gauss_plume_cartesian_func=1./(2.*pi*sig_y*sig_z)*exp(-y*y/2./sig_y/sig_y) &
                *(exp(-(z_r-z_s)*(z_r-z_s)/2./sig_z/sig_z)+exp(-(z_r+z_s)*(z_r+z_s)/2./sig_z/sig_z))


        endif

    end function gauss_plume_cartesian_func