gauss_plume_cartesian_trajectory_func Function

    function gauss_plume_cartesian_trajectory_func(x,y,z_s,z_r,ay,by,az,bz,sig_y_0,sig_z_0,delta)

        implicit none
        real x,y,z_s,z_r
        real ay,by,az,bz,sig_y_0,sig_z_0,delta
        real gauss_plume_cartesian_trajectory_func
        real sig_y,sig_z
        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_trajectory_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_trajectory_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_trajectory_func