From 355f5ff66081d581191064ca46d8b7c67de7db9e Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Tue, 13 May 2025 15:47:52 -0700 Subject: [PATCH 1/9] Added alpha tuning parameters, one for stable one for unstable regime, to tune the surface layer parameterization --- GFS_layer/GFS_physics_driver.F90 | 5 +- GFS_layer/GFS_typedefs.F90 | 8 +++ gsmphys/sfc_diff_coupled.f | 15 ++++-- gsmphys/sfc_diff_gfdl.f | 84 +++++++++++++++++++++++++------- 4 files changed, 89 insertions(+), 23 deletions(-) diff --git a/GFS_layer/GFS_physics_driver.F90 b/GFS_layer/GFS_physics_driver.F90 index 90ead827..a34f1daa 100644 --- a/GFS_layer/GFS_physics_driver.F90 +++ b/GFS_layer/GFS_physics_driver.F90 @@ -1202,7 +1202,8 @@ subroutine GFS_physics_driver & Sfcprop%ffmm, Sfcprop%ffhh, Sfcprop%uustar, & wind, Tbd%phy_f2d(1,Model%num_p2d), fm10, fh2, & sigmaf, vegtype, Sfcprop%shdmax, Model%ivegsrc, & - tsurf, flag_iter) !, Model%redrag, Model%z0s_max, & + tsurf, flag_iter, Model%alpha_stable, Model%alpha_unstable) +! , Model%redrag, Model%z0s_max, & !Model%do_z0_moon, Model%do_z0_hwrf15, & !Model%do_z0_hwrf17, Model%do_z0_hwrf17_hwonly, & !Model%wind_th_hwrf) @@ -1225,7 +1226,7 @@ subroutine GFS_physics_driver & tsurf, flag_iter, Model%redrag, Model%z0s_max, & Model%do_z0_moon, Model%do_z0_hwrf15, & Model%do_z0_hwrf17, Model%do_z0_hwrf17_hwonly, & - Model%wind_th_hwrf) + Model%wind_th_hwrf, Model%alpha_stable, Model%alpha_unstable) else ! GFS original sfc_diff modified by kgao call sfc_diff (im,Statein%pgr, Statein%ugrs, Statein%vgrs,& diff --git a/GFS_layer/GFS_typedefs.F90 b/GFS_layer/GFS_typedefs.F90 index c75ed8b5..b41e1b13 100644 --- a/GFS_layer/GFS_typedefs.F90 +++ b/GFS_layer/GFS_typedefs.F90 @@ -680,6 +680,8 @@ module GFS_typedefs logical :: do_z0_hwrf17 !< flag for using z0 scheme in 2017 HWRF (kgao) logical :: do_z0_hwrf17_hwonly !< flag for using z0 scheme in 2017 HWRF only under high wind (kgao) real(kind=kind_phys) :: wind_th_hwrf !< wind speed threshold when z0 level off as in HWRF (kgao) + real(kind=kind_phys) :: alpha_stable !< description TBD + real(kind=kind_phys) :: alpha_unstable !< description TBD logical :: hybedmf !< flag for hybrid edmf pbl scheme logical :: myj_pbl !< flag for NAM MYJ tke scheme logical :: ysupbl !< flag for ysu pbl scheme (version in WRFV3.8) @@ -2401,6 +2403,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & logical :: do_z0_hwrf17 = .false. !< flag for using z0 scheme in 2017 HWRF logical :: do_z0_hwrf17_hwonly = .false. !< flag for using z0 scheme in 2017 HWRF only under high wind real(kind=kind_phys) :: wind_th_hwrf = 33. !< wind speed threshold when z0 level off as in HWRF + real(kind=kind_phys) :: alpha_stable = 5. !< TBD + real(kind=kind_phys) :: alpha_unstable = 16. !< TBD logical :: hybedmf = .false. !< flag for hybrid edmf pbl scheme logical :: myj_pbl = .false. !< flag for NAM MYJ tke-based scheme logical :: ysupbl = .false. !< flag for hybrid edmf pbl scheme @@ -2878,6 +2882,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%do_z0_hwrf17 = do_z0_hwrf17 Model%do_z0_hwrf17_hwonly = do_z0_hwrf17_hwonly Model%wind_th_hwrf = wind_th_hwrf + Model%alpha_stable = alpha_stable + Model%alpha_unstable = alpha_unstable Model%hybedmf = hybedmf Model%myj_pbl = myj_pbl Model%ysupbl = ysupbl @@ -3585,6 +3591,8 @@ subroutine control_print(Model) print *, ' do_z0_hwrf17 : ', Model%do_z0_hwrf17 print *, ' do_z0_hwrf17_hwonly : ', Model%do_z0_hwrf17_hwonly print *, ' wind_th_hwrf : ', Model%wind_th_hwrf + print *, ' alpha_stable : ', Model%alpha_stable + print *, ' alpha_unstable : ', Model%alpha_unstable print *, ' hybedmf : ', Model%hybedmf print *, ' myj_pbl : ', Model%myj_pbl print *, ' ysupbl : ', Model%ysupbl diff --git a/gsmphys/sfc_diff_coupled.f b/gsmphys/sfc_diff_coupled.f index dfe08bcb..03b6d372 100644 --- a/gsmphys/sfc_diff_coupled.f +++ b/gsmphys/sfc_diff_coupled.f @@ -4,7 +4,8 @@ subroutine sfc_diff_coupled(im,ps,u1,v1,t1,q1,z1, & stress,fm,fh, & ustar,wind,ddvel,fm10,fh2, & sigmaf,vegtype,shdmax,ivegsrc, - & tsurf,flag_iter) !,redrag, + & tsurf,flag_iter, alpha_stable, alpha_unstable) +! ,redrag, ! & z0s_max) ! & do_z0_moon, do_z0_hwrf15, do_z0_hwrf17, ! & do_z0_hwrf17_hwonly, wind_th_hwrf) @@ -26,7 +27,9 @@ subroutine sfc_diff_coupled(im,ps,u1,v1,t1,q1,z1, &, prsl1, prslki, stress &, fm, fh, ustar, wind, ddvel &, fm10, fh2, sigmaf, shdmax - &, tsurf, snwdph + &, tsurf, snwdph, alpha_stable + &, alpha_unstable + integer, dimension(im) ::vegtype, islimsk logical flag_iter(im) @@ -143,7 +146,9 @@ subroutine sfc_diff_coupled(im,ps,u1,v1,t1,q1,z1, ! --- call similarity call monin_obukhov_similarity - & (z1(i), snwdph(i), thv1, wind(i), z0max, ztmax, tvs, + & (islimsk(i), z1(i), snwdph(i), thv1, wind(i), z0max, + & ztmax, tvs, + & alpha_stable, alpha_unstable, & rb(i), fm(i), fh(i), fm10(i), fh2(i), & cm(i), ch(i), stress(i), ustar(i)) @@ -165,7 +170,9 @@ subroutine sfc_diff_coupled(im,ps,u1,v1,t1,q1,z1, ! --- call similarity ! kgao: use z0/zt to do sfc diagnosis call monin_obukhov_similarity - & (z1(i), snwdph(i), thv1, wind(i), z0max, ztmax, tvs, + & (islimsk(i), z1(i), snwdph(i), thv1, wind(i), z0max, + & ztmax, tvs, + & alpha_stable, alpha_unstable, & rb(i), fm(i), fh(i), fm10(i), fh2(i), & cm(i), ch(i), tem1, tem2) !stress(i), ustar(i)) diff --git a/gsmphys/sfc_diff_gfdl.f b/gsmphys/sfc_diff_gfdl.f index 87aeaf7f..5ef69f64 100644 --- a/gsmphys/sfc_diff_gfdl.f +++ b/gsmphys/sfc_diff_gfdl.f @@ -22,7 +22,8 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, & tsurf,flag_iter,redrag, & z0s_max, & do_z0_moon, do_z0_hwrf15, do_z0_hwrf17, - & do_z0_hwrf17_hwonly, wind_th_hwrf) + & do_z0_hwrf17_hwonly, wind_th_hwrf, + & alpha_stable, alpha_unstable) ! oct 2019 - a clean and updated version by Kun Gao at GFDL (Kun.Gao@noaa.gov) @@ -50,7 +51,7 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, &, fm10, fh2, sigmaf, shdmax &, tsurf, snwdph &, fm_neutral, fm10_neutral ! Sofar added Spring 2023 - real(kind=kind_phys) :: ws1, ws10n ! Sofar added 10/19/23 + real(kind=kind_phys) :: ws1, ws10n, alpha_stable, alpha_unstable ! Sofar added 10/19/23 integer, dimension(im) ::vegtype, islimsk logical flag_iter(im) @@ -171,7 +172,9 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, ! --- call similarity call monin_obukhov_similarity - & (z1(i), snwdph(i), thv1, wind(i), z0max, ztmax, tvs, + & (islimsk(i), z1(i), snwdph(i), thv1, wind(i), z0max, + & ztmax, tvs, + & alpha_stable, alpha_unstable, & rb(i), fm(i), fh(i), fm10(i), fh2(i), & fm_neutral(i), fm10_neutral(i), !(ADDED by Sofar) & cm(i), ch(i), stress(i), ustar(i)) @@ -202,7 +205,9 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, ! --- call similarity call monin_obukhov_similarity - & (z1(i), snwdph(i), thv1, wind(i), z0max, ztmax, tvs, + & (islimsk(i), z1(i), snwdph(i), thv1, wind(i), z0max, + & ztmax, tvs, + & alpha_stable, alpha_unstable, & rb(i), fm(i), fh(i), fm10(i), fh2(i), & fm_neutral(i), fm10_neutral(i), !(ADDED by Sofar) & cm(i), ch(i), stress(i), ustar(i)) @@ -302,7 +307,9 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, ! --- call similarity call monin_obukhov_similarity - & (z1(i), snwdph(i), thv1, wind(i), z0max, ztmax, tvs, + & (islimsk(i), z1(i), snwdph(i), thv1, wind(i), z0max, + & ztmax, tvs, + & alpha_stable, alpha_unstable, & rb(i), fm(i), fh(i), fm10(i), fh2(i), & fm_neutral(i), fm10_neutral(i), !(ADDED by Sofar) & cm(i), ch(i), stress(i), ustar(i)) @@ -508,8 +515,9 @@ end subroutine cal_z0_moon ! ======================================================================= subroutine monin_obukhov_similarity - & ( z1, snwdph, thv1, wind, z0max, ztmax, tvs, - & rb, fm, fh, fm10, fh2, + & ( ilsimask, z1, snwdph, thv1, wind, z0max, ztmax, tvs, + & alpha_stable, alpha_unstable, + & rb, fm, fh, fm10, fh2, & fm_neutral, fm10_neutral, !(ADDED by Sofar) & cm, ch, stress, ustar) @@ -536,7 +544,10 @@ subroutine monin_obukhov_similarity ! --- inputs: real(kind=kind_phys), intent(in) :: - & z1, snwdph, thv1, wind, z0max, ztmax, tvs + & z1, snwdph, thv1, wind, z0max, ztmax, tvs, alpha_stable, + & alpha_unstable + + integer, intent(in) :: ilsimask ! --- outputs: real(kind=kind_phys), intent(out) :: @@ -547,7 +558,8 @@ subroutine monin_obukhov_similarity ! --- locals: real(kind=kind_phys), parameter :: alpha=5., a0=-3.975 - &, a1=12.32, alpha4=4.0*alpha + &, a1=12.32, alpha4_nonocean=4.0*alpha + &, alpha4_ocean = 4.*alpha_stable &, b1=-7.755, b2=6.041, alpha2=alpha+alpha, beta=1.0 &, a0p=-7.941, a1p=24.75, b1p=-8.705, b2p=7.899 &, ztmin1=-999.0, ca=.4 @@ -557,10 +569,16 @@ subroutine monin_obukhov_similarity & z1i, & fms, fhs, hl0, hl0inf, hlinf, & hl110, hlt, hltinf, olinf, - & tem1, tem2, ztmax1 + & tem1, tem2, ztmax1, alpha4 z1i = 1.0 / z1 + if (ilsimask == 0) then + alpha4 = alpha4_ocean + else + alpha4 = alpha4_nonocean + end if + tem1 = z0max/z1 if (abs(1.0-tem1) > 1.0e-6) then ztmax1 = - beta*log(tem1)/(alpha2*(1.-tem1)) @@ -642,32 +660,49 @@ subroutine monin_obukhov_similarity ! ! get pm and ph ! + if (hlinf >= -0.5) then hl1 = hlinf - pm = (a0 + a1*hl1) * hl1 / (1.+ (b1+b2*hl1) *hl1) - ph = (a0p + a1p*hl1) * hl1 / (1.+ (b1p+b2p*hl1)*hl1) hl110 = hl1 * 10. * z1i hl110 = min(max(hl110, ztmin1), ztmax1) - pm10 = (a0 + a1*hl110) * hl110 / (1.+(b1+b2*hl110)*hl110) + if (ilsimask == 0) then + call get_psim_unstable(hlinf, z1i, z0max, + & alpha_unstable, pm) + call get_psim_unstable(hl110, 1./10., z0max, + & alpha_unstable, pm10) + else + pm = (a0 + a1*hl1) * hl1 / (1.+ (b1+b2*hl1) *hl1) + pm10 = (a0 + a1*hl110) * hl110 / + & (1.+(b1+b2*hl110)*hl110) + endif + ph = (a0p + a1p*hl1) * hl1 / (1.+ (b1p+b2p*hl1)*hl1) hl12 = (hl1+hl1) * z1i hl12 = min(max(hl12, ztmin1), ztmax1) ph2 = (a0p + a1p*hl12) * hl12 / (1.+(b1p+b2p*hl12)*hl12) else ! hlinf < 0.05 hl1 = -hlinf + hl110 = hl1 * 10. * z1i + hl110 = min(max(hl110, ztmin1), ztmax1) tem1 = 1.0 / sqrt(hl1) - pm = log(hl1) + 2. * sqrt(tem1) - .8776 + if (ilsimask == 0) then + call get_psim_unstable(hlinf, z1i, z0max, + & alpha_unstable, pm) + call get_psim_unstable(-hl110, 1./10., z0max, + & alpha_unstable, pm10) + else + pm = log(hl1) + 2. * sqrt(tem1) - .8776 + pm10 = log(hl110) + 2.0 / sqrt(sqrt(hl110)) - .8776 + endif ph = log(hl1) + .5 * tem1 + 1.386 ! pm = log(hl1) + 2.0 * hl1 ** (-.25) - .8776 ! ph = log(hl1) + 0.5 * hl1 ** (-.5) + 1.386 - hl110 = hl1 * 10. * z1i - hl110 = min(max(hl110, ztmin1), ztmax1) - pm10 = log(hl110) + 2.0 / sqrt(sqrt(hl110)) - .8776 ! pm10 = log(hl110) + 2. * hl110 ** (-.25) - .8776 hl12 = (hl1+hl1) * z1i hl12 = min(max(hl12, ztmin1), ztmax1) ph2 = log(hl12) + 0.5 / sqrt(hl12) + 1.386 ! ph2 = log(hl12) + .5 * hl12 ** (-.5) + 1.386 endif + endif ! end of if (dtv >= 0 ) then loop ! @@ -687,3 +722,18 @@ subroutine monin_obukhov_similarity return end subroutine monin_obukhov_similarity + + subroutine get_psim_unstable(zeta, z1i, z0, alpha, psi_m) + use machine , only : kind_phys + implicit none + real(kind=kind_phys),intent(in) :: zeta, z1i, z0, alpha + real(kind=kind_phys),intent(out) :: psi_m + real(kind=kind_phys) :: tem1, zeta0, Rz, R0 + tem1 = zeta * z1i + zeta0 = z0 * tem1 + Rz = (1. - alpha * zeta)**0.25 + R0 = (1. - alpha * zeta0)**0.25 + psi_m = log((1. + Rz)**2. * (1. + Rz**2.) / ((1. + R0)**2. + & *(1. + R0**2.))) + 2. * np.arctan(R0) - 2. * np.arctan(Rz) + return + end subroutine From 21a177db6a2cd3ae7b78c534c721a8741ea7e981 Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Tue, 13 May 2025 17:00:23 -0700 Subject: [PATCH 2/9] Added parameter tune_ocean_surface_layer to switch on/off the tuning of the SL layer over oceans --- GFS_layer/GFS_typedefs.F90 | 7 ++++++- gsmphys/sfc_diff_coupled.f | 9 +++++---- gsmphys/sfc_diff_gfdl.f | 20 +++++++++++--------- 3 files changed, 22 insertions(+), 14 deletions(-) diff --git a/GFS_layer/GFS_typedefs.F90 b/GFS_layer/GFS_typedefs.F90 index b41e1b13..1170dfc3 100644 --- a/GFS_layer/GFS_typedefs.F90 +++ b/GFS_layer/GFS_typedefs.F90 @@ -682,6 +682,7 @@ module GFS_typedefs real(kind=kind_phys) :: wind_th_hwrf !< wind speed threshold when z0 level off as in HWRF (kgao) real(kind=kind_phys) :: alpha_stable !< description TBD real(kind=kind_phys) :: alpha_unstable !< description TBD + logical :: tune_ocean_surface_layer !< if true, use alphas above to modify the SL similarity profiles over oceans logical :: hybedmf !< flag for hybrid edmf pbl scheme logical :: myj_pbl !< flag for NAM MYJ tke scheme logical :: ysupbl !< flag for ysu pbl scheme (version in WRFV3.8) @@ -2405,6 +2406,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & real(kind=kind_phys) :: wind_th_hwrf = 33. !< wind speed threshold when z0 level off as in HWRF real(kind=kind_phys) :: alpha_stable = 5. !< TBD real(kind=kind_phys) :: alpha_unstable = 16. !< TBD + logical :: tune_ocean_surface_layer = .false. !< if true, use alphas above to modify the SL similarity profiles over oceans logical :: hybedmf = .false. !< flag for hybrid edmf pbl scheme logical :: myj_pbl = .false. !< flag for NAM MYJ tke-based scheme logical :: ysupbl = .false. !< flag for hybrid edmf pbl scheme @@ -2637,7 +2639,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & cscnv, cal_pre, do_aw, do_shoc, shocaftcnv, shoc_cld, & h2o_phys, pdfcld, shcnvcw, redrag, sfc_gfdl, z0s_max, & sfc_coupled, do_z0_moon, do_z0_hwrf15, do_z0_hwrf17, & - do_z0_hwrf17_hwonly, wind_th_hwrf, & + do_z0_hwrf17_hwonly, wind_th_hwrf, alpha_stable, & + alpha_unstable, tune_ocean_surface_layer, & & hybedmf, dspheat, lheatstrg, hour_canopy, afac_canopy, & cnvcld, no_pbl, xkzm_lim, xkzm_fac, xkgdx, & rlmn, rlmx, zolcru, cs0, & @@ -2884,6 +2887,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%wind_th_hwrf = wind_th_hwrf Model%alpha_stable = alpha_stable Model%alpha_unstable = alpha_unstable + Model%tune_ocean_surface_layer = tune_ocean_surface_layer Model%hybedmf = hybedmf Model%myj_pbl = myj_pbl Model%ysupbl = ysupbl @@ -3593,6 +3597,7 @@ subroutine control_print(Model) print *, ' wind_th_hwrf : ', Model%wind_th_hwrf print *, ' alpha_stable : ', Model%alpha_stable print *, ' alpha_unstable : ', Model%alpha_unstable + print *, ' tune_ocean_surface_layer : ', Model%tune_ocean_surface_layer print *, ' hybedmf : ', Model%hybedmf print *, ' myj_pbl : ', Model%myj_pbl print *, ' ysupbl : ', Model%ysupbl diff --git a/gsmphys/sfc_diff_coupled.f b/gsmphys/sfc_diff_coupled.f index 03b6d372..e8eb2d17 100644 --- a/gsmphys/sfc_diff_coupled.f +++ b/gsmphys/sfc_diff_coupled.f @@ -4,7 +4,8 @@ subroutine sfc_diff_coupled(im,ps,u1,v1,t1,q1,z1, & stress,fm,fh, & ustar,wind,ddvel,fm10,fh2, & sigmaf,vegtype,shdmax,ivegsrc, - & tsurf,flag_iter, alpha_stable, alpha_unstable) + & tsurf,flag_iter, alpha_stable, alpha_unstable, + & tune_ocean_surface_layer) ! ,redrag, ! & z0s_max) ! & do_z0_moon, do_z0_hwrf15, do_z0_hwrf17, @@ -32,7 +33,7 @@ subroutine sfc_diff_coupled(im,ps,u1,v1,t1,q1,z1, integer, dimension(im) ::vegtype, islimsk - logical flag_iter(im) + logical flag_iter(im), tune_ocean_surface_layer ! logical redrag ! logical do_z0_moon, do_z0_hwrf15, do_z0_hwrf17 ! kgao ! &, do_z0_hwrf17_hwonly ! kgao @@ -148,7 +149,7 @@ subroutine sfc_diff_coupled(im,ps,u1,v1,t1,q1,z1, call monin_obukhov_similarity & (islimsk(i), z1(i), snwdph(i), thv1, wind(i), z0max, & ztmax, tvs, - & alpha_stable, alpha_unstable, + & alpha_stable, alpha_unstable, tune_ocean_surface_layer, & rb(i), fm(i), fh(i), fm10(i), fh2(i), & cm(i), ch(i), stress(i), ustar(i)) @@ -172,7 +173,7 @@ subroutine sfc_diff_coupled(im,ps,u1,v1,t1,q1,z1, call monin_obukhov_similarity & (islimsk(i), z1(i), snwdph(i), thv1, wind(i), z0max, & ztmax, tvs, - & alpha_stable, alpha_unstable, + & alpha_stable, alpha_unstable, tune_ocean_surface_layer, & rb(i), fm(i), fh(i), fm10(i), fh2(i), & cm(i), ch(i), tem1, tem2) !stress(i), ustar(i)) diff --git a/gsmphys/sfc_diff_gfdl.f b/gsmphys/sfc_diff_gfdl.f index 5ef69f64..95b8d2fb 100644 --- a/gsmphys/sfc_diff_gfdl.f +++ b/gsmphys/sfc_diff_gfdl.f @@ -23,7 +23,8 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, & z0s_max, & do_z0_moon, do_z0_hwrf15, do_z0_hwrf17, & do_z0_hwrf17_hwonly, wind_th_hwrf, - & alpha_stable, alpha_unstable) + & alpha_stable, alpha_unstable, + & tune_ocean_surface_layer) ! oct 2019 - a clean and updated version by Kun Gao at GFDL (Kun.Gao@noaa.gov) @@ -54,7 +55,7 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, real(kind=kind_phys) :: ws1, ws10n, alpha_stable, alpha_unstable ! Sofar added 10/19/23 integer, dimension(im) ::vegtype, islimsk - logical flag_iter(im) + logical flag_iter(im), tune_ocean_surface_layer logical redrag logical do_z0_moon, do_z0_hwrf15, do_z0_hwrf17 ! kgao &, do_z0_hwrf17_hwonly ! kgao @@ -174,7 +175,7 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, call monin_obukhov_similarity & (islimsk(i), z1(i), snwdph(i), thv1, wind(i), z0max, & ztmax, tvs, - & alpha_stable, alpha_unstable, + & alpha_stable, alpha_unstable, tune_ocean_surface_layer, & rb(i), fm(i), fh(i), fm10(i), fh2(i), & fm_neutral(i), fm10_neutral(i), !(ADDED by Sofar) & cm(i), ch(i), stress(i), ustar(i)) @@ -207,7 +208,7 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, call monin_obukhov_similarity & (islimsk(i), z1(i), snwdph(i), thv1, wind(i), z0max, & ztmax, tvs, - & alpha_stable, alpha_unstable, + & alpha_stable, alpha_unstable, tune_ocean_surface_layer, & rb(i), fm(i), fh(i), fm10(i), fh2(i), & fm_neutral(i), fm10_neutral(i), !(ADDED by Sofar) & cm(i), ch(i), stress(i), ustar(i)) @@ -309,7 +310,7 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, call monin_obukhov_similarity & (islimsk(i), z1(i), snwdph(i), thv1, wind(i), z0max, & ztmax, tvs, - & alpha_stable, alpha_unstable, + & alpha_stable, alpha_unstable, tune_ocean_surface_layer, & rb(i), fm(i), fh(i), fm10(i), fh2(i), & fm_neutral(i), fm10_neutral(i), !(ADDED by Sofar) & cm(i), ch(i), stress(i), ustar(i)) @@ -516,7 +517,7 @@ end subroutine cal_z0_moon subroutine monin_obukhov_similarity & ( ilsimask, z1, snwdph, thv1, wind, z0max, ztmax, tvs, - & alpha_stable, alpha_unstable, + & alpha_stable, alpha_unstable, tune_ocean_surface_layer, & rb, fm, fh, fm10, fh2, & fm_neutral, fm10_neutral, !(ADDED by Sofar) & cm, ch, stress, ustar) @@ -548,6 +549,7 @@ subroutine monin_obukhov_similarity & alpha_unstable integer, intent(in) :: ilsimask + logical, intent(in) :: tune_ocean_surface_layer ! --- outputs: real(kind=kind_phys), intent(out) :: @@ -573,7 +575,7 @@ subroutine monin_obukhov_similarity z1i = 1.0 / z1 - if (ilsimask == 0) then + if (ilsimask == 0 .and. tune_ocean_surface_layer) then alpha4 = alpha4_ocean else alpha4 = alpha4_nonocean @@ -665,7 +667,7 @@ subroutine monin_obukhov_similarity hl1 = hlinf hl110 = hl1 * 10. * z1i hl110 = min(max(hl110, ztmin1), ztmax1) - if (ilsimask == 0) then + if (ilsimask == 0 .and. tune_ocean_surface_layer) then call get_psim_unstable(hlinf, z1i, z0max, & alpha_unstable, pm) call get_psim_unstable(hl110, 1./10., z0max, @@ -684,7 +686,7 @@ subroutine monin_obukhov_similarity hl110 = hl1 * 10. * z1i hl110 = min(max(hl110, ztmin1), ztmax1) tem1 = 1.0 / sqrt(hl1) - if (ilsimask == 0) then + if (ilsimask == 0 .and. tune_ocean_surface_layer) then call get_psim_unstable(hlinf, z1i, z0max, & alpha_unstable, pm) call get_psim_unstable(-hl110, 1./10., z0max, From ccd9ba3608582addeb48e8f2771f90228312e18e Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Wed, 14 May 2025 09:41:04 -0700 Subject: [PATCH 3/9] Add missing description of new parameters --- GFS_layer/GFS_typedefs.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/GFS_layer/GFS_typedefs.F90 b/GFS_layer/GFS_typedefs.F90 index 1170dfc3..c2f53b74 100644 --- a/GFS_layer/GFS_typedefs.F90 +++ b/GFS_layer/GFS_typedefs.F90 @@ -680,9 +680,9 @@ module GFS_typedefs logical :: do_z0_hwrf17 !< flag for using z0 scheme in 2017 HWRF (kgao) logical :: do_z0_hwrf17_hwonly !< flag for using z0 scheme in 2017 HWRF only under high wind (kgao) real(kind=kind_phys) :: wind_th_hwrf !< wind speed threshold when z0 level off as in HWRF (kgao) - real(kind=kind_phys) :: alpha_stable !< description TBD - real(kind=kind_phys) :: alpha_unstable !< description TBD - logical :: tune_ocean_surface_layer !< if true, use alphas above to modify the SL similarity profiles over oceans + real(kind=kind_phys) :: alpha_stable !< tuning parameter for the dimensionless momentum and scalar gradient function in the surface layer + real(kind=kind_phys) :: alpha_unstable !< tuning parameter for the dimensionless momentum gradient function in the surface layer + logical :: tune_ocean_surface_layer !< if true, use alphas above to modify the surface layer dimensionless gradients over oceans logical :: hybedmf !< flag for hybrid edmf pbl scheme logical :: myj_pbl !< flag for NAM MYJ tke scheme logical :: ysupbl !< flag for ysu pbl scheme (version in WRFV3.8) @@ -2404,9 +2404,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & logical :: do_z0_hwrf17 = .false. !< flag for using z0 scheme in 2017 HWRF logical :: do_z0_hwrf17_hwonly = .false. !< flag for using z0 scheme in 2017 HWRF only under high wind real(kind=kind_phys) :: wind_th_hwrf = 33. !< wind speed threshold when z0 level off as in HWRF - real(kind=kind_phys) :: alpha_stable = 5. !< TBD - real(kind=kind_phys) :: alpha_unstable = 16. !< TBD - logical :: tune_ocean_surface_layer = .false. !< if true, use alphas above to modify the SL similarity profiles over oceans + real(kind=kind_phys) :: alpha_stable = 5. !< tuning parameter for the dimensionless momentum and scalar gradient function in the surface layer + real(kind=kind_phys) :: alpha_unstable = 16. !< tuning parameter for the dimensionless momentum gradient function in the surface layer + logical :: tune_ocean_surface_layer = .false. !< if true, use alphas above to modify the surface layer dimensionless gradients over oceans logical :: hybedmf = .false. !< flag for hybrid edmf pbl scheme logical :: myj_pbl = .false. !< flag for NAM MYJ tke-based scheme logical :: ysupbl = .false. !< flag for hybrid edmf pbl scheme From 2fb14099661322e769274832a0acaab2a426b88a Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Wed, 14 May 2025 09:57:39 -0700 Subject: [PATCH 4/9] Added what was missing because of auto-save deactivated --- GFS_layer/GFS_physics_driver.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/GFS_layer/GFS_physics_driver.F90 b/GFS_layer/GFS_physics_driver.F90 index a34f1daa..0e266cb3 100644 --- a/GFS_layer/GFS_physics_driver.F90 +++ b/GFS_layer/GFS_physics_driver.F90 @@ -1202,7 +1202,8 @@ subroutine GFS_physics_driver & Sfcprop%ffmm, Sfcprop%ffhh, Sfcprop%uustar, & wind, Tbd%phy_f2d(1,Model%num_p2d), fm10, fh2, & sigmaf, vegtype, Sfcprop%shdmax, Model%ivegsrc, & - tsurf, flag_iter, Model%alpha_stable, Model%alpha_unstable) + tsurf, flag_iter, & + Model%alpha_stable, Model%alpha_unstable, Model%tune_ocean_surface_layer) ! , Model%redrag, Model%z0s_max, & !Model%do_z0_moon, Model%do_z0_hwrf15, & !Model%do_z0_hwrf17, Model%do_z0_hwrf17_hwonly, & @@ -1226,7 +1227,8 @@ subroutine GFS_physics_driver & tsurf, flag_iter, Model%redrag, Model%z0s_max, & Model%do_z0_moon, Model%do_z0_hwrf15, & Model%do_z0_hwrf17, Model%do_z0_hwrf17_hwonly, & - Model%wind_th_hwrf, Model%alpha_stable, Model%alpha_unstable) + Model%wind_th_hwrf, & + Model%alpha_stable, Model%alpha_unstable, Model%tune_ocean_surface_layer) else ! GFS original sfc_diff modified by kgao call sfc_diff (im,Statein%pgr, Statein%ugrs, Statein%vgrs,& From 34c9d183e1f3b1b88e28ab5aa3f95f139e127c02 Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Wed, 14 May 2025 10:06:31 -0700 Subject: [PATCH 5/9] Added description of new input arguments to MO subroutine --- gsmphys/sfc_diff_gfdl.f | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/gsmphys/sfc_diff_gfdl.f b/gsmphys/sfc_diff_gfdl.f index 95b8d2fb..305668d9 100644 --- a/gsmphys/sfc_diff_gfdl.f +++ b/gsmphys/sfc_diff_gfdl.f @@ -523,11 +523,16 @@ subroutine monin_obukhov_similarity & cm, ch, stress, ustar) ! --- input +! ilsimask - land/sea/ice mask ! z1 - lowest model level height ! snwdph - surface snow thickness ! wind - wind speed at lowest model layer ! thv1 - virtual potential temp at lowest model layer ! tvs - surface temp +! alpha_stable - tuning parameter (inverse crit. Richarson number) in dimensinless shear +! and scalar gradient; stable case +! alpha_unstable - tuning parameter in dimensionless wind shear; unstable case +! tune_ocean_surface_layer - if true, use alphas to tune surface layer over oceans ! z0max - surface roughness length for momentum ! ztmax - surface roughness length for heat ! From 4571c80edd5d3c42a09a1a1ba32e0f268d3544fd Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Wed, 14 May 2025 15:31:31 -0700 Subject: [PATCH 6/9] Make the PBL schemes (not all of them) consistent with the implemented tuning by using the modified similarity functions to define phim and phih in satmedmfvdiff, satmedmfvdiffq, and moninedmf; a subroutine to compute the similarity functions has been created to avoid redundancy --- GFS_layer/GFS_physics_driver.F90 | 10 ++-- gsmphys/moninedmf.f | 30 ++++------- gsmphys/satmedmfvdiff.f | 88 ++++++++++++++++++++++++-------- gsmphys/satmedmfvdifq.f | 30 ++++------- 4 files changed, 93 insertions(+), 65 deletions(-) diff --git a/GFS_layer/GFS_physics_driver.F90 b/GFS_layer/GFS_physics_driver.F90 index 0e266cb3..b91a22c0 100644 --- a/GFS_layer/GFS_physics_driver.F90 +++ b/GFS_layer/GFS_physics_driver.F90 @@ -1712,7 +1712,8 @@ subroutine GFS_physics_driver & Model%dspheat, dusfc1, dvsfc1, dtsfc1, dqsfc1, Diag%hpbl,& gamt, gamq, dkt, kinver, Model%xkzm_m, Model%xkzm_h, & Model%xkzm_s, lprnt, ipr, & - Model%xkzminv, Model%moninq_fac) + Model%xkzminv, Model%moninq_fac, Model%alpha_stable, & + Model%alpha_unstable, Model%tune_ocean_surface_layer) ! if (lprnt) write(0,*)' dtdtm=',(dtdt(ipr,k),k=1,15) ! if (lprnt) write(0,*)' dqdtm=',(dqdt(ipr,k,1),k=1,15) @@ -1734,7 +1735,9 @@ subroutine GFS_physics_driver & Model%xkzm_ml, Model%xkzm_hl, Model%xkzm_mi, Model%xkzm_hi, & Model%xkzm_s, Model%xkzminv, Model%do_dk_hb19, & Model%xkzm_lim, Model%xkgdx, & - Model%rlmn, Model%rlmx, Model%cap_k0_land, dkt) + Model%rlmn, Model%rlmx, Model%cap_k0_land, dkt, & + Model%alpha_stable, Model%alpha_unstable, & + Model%tune_ocean_surface_layer) elseif (Model%isatmedmf == 1) then do i=1,im @@ -1761,7 +1764,8 @@ subroutine GFS_physics_driver & Model%xkzm_s, Model%xkzminv, Model%rlmx, Model%zolcru, & Model%cs0, Model%do_dk_hb19, Model%xkgdx, & Model%dspfac, Model%bl_upfr, Model%bl_dnfr, dkt, & - flux_cg, flux_en) !cg as up and en as down + flux_cg, flux_en, Model%alpha_stable, Model%alpha_unstable, & + Model%tune_ocean_surface_layer) !cg as up and en as down endif elseif (Model%ysupbl) then diff --git a/gsmphys/moninedmf.f b/gsmphys/moninedmf.f index 1c42c6d7..28e010b5 100755 --- a/gsmphys/moninedmf.f +++ b/gsmphys/moninedmf.f @@ -93,7 +93,8 @@ subroutine moninedmf(ix,im,km,ntrac,ntcw,dv,du,tau,rtg, & & prsi,del,prsl,prslk,phii,phil,delt,dspheat, & & dusfc,dvsfc,dtsfc,dqsfc,hpbl,hgamt,hgamq,dkt, & & kinver,xkzm_m,xkzm_h,xkzm_s,lprnt,ipr, & - & xkzminv,moninq_fac) + & xkzminv,moninq_fac, alpha_stable, alpha_unstable, & + & tune_ocean_surface_layer) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -103,11 +104,12 @@ subroutine moninedmf(ix,im,km,ntrac,ntcw,dv,du,tau,rtg, & ! ! arguments ! - logical lprnt + logical lprnt, tune_ocean_surface_layer integer ipr integer ix, im, km, ntrac, ntcw, kpbl(im), kinver(im) ! - real(kind=kind_phys) delt, xkzm_m, xkzm_h, xkzm_s + real(kind=kind_phys) delt, xkzm_m, xkzm_h, xkzm_s, alpha_stable, + & alpha_unstable real(kind=kind_phys) dv(im,km), du(im,km), & & tau(im,km), rtg(im,km,ntrac), & & u1(ix,km), v1(ix,km), & @@ -512,24 +514,12 @@ subroutine moninedmf(ix,im,km,ntrac,ntcw,dv,du,tau,rtg, & !! \f[ !! w_s = (u_*^3 + 7\epsilon k w_*^3)^{1/3} !! \f] + call get_similarity_parameters(im, rbsoil, fm, fh, rimin, + & sfcflg, zfmin, sfcfrac, hpbl, zl(1:im,1), + & aphi16, aphi5, alpha_unstable, alpha_stable, + & tune_ocean_surface_layer, phim, phih, zol) + do i=1,im - zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin) - if(sfcflg(i)) then - zol(i) = min(zol(i),-zfmin) - else - zol(i) = max(zol(i),zfmin) - endif - zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i,1) - if(sfcflg(i)) then -! phim(i) = (1.-aphi16*zol1)**(-1./4.) -! phih(i) = (1.-aphi16*zol1)**(-1./2.) - tem = 1.0 / (1. - aphi16*zol1) - phih(i) = sqrt(tem) - phim(i) = sqrt(phih(i)) - else - phim(i) = 1. + aphi5*zol1 - phih(i) = phim(i) - endif wscale(i) = ustar(i)/phim(i) ustmin(i) = ustar(i)/aphi5 wscale(i) = max(wscale(i),ustmin(i)) diff --git a/gsmphys/satmedmfvdiff.f b/gsmphys/satmedmfvdiff.f index 4c7e32f4..259d0f43 100644 --- a/gsmphys/satmedmfvdiff.f +++ b/gsmphys/satmedmfvdiff.f @@ -36,7 +36,8 @@ subroutine satmedmfvdif(ix,im,km,ntrac,ntcw,ntiw,ntke, & dspheat,dusfc,dvsfc,dtsfc,dqsfc,hpbl, & kinver,xkzm_mo,xkzm_ho,xkzm_ml,xkzm_hl,xkzm_mi,xkzm_hi, & xkzm_s,xkzinv,do_dk_hb19,xkzm_lim,xkgdx, - & rlmn, rlmx, cap_k0_land, dkt_out) + & rlmn, rlmx, cap_k0_land, dkt_out, + & alpha_stable, alpha_unstable, tune_ocean_surface_layer) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -53,7 +54,8 @@ subroutine satmedmfvdif(ix,im,km,ntrac,ntcw,ntiw,ntke, ! real(kind=kind_phys) delt, xkzm_s, xkzm_lim, & xkzm_mo, xkzm_ho, xkzm_ml, xkzm_hl, - & xkzm_mi, xkzm_hi + & xkzm_mi, xkzm_hi, + & alpha_stable, alpha_unstable real(kind=kind_phys) dv(im,km), du(im,km), & tdt(im,km), rtg(im,km,ntrac), & u1(ix,km), v1(ix,km), @@ -77,7 +79,7 @@ subroutine satmedmfvdif(ix,im,km,ntrac,ntcw,ntiw,ntke, ! kgao note - q1 and rtg are local var now ! - logical dspheat, cap_k0_land, do_dk_hb19 + logical dspheat, cap_k0_land, do_dk_hb19, tune_ocean_surface_layer ! flag for tke dissipative heating real(kind=kind_phys),dimension(1:im,1:km),intent(OUT)::dkt_out @@ -592,24 +594,10 @@ subroutine satmedmfvdif(ix,im,km,ntrac,ntcw,ntiw,ntke, ! ! compute similarity parameters ! - do i=1,im - zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin) - if(sfcflg(i)) then - zol(i) = min(zol(i),-zfmin) - else - zol(i) = max(zol(i),zfmin) - endif -! - zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i,1) - if(sfcflg(i)) then - tem = 1.0 / (1. - aphi16*zol1) - phih(i) = sqrt(tem) - phim(i) = sqrt(phih(i)) - else - phim(i) = 1. + aphi5*zol1 - phih(i) = phim(i) - endif - enddo + call get_similarity_parameters(im, rbsoil, fm, fh, rimin, + & sfcflg, zfmin, sfcfrac, hpbl, zl(1:im,1), + & aphi16, aphi5, alpha_unstable, alpha_stable, + & tune_ocean_surface_layer, phim, phih, zol) ! do i=1,im if(pblflg(i)) then @@ -1564,3 +1552,61 @@ subroutine tridit(l,n,nt,cl,cm,cu,rt,au,at) c----------------------------------------------------------------------- return end + + subroutine get_similarity_parameters( + & im, rbsoil, fm, fh, rimin, + & sfcflg, zfmin, sfcfrac, hpbl, zl, + & aphi16, aphi5, alpha_unstable, alpha_stable, + & tune_ocean_surface_layer, phim, phih, zol + & ) + + use machine, only : kind_phys + + implicit none + ! input + integer, intent(in) :: im + real(kind=kind_phys), dimension(1:im), intent(in) :: + & rbsoil, fm, fh, hpbl, zl + real(kind=kind_phys), intent(in) :: rimin, zfmin, aphi16, aphi5, + & alpha_unstable, alpha_stable, sfcfrac + logical, dimension(1:im), intent(in) :: sfcflg + logical, intent(in) :: tune_ocean_surface_layer + + ! output + real(kind=kind_phys), dimension(1:im), intent(out) :: phim, phih, + & zol + + ! local + integer i + real :: zol1 + + do i=1,im + zol(i) = max(rbsoilv*fm(i)*fm(i)/fh(i),rimin) + if(sfcflg(i)) then + zol(i) = min(zol(i),-zfmin) + else + zol(i) = max(zol(i),zfmin) + endif + ! + zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i) + if(sfcflg(i)) then + tem = 1.0 / (1. - aphi16*zol1) + phih(i) = sqrt(tem) + if (tune_ocean_surface_layer) then + tem = 1.0 / (1. - alpha_unstable*zol1) + phim(i) = sqrt(sqrt(tem)) + else + phim(i) = sqrt(phih(i)) + end if + else + if (tune_ocean_surface_layer) then + phim(i) = 1. + alpha_stable*zol1 + else + phim(i) = 1. + aphi5*zol1 + end if + phih(i) = phim(i) + endif + enddo + return + end subroutine + diff --git a/gsmphys/satmedmfvdifq.f b/gsmphys/satmedmfvdifq.f index 4b233c04..4f99b6f2 100644 --- a/gsmphys/satmedmfvdifq.f +++ b/gsmphys/satmedmfvdifq.f @@ -56,7 +56,8 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke, & kinver,xkzm_mo,xkzm_ho,xkzm_ml,xkzm_hl,xkzm_mi,xkzm_hi, & xkzm_s,xkzinv,rlmx,zolcru,cs0, & do_dk_hb19,xkgdx,dspfac,bl_upfr,bl_dnfr,dkt_out, - & flux_up, flux_dn) + & flux_up, flux_dn, + & alpha_stable, alpha_unstable, tune_ocean_surface_layer) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -73,7 +74,8 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke, ! real(kind=kind_phys) delt, xkzm_mo, xkzm_ho, xkzm_s, dspfac, & bl_upfr, bl_dnfr, xkzm_ml, xkzm_hl, - & xkzm_mi, xkzm_hi + & xkzm_mi, xkzm_hi, alpha_stable, + & alpha_unstable real(kind=kind_phys) dv(im,km), du(im,km), & tdt(im,km), rtg(im,km,ntrac), & u1(ix,km), v1(ix,km), @@ -97,7 +99,7 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke, & rtg_in(im,km,ntrac) ! kgao note - q1 and rtg are local var now ! - logical dspheat, do_dk_hb19 + logical dspheat, do_dk_hb19, tune_ocean_surface_layer ! flag for tke dissipative heating real(kind=kind_phys)::dkt_out(im,km),flux_up(im,km),flux_dn(im,km) ! @@ -619,24 +621,10 @@ subroutine satmedmfvdifq(ix,im,km,ntrac,ntcw,ntiw,ntke, ! ! compute similarity parameters ! - do i=1,im - zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin) - if(sfcflg(i)) then - zol(i) = min(zol(i),-zfmin) - else - zol(i) = max(zol(i),zfmin) - endif -! - zol1 = zol(i)*sfcfrac*hpbl(i)/zl(i,1) - if(sfcflg(i)) then - tem = 1.0 / (1. - aphi16*zol1) - phih(i) = sqrt(tem) - phim(i) = sqrt(phih(i)) - else - phim(i) = 1. + aphi5*zol1 - phih(i) = phim(i) - endif - enddo + call get_similarity_parameters(im, rbsoil, fm, fh, rimin, + & sfcflg, zfmin, sfcfrac, hpbl, zl(1:im,1), + & aphi16, aphi5, alpha_unstable, alpha_stable, + & tune_ocean_surface_layer, phim, phih, zol) ! do i=1,im if(pblflg(i)) then From 87358fab0480e40731eb4055b9f13e2d5de73603 Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Thu, 15 May 2025 07:26:03 -0700 Subject: [PATCH 7/9] Fixed the atan functions from python lingo to fortran --- gsmphys/satmedmfvdiff.f | 2 +- gsmphys/sfc_diff_gfdl.f | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/gsmphys/satmedmfvdiff.f b/gsmphys/satmedmfvdiff.f index 259d0f43..e2452de2 100644 --- a/gsmphys/satmedmfvdiff.f +++ b/gsmphys/satmedmfvdiff.f @@ -1578,7 +1578,7 @@ subroutine get_similarity_parameters( ! local integer i - real :: zol1 + real(kind=kind_phys) :: zol1 do i=1,im zol(i) = max(rbsoilv*fm(i)*fm(i)/fh(i),rimin) diff --git a/gsmphys/sfc_diff_gfdl.f b/gsmphys/sfc_diff_gfdl.f index 305668d9..df9a5217 100644 --- a/gsmphys/sfc_diff_gfdl.f +++ b/gsmphys/sfc_diff_gfdl.f @@ -741,6 +741,6 @@ subroutine get_psim_unstable(zeta, z1i, z0, alpha, psi_m) Rz = (1. - alpha * zeta)**0.25 R0 = (1. - alpha * zeta0)**0.25 psi_m = log((1. + Rz)**2. * (1. + Rz**2.) / ((1. + R0)**2. - & *(1. + R0**2.))) + 2. * np.arctan(R0) - 2. * np.arctan(Rz) + & *(1. + R0**2.))) + 2. * atan(R0) - 2. * atan(Rz) return end subroutine From c78cc6d1f151c1a540e437b7f93466340650b9c6 Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Thu, 15 May 2025 20:00:16 -0700 Subject: [PATCH 8/9] Fixed some up --- gsmphys/satmedmfvdiff.f | 4 ++-- gsmphys/sfc_diff_gfdl.f | 3 +-- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/gsmphys/satmedmfvdiff.f b/gsmphys/satmedmfvdiff.f index e2452de2..dbb31c8d 100644 --- a/gsmphys/satmedmfvdiff.f +++ b/gsmphys/satmedmfvdiff.f @@ -1578,10 +1578,10 @@ subroutine get_similarity_parameters( ! local integer i - real(kind=kind_phys) :: zol1 + real(kind=kind_phys) :: zol1, tem do i=1,im - zol(i) = max(rbsoilv*fm(i)*fm(i)/fh(i),rimin) + zol(i) = max(rbsoil(i)*fm(i)*fm(i)/fh(i),rimin) if(sfcflg(i)) then zol(i) = min(zol(i),-zfmin) else diff --git a/gsmphys/sfc_diff_gfdl.f b/gsmphys/sfc_diff_gfdl.f index df9a5217..373d8f1f 100644 --- a/gsmphys/sfc_diff_gfdl.f +++ b/gsmphys/sfc_diff_gfdl.f @@ -566,7 +566,6 @@ subroutine monin_obukhov_similarity real(kind=kind_phys), parameter :: alpha=5., a0=-3.975 &, a1=12.32, alpha4_nonocean=4.0*alpha - &, alpha4_ocean = 4.*alpha_stable &, b1=-7.755, b2=6.041, alpha2=alpha+alpha, beta=1.0 &, a0p=-7.941, a1p=24.75, b1p=-8.705, b2p=7.899 &, ztmin1=-999.0, ca=.4 @@ -581,7 +580,7 @@ subroutine monin_obukhov_similarity z1i = 1.0 / z1 if (ilsimask == 0 .and. tune_ocean_surface_layer) then - alpha4 = alpha4_ocean + alpha4 = 4. * alpha_stable else alpha4 = alpha4_nonocean end if From 8c5be0edaed6c6248f0207234e2dbdf2caa0201a Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Fri, 16 May 2025 12:05:27 -0700 Subject: [PATCH 9/9] removed extra anpersand --- GFS_layer/GFS_typedefs.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GFS_layer/GFS_typedefs.F90 b/GFS_layer/GFS_typedefs.F90 index c2f53b74..561ab78c 100644 --- a/GFS_layer/GFS_typedefs.F90 +++ b/GFS_layer/GFS_typedefs.F90 @@ -2640,7 +2640,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & h2o_phys, pdfcld, shcnvcw, redrag, sfc_gfdl, z0s_max, & sfc_coupled, do_z0_moon, do_z0_hwrf15, do_z0_hwrf17, & do_z0_hwrf17_hwonly, wind_th_hwrf, alpha_stable, & - alpha_unstable, tune_ocean_surface_layer, & & + alpha_unstable, tune_ocean_surface_layer, & hybedmf, dspheat, lheatstrg, hour_canopy, afac_canopy, & cnvcld, no_pbl, xkzm_lim, xkzm_fac, xkgdx, & rlmn, rlmx, zolcru, cs0, &