uEMEP_set_dispersion_sigma_Kz Subroutine

    subroutine uEMEP_set_dispersion_sigma_Kz(Kz_scheme_in,x_in,sig_z00,sig_y00,sigy_0_subgid_width_scale,sig_z_in,z_emis_loc,h_mix_loc,invL,u_val,z_val,logz0,subgrid_delta,u_star0_in,average_zc_h_in_Kz_flag,n_kz_iterations,sig_y_scaling_factor,sig_z,sig_y,u_zc)

        implicit none

        real, intent(in) :: x_in,sig_z00,sig_y00,sigy_0_subgid_width_scale,sig_z_in,z_emis_loc,h_mix_loc,invL,u_val,z_val,logz0,subgrid_delta(2),u_star0_in,sig_y_scaling_factor
        integer, intent(in) :: Kz_scheme_in,n_kz_iterations
        logical, intent(in) :: average_zc_h_in_Kz_flag
        real, intent(out) :: sig_z,sig_y,u_zc

        integer n_loop,j

        real K_z,L
        real :: z_tau_min=2.,z_tau_max=100.
        real z0,zc
        real :: K_min=0.001
        real l_t,f_t
        real u_star0,u_star0_val,tau
        real min_xy,x
        real phih_temp

        n_loop=n_kz_iterations

        !Limit stable L to this positive L value
        !L=1.e6
        !L=lowest_stable_L
        !if (abs(invL).gt.1./L) L=1./invL
        L=1./invL
        !Why does it not explode? when invL=0?
        if (abs(invL) < epsilon0) L = 1.0e6

        z0=exp(logz0)
        min_xy=(subgrid_delta(1)+subgrid_delta(2))/4.
        !min_x=1.
        !x=max(x_in,min_xy)
        !Set x to this value because it simulates that it has already travelled half a grid to get this sig_z0
        !Which is why sg_z is added only to sig_z00
        x=x_in+min_xy

        !Initialise sig_z
        sig_z=sig_z_in

        !If the emission is above the boundary layer height then set the initial plume guess to its low turbulence form
        !Same values used in the emulator (1/1000 slender plume)
        if (z_emis_loc/h_mix_loc.ge.1.0) sig_z=sig_z00+0.001*exp(1.0*log(x))

        !Set ustar0 for K_z to the value from EMEP
        u_star0=u_star0_in

        zc=z_emis_loc

        !Set zc and K_z at start of plume based on sig_z0
        !call z_centremass_gauss_func(sig_z0,z_emis_loc,h_mix_loc,zc_start)
        !call Kz_func(h_mix_loc,L,u_star0,zc_start,K_min,K_z_start)
        !call u_profile_val_func(zc_start,L,u_val,z_val,h_mix_loc,z0,u_zc_start,u_star0,u_hmix)

        !Put bug back in
        u_zc=u_val

        !Calculate the Lagrangian time scale before the iteration loop using a minimum distance for this to make it non zero on the grid
        !l_t=max(x,min_xy)/u_zc
        tau=0.6*max(min(z_tau_max,z_emis_loc),z_tau_min)/u_star0
        !f_t=1.+tau/l_t*(exp(-l_t/tau)-1.)
        !if (x.lt.200) then
        !write(*,'(a,i,5f)') 'S: ',j,tau,l_t,f_t,u_zc,sig_z
        !write(*,'(a,5f)') 'S: ',tau,u_star0,z_tau_max,z_emis_loc,z_tau_min
        !endif

        !All functions commented out 4 secs
        !All functions included 15 secs
        !Without u_profile 11 secs
        !Commenting out all functions below gives 7 secs
        do j=1,n_loop

            !Calculate centre of mass for the emission height.
            call z_centremass_gauss_func(sig_z,z_emis_loc,h_mix_loc,zc) !4.5 sec

            if (average_zc_h_in_Kz_flag) zc=(z_emis_loc+zc)/2.

            !write(*,'(i,4f)') j,sig_z,z_emis_loc,h_mix_loc,zc

            !Calculate the wind profile
            !call u_profile_val_func(zc,L,u_val,z_val,h_mix_loc,z0,u_zc,u_star0_val,u_hmix) !7secs. Without phim_func calls 5 secs
            call u_profile_neutral_val_func(zc,u_val,z_val,h_mix_loc,z0,u_zc,u_star0_val)

            !write(*,'(i,9f)') j,zc,L,u_val,z_val,h_mix_loc,z0,u_zc,u_star0,u_hmix
            !u_zc=(u_zc+u_zc_start)/2.
            !write(*,'(i,8f)') j,x,u_star0_in,u_star0,zc,z_val,u_zc,u_zc_start,u_val

            !Use the calculated u_star for the dispersion, not the EMEP input
            !u_star0=u_star0_val

            !Calculate K_z at the centre of mass
            if (Kz_scheme_in.eq.2) then
                call TROENKz(zc,h_mix_loc,u_star0,invL,K_min,K_z,phih_temp)
            else
                call Kz_func(h_mix_loc,L,u_star0,zc,K_min,K_z) !5.5 secs
            endIf

            !write(*,'(i,1f)') j,K_z
            !Take average of K_z_start(x=0) and K_z(x)
            !K_z=(K_z+K_z_start)/2.

            !calculate l_t based on the centre of mass wind speed
            l_t=max(x,min_xy)/u_zc
            f_t=1.+tau/l_t*(exp(-l_t/tau)-1.)

            !Calculate sig_z for the next iteration, using the sig_z0 from the neutral plume approximation
            sig_z=sig_z00+sqrt(2.*K_z*l_t*f_t)
            !write(*,'(i,1f)') j,sig_z
            !write(*,*) K_z,l_t,f_t,sig_z,u_zc
            !write(*,*) j,z_emis_loc,zc,h_mix_loc*0.25,sig_z

            !if (x.lt.200) then
            !write(*,'(a,i,5f)') 'S: ',j,tau,l_t,f_t,u_zc,sig_z
            !endif

        enddo

        !Calculate sigma_y at the maximum K, around the average of the emission height and 0.25 of the boundary layer height
        !THis is new 27.10.2018 and not tested. It will reduce sig_y which is OK
        !h_y=(h_mix_loc*0.25+z_emis_loc)/2.
        !call Kz_func(h_mix_loc,L,u_star0,h_y,K_min,K_y)
        !sig_y=sig_y00+min_xy*sigy_0_subgid_width_scale+sqrt(2.*K_y*l_t*f_t)

        !Should change this to what is documented, i.e. 2*sig_z. Need to test
        sig_y=sig_y00+min_xy*sigy_0_subgid_width_scale+(sig_z-sig_z00)*sig_y_scaling_factor

    end subroutine uEMEP_set_dispersion_sigma_Kz