subroutine uEMEP_set_dispersion_sigma_PG(invL_in,logz0,sig_z_00,sig_y_00,sigy_0_subgid_width_scale,subgrid_delta,delta_wind,x,sig_z,sig_y,sig_z_0,sig_y_0)
implicit none
real, intent(in) :: invL_in,logz0,sig_z_00,sig_y_00,sigy_0_subgid_width_scale,subgrid_delta(2),delta_wind,x
real, intent(out) :: sig_z,sig_y,sig_z_0,sig_y_0
integer i_bot,i_top,i
real weight
real ay,by,az,bz
real min_xy
real ay_pg(6),by_pg(6),az_pg(6),bz_pg(6)
real a_invL(6),b_invL(6),invL(6)
!Use ASME parameters
data ay_pg /0.40,0.36,0.34,0.32,0.315,0.31/
data by_pg /0.91,0.86,0.82,0.78,0.75,0.71/
data az_pg /0.40,0.33,0.27,0.22,0.14,0.06/
data bz_pg /0.91,0.86,0.82,0.78,0.75,0.71/
!Conversion of classes to L
data a_invL /-0.096,-0.037,-0.002,0.0,0.004,0.035/
data b_invL /0.029,0.029,0.018,0.0,-0.018,-0.036/
invL=a_invL+b_invL*logz0
min_xy=(subgrid_delta(1)+subgrid_delta(2))/4.
!Find and interpolate the stability class based on input invL
if (invL_in.le.invL(1)) then
i_bot=1
i_top=1
weight=0.
elseif (invL_in.gt.invL(6)) then
i_bot=6
i_top=6
weight=1.
else
do i=1,5
if (invL_in.gt.invL(i).and.invL_in.le.invL(i+1)) then
i_bot=i
i_top=i+1
weight=(invL_in-invL(i))/(invL(i+1)-invL(i))
exit
endif
enddo
endif
ay=ay_pg(i_bot)*(1.-weight)+ay_pg(i_top)*weight
by=by_pg(i_bot)*(1.-weight)+by_pg(i_top)*weight
az=az_pg(i_bot)*(1.-weight)+az_pg(i_top)*weight
bz=bz_pg(i_bot)*(1.-weight)+bz_pg(i_top)*weight
!Set sig_y_0 to be half of the average x,y grid size
!Add this here ay*exp(by*log(min_xy)) to be the same as sig_z and the same as the Kz calculation
!Does not mean it is correct, just closer to the Kz which is perhaps not so correct
sig_y_0=sig_y_00+min_xy*sigy_0_subgid_width_scale+ay*exp(by*log(min_xy))
!Set sig_z_0 to be the size of the plume after travelling half of the grid size
sig_z_0=sig_z_00+az*exp(bz*log(min_xy))
!Set sig_y and sig_z = sig_0 + a*x^b +x*delata_wind
!Changed this because of error when x=0, check the other similar routines
if (x.le.0) then
sig_y=sig_y_0
sig_z=sig_z_0
else
sig_y=sig_y_0+ay*exp(by*log(x))+x*abs(delta_wind)
sig_z=sig_z_0+az*exp(bz*log(x))
endif
!write(*,'(i,6f)') i,weight,sig_y,sig_z,az,bz,x
end subroutine uEMEP_set_dispersion_sigma_PG