z_centremass_gauss_array_func Subroutine

    subroutine z_centremass_gauss_array_func(sig_norm,h_norm,n_array,zc_array)

        implicit none
        real, intent(in) :: sig_norm,h_norm,n_array
        real, intent(out) :: zc_array(int(n_array))
        real z_loop(5)
        real H_c
        real z_c
        integer i_loop,i,k
        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_norm;z_loop(2)=-h_norm;z_loop(3)=2.-h_norm;z_loop(4)=2.+h_norm;z_loop(5)=-2.+h_norm
        H_c=1.

        !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_norm.gt.1.) then
            H_c=1.e4
            i_loop=2
        endif

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

        !Remove this after finished testing
        i_loop=5

        do k=1, int(n_array)

            z_c=0.

            do i=1,i_loop

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

            enddo
            zc_array(k)=z_c

        enddo

    end subroutine z_centremass_gauss_array_func