z_centremass_gauss_func Subroutine

    subroutine z_centremass_gauss_func(sigma,h,z_pbl,z_c)

        implicit none
        real, intent(in) :: sigma,h,z_pbl
        real, intent(out) :: z_c
        real z_loop(5)
        real H_c
        integer i_loop,i
        real sqrt_2pi,sqrt_2
        real pi
        parameter (pi=3.141592653589793)

        sqrt_2pi=sqrt(2.*pi)
        sqrt_2=sqrt(2.)

        i_loop=5
        z_loop(1)=h;z_loop(2)=-h;z_loop(3)=2.*z_pbl-h;z_loop(4)=2.*z_pbl+h;z_loop(5)=-2.*z_pbl+h
        z_c=0.
        !c_z=0.;c_av=0.
        H_c=z_pbl

        !If the emission height h is greater than the boundary layer height then only allow reflection from the surface
        !and set the top of the integration H_c to infinity
        if (h.gt.z_pbl) then
            H_c=1.e16
            i_loop=2
        endif

        !Reduce the loop size when the reflection from the boundary layer is not important
        if (sigma+h.lt.z_pbl/3.) then
            i_loop=2
        endif

        !Remove this after finished testing
        !i_loop=5

        do i=1,i_loop

            z_c=z_c+sigma/sqrt_2pi*(exp(-0.5*(z_loop(i)/sigma)*(z_loop(i)/sigma))-exp(-0.5*((H_c-z_loop(i))/sigma)*((H_c-z_loop(i))/sigma))) &
                +z_loop(i)/2*(erf((H_c-z_loop(i))/sqrt_2/sigma)+erf((z_loop(i))/sqrt_2/sigma))

            !c_av=c_av+0.5*(erf((z_pbl-z_loop(i))/sqrt_2/sigma)+erf((z_loop(i))/sqrt_2/sigma))/z_pbl

            !c_z=c_z+1./sigma/sqrt_2pi*exp(-0.5*((z-z_loop(i))/sigma)*((z-z_loop(i))/sigma))

        enddo

        !if (sigma.gt.0.9*z_pbl) then
        !c_z=1./z_pbl
        !c_av=1./z_pbl
        !endif

    end subroutine z_centremass_gauss_func