From b5f8be6b4647586d4e22ea7ac218a5ad3f8cfad9 Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Mon, 14 Jul 2025 20:47:22 -0700 Subject: [PATCH 01/20] Fixed use of new parameters in gsmphys/satmedmfvdiff.f --- gsmphys/satmedmfvdiff.f | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/gsmphys/satmedmfvdiff.f b/gsmphys/satmedmfvdiff.f index 1b3c8698..be5b5f2a 100644 --- a/gsmphys/satmedmfvdiff.f +++ b/gsmphys/satmedmfvdiff.f @@ -797,18 +797,23 @@ subroutine satmedmfvdif(ix,im,km,ntrac,ntcw,ntiw,ntke, ! if (islimsk(i) == 0) then ckz(i,k) = ck1 + (ck0_o-ck1)*exp(ptem) + ckz(i,k) = min(ckz(i,k),ck0_o) + ckz(i,k) = max(ckz(i,k),ck1) else - ckz(i,k) = ck1 + (ck0_o-ck1)*exp(ptem) + ckz(i,k) = ck1 + (ck0-ck1)*exp(ptem) + ckz(i,k) = min(ckz(i,k),ck0) + ckz(i,k) = max(ckz(i,k),ck1) end if - ckz(i,k) = min(ckz(i,k),ck0) - ckz(i,k) = max(ckz(i,k),ck1) + if (islimsk(i) == 0) then chz(i,k) = ch1 + (ch0_o-ch1)*exp(ptem) + chz(i,k) = min(chz(i,k),ch0_o) + chz(i,k) = max(chz(i,k),ch1) else chz(i,k) = ch1 + (ch0-ch1)*exp(ptem) + chz(i,k) = min(chz(i,k),ch0) + chz(i,k) = max(chz(i,k),ch1) end if - chz(i,k) = min(chz(i,k),ch0) - chz(i,k) = max(chz(i,k),ch1) endif enddo enddo From b2435297c3470d95831b7aec2af1986def7a9c53 Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Tue, 15 Jul 2025 10:11:20 -0700 Subject: [PATCH 02/20] Fixed bug, which casued plume u flux, v flux, and eddy viscosity to be written out incocrrectly --- gsmphys/satmedmfvdiff.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gsmphys/satmedmfvdiff.f b/gsmphys/satmedmfvdiff.f index be5b5f2a..5773b823 100644 --- a/gsmphys/satmedmfvdiff.f +++ b/gsmphys/satmedmfvdiff.f @@ -36,7 +36,7 @@ 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,ck0_o,ce0_o, afrac_o,do_dk_hb19,xkzm_lim,xkgdx, - & rlmn, rlmx, cap_k0_land, dkt_out,flux_uup,flux_vup,dku_out, + & rlmn, rlmx, cap_k0_land, dkt_out, dku_out,flux_uup,flux_vup, & alpha_stable, alpha_unstable, tune_ocean_surface_layer) ! use machine , only : kind_phys From deb7add43277de82ebd8408dbaa15eda1b019c2f Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Wed, 6 Aug 2025 13:15:40 -0700 Subject: [PATCH 03/20] Update diagnost variable ZOL in code in case it changes in an iteration; add diagnostics for drag coefficient at height of first model level named cd; added fm10 from SL scheme as diagnostic output --- FV3GFS/FV3GFS_io.F90 | 35 +++++++++++++++++++++++++++----- GFS_layer/GFS_physics_driver.F90 | 9 ++++---- GFS_layer/GFS_typedefs.F90 | 4 ++++ gsmphys/sfc_diff_gfdl.f | 2 ++ 4 files changed, 41 insertions(+), 9 deletions(-) diff --git a/FV3GFS/FV3GFS_io.F90 b/FV3GFS/FV3GFS_io.F90 index aa178e60..dc186323 100644 --- a/FV3GFS/FV3GFS_io.F90 +++ b/FV3GFS/FV3GFS_io.F90 @@ -4792,6 +4792,31 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%zol(:) enddo + + idx = idx + 1 + Diag(idx)%axes = 2 + Diag(idx)%name = 'cd' + Diag(idx)%desc = 'Momentum drag coefficient' + Diag(idx)%unit = '1' + Diag(idx)%mod_name = 'gfs_phys' + Diag(idx)%intpl_method = 'bilinear' + allocate (Diag(idx)%data(nblks)) + do nb = 1,nblks + Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%cd(:) + enddo + + + idx = idx + 1 + Diag(idx)%axes = 2 + Diag(idx)%name = 'fm10' + Diag(idx)%desc = 'fm10 parameter from SL scheme at 10 meters' + Diag(idx)%unit = '1' + Diag(idx)%mod_name = 'gfs_sfc' + allocate (Diag(idx)%data(nblks)) + do nb = 1,nblks + Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%fm10(:) + enddo + idx = idx + 1 Diag(idx)%axes = 2 Diag(idx)%name = 'dpt2m' @@ -7000,7 +7025,7 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block Diag(idx)%axes = 2 Diag(idx)%name = 'f10m' Diag(idx)%desc = '10-meter wind speed divided by lowest model wind speed' - Diag(idx)%unit = 'N/A' + Diag(idx)%unit = '1' Diag(idx)%mod_name = 'gfs_sfc' allocate (Diag(idx)%data(nblks)) do nb = 1,nblks @@ -7033,7 +7058,7 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block Diag(idx)%axes = 2 Diag(idx)%name = 'ffhh' Diag(idx)%desc = 'fh parameter from PBL scheme' - Diag(idx)%unit = 'XXX' + Diag(idx)%unit = '1' Diag(idx)%mod_name = 'gfs_sfc' allocate (Diag(idx)%data(nblks)) do nb = 1,nblks @@ -7043,8 +7068,8 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block idx = idx + 1 Diag(idx)%axes = 2 Diag(idx)%name = 'ffmm' - Diag(idx)%desc = 'fm parameter from PBL scheme' - Diag(idx)%unit = 'XXX' + Diag(idx)%desc = 'fm parameter from SL scheme' + Diag(idx)%unit = '1' Diag(idx)%mod_name = 'gfs_sfc' allocate (Diag(idx)%data(nblks)) do nb = 1,nblks @@ -7055,7 +7080,7 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block Diag(idx)%axes = 2 Diag(idx)%name = 'uustar' Diag(idx)%desc = 'uustar surface frictional wind' - Diag(idx)%unit = 'XXX' + Diag(idx)%unit = 'm/s' Diag(idx)%mod_name = 'gfs_sfc' allocate (Diag(idx)%data(nblks)) do nb = 1,nblks diff --git a/GFS_layer/GFS_physics_driver.F90 b/GFS_layer/GFS_physics_driver.F90 index b5ec650d..4100203c 100644 --- a/GFS_layer/GFS_physics_driver.F90 +++ b/GFS_layer/GFS_physics_driver.F90 @@ -1202,7 +1202,7 @@ subroutine GFS_physics_driver & Sfcprop%charnock, & fm10_neutral, & Sfcprop%uustar, & - wind, Tbd%phy_f2d(1,Model%num_p2d), fm10, fh2, & + wind, Tbd%phy_f2d(1,Model%num_p2d), Diag%fm10, fh2, & sigmaf, vegtype, Sfcprop%shdmax, Model%ivegsrc, & tsurf, flag_iter, Model%redrag, Model%z0s_max, & Model%do_z0_moon, Model%do_z0_hwrf15, & @@ -1217,7 +1217,7 @@ subroutine GFS_physics_driver & Sfcprop%snowd, Sfcprop%tsfc, Sfcprop%zorl, cd, & cdq, rb, Statein%prsl(1,1), work3, islmsk, stress, & Sfcprop%ffmm, Sfcprop%ffhh, Sfcprop%uustar, & - wind, Tbd%phy_f2d(1,Model%num_p2d), fm10, fh2, & + wind, Tbd%phy_f2d(1,Model%num_p2d), Diag%fm10, fh2, & sigmaf, vegtype, Sfcprop%shdmax, Model%ivegsrc, & tsurf, flag_iter, Model%redrag, Model%czil_sfc, & Model%z0s_max, & @@ -1226,6 +1226,7 @@ subroutine GFS_physics_driver & Model%wind_th_hwrf, Diag%zol) endif Diag%rb = rb + Diag%cd = cd !endif !endif @@ -1484,7 +1485,7 @@ subroutine GFS_physics_driver & Statein%tgrs, Statein%qgrs, Sfcprop%tsfc, qss, & Sfcprop%f10m, Diag%u10m, Diag%v10m, & Sfcprop%t2m, Sfcprop%q2m, work3, evap, & - Sfcprop%ffmm, Sfcprop%ffhh, fm10, fh2, & + Sfcprop%ffmm, Sfcprop%ffhh, Diag%fm10, fh2, & Sfcprop%u10m, Sfcprop%v10m,& fm10_neutral, Diag%u10n, Diag%v10n, & Sfcprop%u10n, Sfcprop%v10n, Sfcprop%rhoa) @@ -3850,7 +3851,7 @@ subroutine GFS_physics_driver & Stateout%gt0, Stateout%gq0, Sfcprop%tsfc, qss, & Sfcprop%f10m, Diag%u10m, Diag%v10m, Sfcprop%t2m, & Sfcprop%q2m, work3, evap, Sfcprop%ffmm, & - Sfcprop%ffhh, fm10, fh2, & + Sfcprop%ffhh, Diag%fm10, fh2, & Sfcprop%u10m, Sfcprop%v10m,& fm10_neutral, Diag%u10n, Diag%v10n, & Sfcprop%u10n, Sfcprop%v10n, Sfcprop%rhoa) diff --git a/GFS_layer/GFS_typedefs.F90 b/GFS_layer/GFS_typedefs.F90 index 27451d63..4f74e990 100644 --- a/GFS_layer/GFS_typedefs.F90 +++ b/GFS_layer/GFS_typedefs.F90 @@ -4071,6 +4071,8 @@ subroutine diag_create (Diag, IM, Model) allocate (Diag%spfhmin (IM)) allocate (Diag%spfhmax (IM)) allocate (Diag%u10mmax (IM)) + allocate (Diag%fm10 (IM)) + allocate (Diag%cd (IM)) allocate (Diag%v10mmax (IM)) allocate (Diag%wind10mmax (IM)) allocate (Diag%rain (IM)) @@ -4376,6 +4378,8 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center) Diag%spfhmin = huge Diag%spfhmax = zero Diag%u10mmax = zero + Diag%fm10 = zero + Diag%cd = zero Diag%v10mmax = zero Diag%wind10mmax = zero Diag%rain = zero diff --git a/gsmphys/sfc_diff_gfdl.f b/gsmphys/sfc_diff_gfdl.f index a2c1732f..2d4fb009 100644 --- a/gsmphys/sfc_diff_gfdl.f +++ b/gsmphys/sfc_diff_gfdl.f @@ -617,6 +617,7 @@ subroutine monin_obukhov_similarity fhs = fh - ph hl1 = fms * fms * rb / fhs hl1 = min(max(hl1, ztmin1), ztmax1) + zol = hl1 endif ! ! second iteration @@ -648,6 +649,7 @@ subroutine monin_obukhov_similarity if(abs(olinf) <= tem1) then hlinf = -z1 / tem1 hlinf = min(max(hlinf,ztmin1),ztmax1) + zol = hlinf endif ! ! get pm and ph From d42b343fbcb7009c637c821bc536c2ea56e2c7dc Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Wed, 6 Aug 2025 14:04:46 -0700 Subject: [PATCH 04/20] Added pointers --- GFS_layer/GFS_typedefs.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/GFS_layer/GFS_typedefs.F90 b/GFS_layer/GFS_typedefs.F90 index 4f74e990..3b8ad332 100644 --- a/GFS_layer/GFS_typedefs.F90 +++ b/GFS_layer/GFS_typedefs.F90 @@ -1349,6 +1349,8 @@ module GFS_typedefs real (kind=kind_phys), pointer :: v1 (:) => null() !< layer 1 merdional wind (m/s) real (kind=kind_phys), pointer :: chh (:) => null() !< thermal exchange coefficient real (kind=kind_phys), pointer :: cmm (:) => null() !< momentum exchange coefficient + real (kind=kind_phys), pointer :: fm10 (:) => null() !< 10 meter non-dimensional wind shear + real (kind=kind_phys), pointer :: cd (:) => null() !< drag coefficient for momentum real (kind=kind_phys), pointer :: dlwsfci(:) => null() !< instantaneous sfc dnwd lw flux ( w/m**2 ) real (kind=kind_phys), pointer :: ulwsfci(:) => null() !< instantaneous sfc upwd lw flux ( w/m**2 ) real (kind=kind_phys), pointer :: dswsfci(:) => null() !< instantaneous sfc dnwd sw flux ( w/m**2 ) From d63ddc75a59243d4dc05ed69813552255f92142a Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Thu, 7 Aug 2025 16:49:07 -0700 Subject: [PATCH 05/20] changed the module for one diagnostic --- FV3GFS/FV3GFS_io.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/FV3GFS/FV3GFS_io.F90 b/FV3GFS/FV3GFS_io.F90 index dc186323..de240bb6 100644 --- a/FV3GFS/FV3GFS_io.F90 +++ b/FV3GFS/FV3GFS_io.F90 @@ -4811,7 +4811,7 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block Diag(idx)%name = 'fm10' Diag(idx)%desc = 'fm10 parameter from SL scheme at 10 meters' Diag(idx)%unit = '1' - Diag(idx)%mod_name = 'gfs_sfc' + Diag(idx)%mod_name = 'gfs_phys' allocate (Diag(idx)%data(nblks)) do nb = 1,nblks Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%fm10(:) From 27bbeac9fd5dc7bcde637f514fddc4c67aa93e42 Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Tue, 19 Aug 2025 13:51:51 -0700 Subject: [PATCH 06/20] Added a gust parameterization and diagnostic output --- FV3GFS/FV3GFS_io.F90 | 26 +++++++++++ GFS_layer/GFS_diagnostics.F90 | 19 ++++++++ GFS_layer/GFS_physics_driver.F90 | 18 +++++++- GFS_layer/GFS_typedefs.F90 | 21 +++++++-- gsmphys/gust_parameterization.f90 | 22 ++++++++++ gsmphys/sfc_diff_gfdl.f | 72 +++++++++++++++++++++++++++---- 6 files changed, 165 insertions(+), 13 deletions(-) create mode 100644 gsmphys/gust_parameterization.f90 diff --git a/FV3GFS/FV3GFS_io.F90 b/FV3GFS/FV3GFS_io.F90 index de240bb6..1e853822 100644 --- a/FV3GFS/FV3GFS_io.F90 +++ b/FV3GFS/FV3GFS_io.F90 @@ -4766,6 +4766,32 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block do nb = 1,nblks Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%v10n(:) enddo + + idx = idx + 1 + Diag(idx)%axes = 2 + Diag(idx)%name = '10mgust' + Diag(idx)%desc = '10-meter wind gust [m/s]' + Diag(idx)%unit = 'm/s' + Diag(idx)%mod_name = 'gfs_phys' + Diag(idx)%intpl_method = 'bilinear' + allocate (Diag(idx)%data(nblks)) + do nb = 1,nblks + Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%gust(:) + enddo + + + idx = idx + 1 + Diag(idx)%axes = 2 + Diag(idx)%name = '10mgustmax' + Diag(idx)%desc = 'maximum 10-meter wind gust [m/s]' + Diag(idx)%unit = 'm/s' + Diag(idx)%mod_name = 'gfs_phys' + Diag(idx)%intpl_method = 'bilinear' + allocate (Diag(idx)%data(nblks)) + do nb = 1,nblks + Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%gustmax(:) + enddo + !Sofar added: end idx = idx + 1 diff --git a/GFS_layer/GFS_diagnostics.F90 b/GFS_layer/GFS_diagnostics.F90 index 2ce2abcd..6362bd4b 100644 --- a/GFS_layer/GFS_diagnostics.F90 +++ b/GFS_layer/GFS_diagnostics.F90 @@ -1480,6 +1480,25 @@ subroutine diag_populate (IPD_Diag, Model, Statein, Stateout, Sfcprop, Coupling, IPD_Diag(idx)%data(nb)%var2p => Diag(nb)%wind10mmax enddo + + !---wind10mmax + idx = idx + 1 + IPD_Diag(idx)%name = 'gust10mmaxn' + IPD_Diag(idx)%output_name = 'gust10mmax' + IPD_Diag(idx)%mod_name = 'physics' + IPD_Diag(idx)%file_name = 'flx' + IPD_Diag(idx)%desc = 'maximum wind gust (m/s) at 10 m above ground' + IPD_Diag(idx)%unit = 'm/s' + IPD_Diag(idx)%type_stat_proc = 'max' + IPD_Diag(idx)%level_type = '10 m above grnd' + IPD_Diag(idx)%level = 1 + IPD_Diag(idx)%cnvfac = cn_one + IPD_Diag(idx)%zhour = Model%zhour + IPD_Diag(idx)%fcst_hour = Model%fhour + do nb = 1,nblks + IPD_Diag(idx)%data(nb)%var2p => Diag(nb)%gustmax + enddo + !---rain idx = idx + 1 IPD_Diag(idx)%name = 'rain' diff --git a/GFS_layer/GFS_physics_driver.F90 b/GFS_layer/GFS_physics_driver.F90 index 4100203c..7c386d00 100644 --- a/GFS_layer/GFS_physics_driver.F90 +++ b/GFS_layer/GFS_physics_driver.F90 @@ -484,7 +484,9 @@ subroutine GFS_physics_driver & dtsfc_cice, dqsfc_cice, dusfc_cice, dvsfc_cice, ulwsfc_cice, & tisfc_cice, tsea_cice, hice_cice, fice_cice, & !--- for CS-convection - wcbmax + wcbmax, & + !--- for gust + monin_obukhov_length logical, dimension(size(Grid%xlon,1)) :: & wet, dry, icy @@ -1209,7 +1211,8 @@ subroutine GFS_physics_driver & Model%do_z0_hwrf17, Model%do_z0_hwrf17_hwonly, & Model%wind_th_hwrf, & Model%alpha_stable, Model%alpha_unstable, & - Model%tune_ocean_surface_layer, Diag%zol) + Model%tune_ocean_surface_layer, Diag%zol, & + Model%add_w_freeconv) else ! GFS original sfc_diff modified by kgao call sfc_diff (im,Statein%pgr, Statein%ugrs, Statein%vgrs,& @@ -1490,6 +1493,17 @@ subroutine GFS_physics_driver & fm10_neutral, Diag%u10n, Diag%v10n, & Sfcprop%u10n, Sfcprop%v10n, Sfcprop%rhoa) !endif + monin_obukhov_length = Diag%zlvl / Diag%zol + call compute_gust(im, Diag%u10m, Diag%v10m, Sfcprop%uustar, & + monin_obukhov_length, Model%gust_parameter, & + Diag%gust) + do i=1, im + !find max wind gust + tem = Diag%gust(i) + if (tem > Diag%gustmax(i)) then + Diag%gustmax(i) = tem + endif + end do Tbd%phy_f2d(:,Model%num_p2d) = 0.0 diff --git a/GFS_layer/GFS_typedefs.F90 b/GFS_layer/GFS_typedefs.F90 index 3b8ad332..374f5749 100644 --- a/GFS_layer/GFS_typedefs.F90 +++ b/GFS_layer/GFS_typedefs.F90 @@ -220,6 +220,8 @@ module GFS_typedefs real (kind=kind_phys), pointer :: v10m (:) => null() !< neutral V10 wind in m/s ! Sofar added 11/17/23 real (kind=kind_phys), pointer :: u10n (:) => null() !< neutral U10 wind in m/s ! Sofar added 9/22/23 real (kind=kind_phys), pointer :: v10n (:) => null() !< neutral V10 wind in m/s ! Sofar added 9/22/23 + real (kind=kind_phys), pointer :: gust (:) => null() !< 10-m wind gust m/s + real (kind=kind_phys), pointer :: gustmax(:) => null() !< 10-m maximum wind gust m/s real (kind=kind_phys), pointer :: ztrl (:) => null() !< surface roughness for t and q in cm real (kind=kind_phys), pointer :: fice (:) => null() !< ice fraction over open water grid real (kind=kind_phys), pointer :: hprim (:) => null() !< topographic standard deviation in m ! @@ -680,6 +682,8 @@ module GFS_typedefs 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 :: add_w_freeconv !< add free convective velocity to wind speed in GFDL surface layer scheme + real(kind=kind_phys) :: gust_parameter !< constant in the parameterization of the 10-m wind gust 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) @@ -2406,7 +2410,9 @@ 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. !< 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 + real(kind=kind_phys) :: gust_parameter = 7.2 !< constant in the gust parameterization logical :: tune_ocean_surface_layer = .false. !< if true, use alphas above to modify the surface layer dimensionless gradients over oceans + logical :: add_w_freeconv = .false. !< add free convective velocity to wind speed in GFDL surface layer scheme 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 @@ -2640,9 +2646,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & ras, trans_trac, old_monin, cnvgwd, mstrat, moist_adj, & cscnv, cal_pre, do_aw, do_shoc, shocaftcnv, shoc_cld, & h2o_phys, pdfcld, shcnvcw, redrag, sfc_gfdl, z0s_max, & - do_z0_moon, do_z0_hwrf15, do_z0_hwrf17, & + 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, add_w_freeconv, & hybedmf, dspheat, lheatstrg, hour_canopy, afac_canopy, & cnvcld, no_pbl, xkzm_lim, xkzm_fac, xkgdx, & rlmn, rlmx, zolcru, cs0, & @@ -2656,6 +2662,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & cs_parm, flgmin, cgwf, ccwf, cdmbgwd, sup, ctei_rm, crtrh, & dlqf,rbcr,mix_precip,orogwd,myj_pbl,ysupbl,satmedmf, & cap_k0_land,do_dk_hb19,cloud_gfdl,gwd_p_crit, & + gust_parameter, & !--- Rayleigh friction prslrd0, ral_ts, & !--- mass flux deep convection @@ -2884,7 +2891,9 @@ 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%gust_parameter = gust_parameter Model%tune_ocean_surface_layer = tune_ocean_surface_layer + Model%add_w_freeconv = add_w_freeconv Model%hybedmf = hybedmf Model%myj_pbl = myj_pbl Model%ysupbl = ysupbl @@ -3595,7 +3604,9 @@ 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 *, ' gust_parameter : ', Model%gust_parameter + print *, ' tune_ocean_surface_layer : ', Model%tune_ocean_surface_layer + print *, ' add_w_freeconv : ', Model%add_w_freeconv print *, ' hybedmf : ', Model%hybedmf print *, ' myj_pbl : ', Model%myj_pbl print *, ' ysupbl : ', Model%ysupbl @@ -4077,6 +4088,8 @@ subroutine diag_create (Diag, IM, Model) allocate (Diag%cd (IM)) allocate (Diag%v10mmax (IM)) allocate (Diag%wind10mmax (IM)) + allocate (Diag%gust (IM)) + allocate (Diag%gustmax (IM)) allocate (Diag%rain (IM)) allocate (Diag%rainc (IM)) allocate (Diag%ice (IM)) @@ -4380,6 +4393,7 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center) Diag%spfhmin = huge Diag%spfhmax = zero Diag%u10mmax = zero + Diag%gustmax = zero Diag%fm10 = zero Diag%cd = zero Diag%v10mmax = zero @@ -4391,6 +4405,7 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center) Diag%graupel = zero !--- Out + Diag%gust = zero Diag%u10m = zero Diag%v10m = zero Diag%u10n = zero ! Sofar added: 10/19/23 diff --git a/gsmphys/gust_parameterization.f90 b/gsmphys/gust_parameterization.f90 new file mode 100644 index 00000000..a7c7c7fc --- /dev/null +++ b/gsmphys/gust_parameterization.f90 @@ -0,0 +1,22 @@ +subroutine compute_gust(im, u_10m, v_10m, ustar, monin_obukhov_length, gust_parameter, gust) + use machine, only : kind_phys + implicit none + + integer, intent(in) :: im + real(kind=kind_phys), intent(in) :: u_10m(im), v_10m(im), ustar(im), & + monin_obukhov_length(im), gust_parameter(im) + real(kind=kind_phys), intent(out) :: gust(im) + + real(kind=kind_phys), parameter :: PBL_HEIGHT = 1000.0_kind_phys ! in meters + + ! Follows the IFS gust implementation (see Eq. 3.109 in physics documentation for cycle Cy49r1) + + gust = sqrt(u_10m**2 + v_10m**2) + gust_parameter * ustar + + ! where Monin-Obukov lenght > 0, multiply by similarity function + where (monin_obukhov_length > 0.0_kind_phys) + gust = gust * & + (1.0_kind_phys - 1.0_kind_phys / 24.0_kind_phys * PBL_HEIGHT / monin_obukhov_length)**(1.0_kind_phys/3.0_kind_phys) + end where + +end subroutine compute_gust \ No newline at end of file diff --git a/gsmphys/sfc_diff_gfdl.f b/gsmphys/sfc_diff_gfdl.f index 2d4fb009..a1500739 100644 --- a/gsmphys/sfc_diff_gfdl.f +++ b/gsmphys/sfc_diff_gfdl.f @@ -21,7 +21,8 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, & do_z0_moon, do_z0_hwrf15, do_z0_hwrf17, & do_z0_hwrf17_hwonly, wind_th_hwrf, & alpha_stable, alpha_unstable, - & tune_ocean_surface_layer, zol) + & tune_ocean_surface_layer, zol, + & add_w_freeconv) ! oct 2019 - a clean and updated version by Kun Gao at GFDL (Kun.Gao@noaa.gov) @@ -42,15 +43,15 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, &, prsl1, prslki, stress &, fm, fh &, charnock ! Sofar added Spring 2023 - &, ustar, wind + &, ustar, wind, w_conv &, ddvel &, fm10, fh2, sigmaf, shdmax &, tsurf, snwdph - &, fm_neutral, fm10_neutral, zol! Sofar added Spring 2023 + &, fm_neutral, fm10_neutral, zol real(kind=kind_phys) :: ws1, ws10n, alpha_stable, alpha_unstable ! Sofar added 10/19/23 integer, dimension(im) ::vegtype, islimsk - logical flag_iter(im), tune_ocean_surface_layer + logical flag_iter(im), tune_ocean_surface_layer, add_w_freeconv logical redrag logical do_z0_moon, do_z0_hwrf15, do_z0_hwrf17 ! kgao &, do_z0_hwrf17_hwonly ! kgao @@ -87,15 +88,27 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, if(flag_iter(i)) then ! --- get variables at model lowest layer and surface (water/ice/land) - - wind(i) = max(sqrt(u1(i)*u1(i) + v1(i)*v1(i)) - & + max(0.0, min(ddvel(i), 30.0)), 1.0) + + wind(i) = sqrt(u1(i)*u1(i) + v1(i)*v1(i)) tem1 = 1.0 + rvrdm1 * max(q1(i),1.e-8) + th1 = t1(i) * prslki(i) thv1 = t1(i) * prslki(i) * tem1 tvs = 0.5 * (tsurf(i)+tskin(i)) * tem1 qs1 = fpvs(t1(i)) qs1 = max(1.0e-8, eps * qs1 / (prsl1(i) + epsm1 * qs1)) + ! update wind speed with free convective velocity + if (add_w_freeconv) then + dt1 = th1 - Sfcprop%tsfc(i) + dq1 = qs1 - Sfcprop%qsfc(i) + call compute_w_freeconv(th1, dt1, dq1, ustar(i), + & fh(i), w_conv(i)) + wind(i) = sqrt(wind(i) * wind(i) + w_conv(i) * w_conv(i)) + end if + + ! update wind speed with convective downdraft velocity + wind(i) = max(wind(i) + max(0.0, min(ddvel(i), 30.0)), 1.0) + !(sea/land/ice mask =0/1/2) if(islimsk(i) == 1 .or. islimsk(i) == 2) then ! over land or sea ice @@ -162,6 +175,7 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, ztmax = max(ztmax,1.0e-6) + ! --- call similarity call monin_obukhov_similarity @@ -169,7 +183,7 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, & ztmax, tvs, & 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) + & fm_neutral(i), fm10_neutral(i), & cm(i), ch(i), stress(i), ustar(i), zol(i)) elseif (islimsk(i) == 0) then ! over water @@ -731,3 +745,45 @@ subroutine get_psim_unstable(zeta, z1i, z0, alpha, psi_m) & *(1. + R0**2.))) + 2. * atan(R0) - 2. * atan(Rz) return end subroutine + + subroutine compute_w_freeconv(t1, dt1, dq1, ustar, fh, w_conv) + ! This routine computes the free convection velocity scale + ! following the equation + ! w_conv = (H g/T * Q)^(1/3) + ! see, e.g., Eq 3.20 in IFS physics doc cycle Cy49r1 + ! with Q the kinematic surface enthalpy flux which + ! is taken from the previous step + + ! inputs: + ! t1: pot. temperature at lowest model level (K) + ! dt1: pot. temperature difference between surface and lowest model level + ! dq1: spec. humidity difference between surface and lowest model level + ! ustar: friction velocity + ! fh: Similarity function for heat at lowest model level + + ! outputs: + ! w_conv: free convection velocity scale (m/s) + + + use machine , only : kind_phys + use physcons, grav => con_g, rvrdm1 => con_fvirt + implicit none + real(kind=kind_phys),intent(in) :: t1, dt1, dq1, ustar, fh + real(kind=kind_phys),intent(out) :: w_conv + real(kind=kind_phys), parameter :: kappa = 0.4, H = 1000. + real(kind_phys) :: tstar, qstar, flux + + ! compute surface flux using velocity, temperature, and humidity scales + ! temp and humidity scales can be expressed as + ! C_H / sqrt(C_M) dt1 and C_H / sqrt(C_M) dq1 which equals + ! kappa / fh * dt1 and kappa / fh * dq1, respectively, + ! with fh, surface T and surface qs from the previous time step. + + tstar = kappa / fh * dt1 + qstar = kappa / fh * dq1 + flux = ustar * tstar + rvrdm1 * t1 * ustar * qstar ! kinematic enthalpy flux (K m/s) + + if (flux < 0.) flux = 0. ! ensure w_conv is zero if not unstable + w_conv = ( H * (grav/t1) * flux )**(1./3.) + + end subroutine From bff205cccfa11a87e60b387614c92376a0f2d646 Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Tue, 19 Aug 2025 14:37:24 -0700 Subject: [PATCH 07/20] Fixed allocation issues --- GFS_layer/GFS_physics_driver.F90 | 2 +- gsmphys/sfc_diff_gfdl.f | 11 ++++++----- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/GFS_layer/GFS_physics_driver.F90 b/GFS_layer/GFS_physics_driver.F90 index 7c386d00..c559056f 100644 --- a/GFS_layer/GFS_physics_driver.F90 +++ b/GFS_layer/GFS_physics_driver.F90 @@ -1198,7 +1198,7 @@ subroutine GFS_physics_driver & ! a new and more flexible version of sfc_diff by kgao call sfc_diff_gfdl(im,Statein%pgr, Statein%ugrs, Statein%vgrs,& Statein%tgrs, Statein%qgrs, Diag%zlvl, Sfcprop%snowd, & - Sfcprop%tsfc, Sfcprop%zorl, Sfcprop%ztrl, cd, & + Sfcprop%tsfc, Sfcprop%qsfc, Sfcprop%zorl, Sfcprop%ztrl, cd, & cdq, rb, Statein%prsl(1,1), work3, islmsk, stress, & Sfcprop%ffmm, Sfcprop%ffhh, & Sfcprop%charnock, & diff --git a/gsmphys/sfc_diff_gfdl.f b/gsmphys/sfc_diff_gfdl.f index a1500739..464f8715 100644 --- a/gsmphys/sfc_diff_gfdl.f +++ b/gsmphys/sfc_diff_gfdl.f @@ -9,7 +9,7 @@ ! Christie Hegermiller, Sofar Ocean subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, - & snwdph,tskin,z0rl,ztrl,cm,ch,rb, + & snwdph,tskin,qss,z0rl,ztrl,cm,ch,rb, & prsl1,prslki,islimsk, & stress,fm,fh, & charnock, @@ -40,7 +40,7 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, real(kind=kind_phys), dimension(im)::ps, u1, v1, t1, q1, z1 &, tskin, z0rl, ztrl, cm, ch, rb - &, prsl1, prslki, stress + &, prsl1, prslki, stress, qss &, fm, fh &, charnock ! Sofar added Spring 2023 &, ustar, wind, w_conv @@ -67,7 +67,8 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, & hl110, hlt, hltinf, olinf, & restar, czilc, tem1, tem2, & u10m, v10m, ws10m, ws10m_moon, !kgao - & z0_1, zt_1, fm1, fh1, ustar_1, ztmax_1 !kgao + & z0_1, zt_1, fm1, fh1, ustar_1, ztmax_1,!kgao + & th1, dt1, dq1 ! real(kind=kind_phys),intent(in ) :: z0s_max, wind_th_hwrf ! kgao @@ -99,8 +100,8 @@ subroutine sfc_diff_gfdl(im,ps,u1,v1,t1,q1,z1, ! update wind speed with free convective velocity if (add_w_freeconv) then - dt1 = th1 - Sfcprop%tsfc(i) - dq1 = qs1 - Sfcprop%qsfc(i) + dt1 = th1 - tskin(i) + dq1 = qs1 - qss(i) call compute_w_freeconv(th1, dt1, dq1, ustar(i), & fh(i), w_conv(i)) wind(i) = sqrt(wind(i) * wind(i) + w_conv(i) * w_conv(i)) From 87f0077c6cc80f7f206a89f29a406c2d7c13e80d Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Tue, 19 Aug 2025 15:24:05 -0700 Subject: [PATCH 08/20] fixed bug --- GFS_layer/GFS_typedefs.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GFS_layer/GFS_typedefs.F90 b/GFS_layer/GFS_typedefs.F90 index 374f5749..dc4fd150 100644 --- a/GFS_layer/GFS_typedefs.F90 +++ b/GFS_layer/GFS_typedefs.F90 @@ -220,8 +220,6 @@ module GFS_typedefs real (kind=kind_phys), pointer :: v10m (:) => null() !< neutral V10 wind in m/s ! Sofar added 11/17/23 real (kind=kind_phys), pointer :: u10n (:) => null() !< neutral U10 wind in m/s ! Sofar added 9/22/23 real (kind=kind_phys), pointer :: v10n (:) => null() !< neutral V10 wind in m/s ! Sofar added 9/22/23 - real (kind=kind_phys), pointer :: gust (:) => null() !< 10-m wind gust m/s - real (kind=kind_phys), pointer :: gustmax(:) => null() !< 10-m maximum wind gust m/s real (kind=kind_phys), pointer :: ztrl (:) => null() !< surface roughness for t and q in cm real (kind=kind_phys), pointer :: fice (:) => null() !< ice fraction over open water grid real (kind=kind_phys), pointer :: hprim (:) => null() !< topographic standard deviation in m ! @@ -1308,6 +1306,8 @@ module GFS_typedefs real (kind=kind_phys), pointer :: runoff (:) => null() !< total water runoff real (kind=kind_phys), pointer :: ep (:) => null() !< potential evaporation real (kind=kind_phys), pointer :: cldwrk (:) => null() !< cloud workfunction (valid only with sas) + real (kind=kind_phys), pointer :: gust (:) => null() !< 10-m wind gust m/s + real (kind=kind_phys), pointer :: gustmax(:) => null() !< 10-m maximum wind gust m/s real (kind=kind_phys), pointer :: dugwd (:) => null() !< vertically integrated u change by OGWD real (kind=kind_phys), pointer :: dvgwd (:) => null() !< vertically integrated v change by OGWD real (kind=kind_phys), pointer :: psmean (:) => null() !< surface pressure (kPa) From eeec9287798a217fd15d95311dffdfe3960b1b79 Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Wed, 20 Aug 2025 11:05:18 -0700 Subject: [PATCH 09/20] compute H/L from zol and z without division; apply limiter to avoid factor going < 1 --- GFS_layer/GFS_physics_driver.F90 | 8 +++----- gsmphys/gust_parameterization.f90 | 17 +++++++++++------ 2 files changed, 14 insertions(+), 11 deletions(-) diff --git a/GFS_layer/GFS_physics_driver.F90 b/GFS_layer/GFS_physics_driver.F90 index c559056f..5a184c0e 100644 --- a/GFS_layer/GFS_physics_driver.F90 +++ b/GFS_layer/GFS_physics_driver.F90 @@ -484,9 +484,7 @@ subroutine GFS_physics_driver & dtsfc_cice, dqsfc_cice, dusfc_cice, dvsfc_cice, ulwsfc_cice, & tisfc_cice, tsea_cice, hice_cice, fice_cice, & !--- for CS-convection - wcbmax, & - !--- for gust - monin_obukhov_length + wcbmax logical, dimension(size(Grid%xlon,1)) :: & wet, dry, icy @@ -1493,9 +1491,9 @@ subroutine GFS_physics_driver & fm10_neutral, Diag%u10n, Diag%v10n, & Sfcprop%u10n, Sfcprop%v10n, Sfcprop%rhoa) !endif - monin_obukhov_length = Diag%zlvl / Diag%zol + call compute_gust(im, Diag%u10m, Diag%v10m, Sfcprop%uustar, & - monin_obukhov_length, Model%gust_parameter, & + Diag%zol, Diag%zlvl, Model%gust_parameter, & Diag%gust) do i=1, im !find max wind gust diff --git a/gsmphys/gust_parameterization.f90 b/gsmphys/gust_parameterization.f90 index a7c7c7fc..fe3f1113 100644 --- a/gsmphys/gust_parameterization.f90 +++ b/gsmphys/gust_parameterization.f90 @@ -1,22 +1,27 @@ -subroutine compute_gust(im, u_10m, v_10m, ustar, monin_obukhov_length, gust_parameter, gust) +subroutine compute_gust(im, u_10m, v_10m, ustar, zol, z1, gust_parameter, gust) use machine, only : kind_phys implicit none integer, intent(in) :: im real(kind=kind_phys), intent(in) :: u_10m(im), v_10m(im), ustar(im), & - monin_obukhov_length(im), gust_parameter(im) + zol(im), z1(im), gust_parameter(im) real(kind=kind_phys), intent(out) :: gust(im) real(kind=kind_phys), parameter :: PBL_HEIGHT = 1000.0_kind_phys ! in meters ! Follows the IFS gust implementation (see Eq. 3.109 in physics documentation for cycle Cy49r1) - gust = sqrt(u_10m**2 + v_10m**2) + gust_parameter * ustar + gust = gust_parameter * ustar - ! where Monin-Obukov lenght > 0, multiply by similarity function - where (monin_obukhov_length > 0.0_kind_phys) + H_o_z = PBL_HEIGHT / z1 + + ! where zol > 0, multiply by similarity function f(H/L) + where (zol > 0.0_kind_phys) gust = gust * & - (1.0_kind_phys - 1.0_kind_phys / 24.0_kind_phys * PBL_HEIGHT / monin_obukhov_length)**(1.0_kind_phys/3.0_kind_phys) + (max(1.0_kind_phys, 1.0_kind_phys - 1.0_kind_phys / 24.0_kind_phys * H_o_z * zol))**(1.0_kind_phys/3.0_kind_phys) end where + + gust = gust + sqrt(u_10m**2 + v_10m**2) + end subroutine compute_gust \ No newline at end of file From 4f176e515c2d26819ab9415557bb8d0715288e87 Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Wed, 20 Aug 2025 11:29:03 -0700 Subject: [PATCH 10/20] fixed something --- gsmphys/gust_parameterization.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/gsmphys/gust_parameterization.f90 b/gsmphys/gust_parameterization.f90 index fe3f1113..d4bddec0 100644 --- a/gsmphys/gust_parameterization.f90 +++ b/gsmphys/gust_parameterization.f90 @@ -8,6 +8,7 @@ subroutine compute_gust(im, u_10m, v_10m, ustar, zol, z1, gust_parameter, gust) real(kind=kind_phys), intent(out) :: gust(im) real(kind=kind_phys), parameter :: PBL_HEIGHT = 1000.0_kind_phys ! in meters + real(kind=kind_phys) :: H_o_z ! Follows the IFS gust implementation (see Eq. 3.109 in physics documentation for cycle Cy49r1) From f8b51d8b8accf05b0f75a359e58ac13f6eb423a7 Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Wed, 20 Aug 2025 11:44:23 -0700 Subject: [PATCH 11/20] fixed dimensions --- gsmphys/gust_parameterization.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gsmphys/gust_parameterization.f90 b/gsmphys/gust_parameterization.f90 index d4bddec0..809de186 100644 --- a/gsmphys/gust_parameterization.f90 +++ b/gsmphys/gust_parameterization.f90 @@ -8,7 +8,7 @@ subroutine compute_gust(im, u_10m, v_10m, ustar, zol, z1, gust_parameter, gust) real(kind=kind_phys), intent(out) :: gust(im) real(kind=kind_phys), parameter :: PBL_HEIGHT = 1000.0_kind_phys ! in meters - real(kind=kind_phys) :: H_o_z + real(kind=kind_phys) :: H_o_z(im) ! Follows the IFS gust implementation (see Eq. 3.109 in physics documentation for cycle Cy49r1) From 817df1caf9763a63a2d4237864dd06bd846660b9 Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Wed, 20 Aug 2025 15:10:25 -0700 Subject: [PATCH 12/20] fixed wrong limiter --- gsmphys/gust_parameterization.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gsmphys/gust_parameterization.f90 b/gsmphys/gust_parameterization.f90 index 809de186..d77ac73b 100644 --- a/gsmphys/gust_parameterization.f90 +++ b/gsmphys/gust_parameterization.f90 @@ -19,7 +19,7 @@ subroutine compute_gust(im, u_10m, v_10m, ustar, zol, z1, gust_parameter, gust) ! where zol > 0, multiply by similarity function f(H/L) where (zol > 0.0_kind_phys) gust = gust * & - (max(1.0_kind_phys, 1.0_kind_phys - 1.0_kind_phys / 24.0_kind_phys * H_o_z * zol))**(1.0_kind_phys/3.0_kind_phys) + (max(0.0_kind_phys, 1.0_kind_phys - 1.0_kind_phys / 24.0_kind_phys * H_o_z * zol))**(1.0_kind_phys/3.0_kind_phys) end where From 5d539c02c639aec1267e6b37ac0a6dcef6438782 Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Wed, 20 Aug 2025 15:34:51 -0700 Subject: [PATCH 13/20] Turned gust computation into loop --- gsmphys/gust_parameterization.f90 | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/gsmphys/gust_parameterization.f90 b/gsmphys/gust_parameterization.f90 index d77ac73b..99bb6dd8 100644 --- a/gsmphys/gust_parameterization.f90 +++ b/gsmphys/gust_parameterization.f90 @@ -9,20 +9,22 @@ subroutine compute_gust(im, u_10m, v_10m, ustar, zol, z1, gust_parameter, gust) real(kind=kind_phys), parameter :: PBL_HEIGHT = 1000.0_kind_phys ! in meters real(kind=kind_phys) :: H_o_z(im) + integer :: i ! Follows the IFS gust implementation (see Eq. 3.109 in physics documentation for cycle Cy49r1) - gust = gust_parameter * ustar + do i = 1,im + gust(i) = gust_parameter * ustar(i) - H_o_z = PBL_HEIGHT / z1 + H_o_z(i) = PBL_HEIGHT / z1(i) - ! where zol > 0, multiply by similarity function f(H/L) - where (zol > 0.0_kind_phys) - gust = gust * & - (max(0.0_kind_phys, 1.0_kind_phys - 1.0_kind_phys / 24.0_kind_phys * H_o_z * zol))**(1.0_kind_phys/3.0_kind_phys) - end where + ! where zol > 0, multiply by similarity function f(H/L) + if (zol(i) > 0.0_kind_phys) then + gust(i) = gust(i) * (max(0.0_kind_phys, 1.0_kind_phys - 1.0_kind_phys / 24.0_kind_phys * H_o_z(i) * zol(i)))**(1.0_kind_phys/3.0_kind_phys) + end if + gust(i) = gust(i) + sqrt(u_10m(i)**2 + v_10m(i)**2) + end do - - gust = gust + sqrt(u_10m**2 + v_10m**2) + end subroutine compute_gust \ No newline at end of file From 3f20bffed1ccbd96da60ce3fd2e176a00538202c Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Wed, 20 Aug 2025 15:58:07 -0700 Subject: [PATCH 14/20] Changed allocation of dimension --- gsmphys/gust_parameterization.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/gsmphys/gust_parameterization.f90 b/gsmphys/gust_parameterization.f90 index 99bb6dd8..ea5912a0 100644 --- a/gsmphys/gust_parameterization.f90 +++ b/gsmphys/gust_parameterization.f90 @@ -3,17 +3,17 @@ subroutine compute_gust(im, u_10m, v_10m, ustar, zol, z1, gust_parameter, gust) implicit none integer, intent(in) :: im - real(kind=kind_phys), intent(in) :: u_10m(im), v_10m(im), ustar(im), & - zol(im), z1(im), gust_parameter(im) - real(kind=kind_phys), intent(out) :: gust(im) + real(kind=kind_phys), dimension(im), intent(in) :: u_10m, v_10m, ustar, & + zol, z1, gust_parameter + real(kind=kind_phys), dimension(im), intent(out) :: gust real(kind=kind_phys), parameter :: PBL_HEIGHT = 1000.0_kind_phys ! in meters - real(kind=kind_phys) :: H_o_z(im) + real(kind=kind_phys), dimension(im) :: H_o_z integer :: i ! Follows the IFS gust implementation (see Eq. 3.109 in physics documentation for cycle Cy49r1) - do i = 1,im + do i=1,im gust(i) = gust_parameter * ustar(i) H_o_z(i) = PBL_HEIGHT / z1(i) From d56fa0b4aea3e3680550616b84b6338cfeb20c7e Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Wed, 20 Aug 2025 17:14:08 -0700 Subject: [PATCH 15/20] fixed bug --- gsmphys/gust_parameterization.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/gsmphys/gust_parameterization.f90 b/gsmphys/gust_parameterization.f90 index ea5912a0..1d34c07b 100644 --- a/gsmphys/gust_parameterization.f90 +++ b/gsmphys/gust_parameterization.f90 @@ -4,7 +4,8 @@ subroutine compute_gust(im, u_10m, v_10m, ustar, zol, z1, gust_parameter, gust) integer, intent(in) :: im real(kind=kind_phys), dimension(im), intent(in) :: u_10m, v_10m, ustar, & - zol, z1, gust_parameter + zol, z1 + real(kind=kind_phys), intent(in) :: gust_parameter real(kind=kind_phys), dimension(im), intent(out) :: gust real(kind=kind_phys), parameter :: PBL_HEIGHT = 1000.0_kind_phys ! in meters From 14dd7899e397839cd54eea97f871ee2943dbf2e5 Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Fri, 5 Sep 2025 10:53:31 -0700 Subject: [PATCH 16/20] First stab at implementing the 100-m wind components to SHiELD physics, compute 100-m inst and max gust --- FV3GFS/FV3GFS_io.F90 | 53 +++++++++++++++++++++++++++++-- GFS_layer/GFS_diagnostics.F90 | 20 +++++++++++- GFS_layer/GFS_physics_driver.F90 | 51 +++++++++++++++++++++++++---- GFS_layer/GFS_typedefs.F90 | 34 ++++++++++++++------ gsmphys/gust_parameterization.f90 | 6 ++-- 5 files changed, 142 insertions(+), 22 deletions(-) diff --git a/FV3GFS/FV3GFS_io.F90 b/FV3GFS/FV3GFS_io.F90 index 1e853822..36e36bda 100644 --- a/FV3GFS/FV3GFS_io.F90 +++ b/FV3GFS/FV3GFS_io.F90 @@ -4742,6 +4742,30 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%v10m(:) enddo + idx = idx + 1 + Diag(idx)%axes = 2 + Diag(idx)%name = 'u100m' + Diag(idx)%desc = '100 meter u wind [m/s]' + Diag(idx)%unit = 'm/s' + Diag(idx)%mod_name = 'gfs_phys' + Diag(idx)%intpl_method = 'vector_bilinear' + allocate (Diag(idx)%data(nblks)) + do nb = 1,nblks + Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%u100m(:) + enddo + + idx = idx + 1 + Diag(idx)%axes = 2 + Diag(idx)%name = 'v100m' + Diag(idx)%desc = '100 meter v wind [m/s]' + Diag(idx)%unit = 'm/s' + Diag(idx)%mod_name = 'gfs_phys' + Diag(idx)%intpl_method = 'vector_bilinear' + allocate (Diag(idx)%data(nblks)) + do nb = 1,nblks + Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%v100m(:) + enddo + !Sofar added: start: (12/20/23) idx = idx + 1 Diag(idx)%axes = 2 @@ -4776,7 +4800,7 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block Diag(idx)%intpl_method = 'bilinear' allocate (Diag(idx)%data(nblks)) do nb = 1,nblks - Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%gust(:) + Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%gust10m(:) enddo @@ -4789,7 +4813,32 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block Diag(idx)%intpl_method = 'bilinear' allocate (Diag(idx)%data(nblks)) do nb = 1,nblks - Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%gustmax(:) + Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%gustmax10m(:) + enddo + + idx = idx + 1 + Diag(idx)%axes = 2 + Diag(idx)%name = '100mgust' + Diag(idx)%desc = '100-meter wind gust [m/s]' + Diag(idx)%unit = 'm/s' + Diag(idx)%mod_name = 'gfs_phys' + Diag(idx)%intpl_method = 'bilinear' + allocate (Diag(idx)%data(nblks)) + do nb = 1,nblks + Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%gust100m(:) + enddo + + + idx = idx + 1 + Diag(idx)%axes = 2 + Diag(idx)%name = '100mgustmax' + Diag(idx)%desc = 'maximum 100-meter wind gust [m/s]' + Diag(idx)%unit = 'm/s' + Diag(idx)%mod_name = 'gfs_phys' + Diag(idx)%intpl_method = 'bilinear' + allocate (Diag(idx)%data(nblks)) + do nb = 1,nblks + Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%gustmax100m(:) enddo !Sofar added: end diff --git a/GFS_layer/GFS_diagnostics.F90 b/GFS_layer/GFS_diagnostics.F90 index 6362bd4b..79302740 100644 --- a/GFS_layer/GFS_diagnostics.F90 +++ b/GFS_layer/GFS_diagnostics.F90 @@ -1496,7 +1496,25 @@ subroutine diag_populate (IPD_Diag, Model, Statein, Stateout, Sfcprop, Coupling, IPD_Diag(idx)%zhour = Model%zhour IPD_Diag(idx)%fcst_hour = Model%fhour do nb = 1,nblks - IPD_Diag(idx)%data(nb)%var2p => Diag(nb)%gustmax + IPD_Diag(idx)%data(nb)%var2p => Diag(nb)%gustmax10m + enddo + + !---wind100mmax + idx = idx + 1 + IPD_Diag(idx)%name = 'gust100mmaxn' + IPD_Diag(idx)%output_name = 'gust100mmax' + IPD_Diag(idx)%mod_name = 'physics' + IPD_Diag(idx)%file_name = 'flx' + IPD_Diag(idx)%desc = 'maximum wind gust (m/s) at 100 m above ground' + IPD_Diag(idx)%unit = 'm/s' + IPD_Diag(idx)%type_stat_proc = 'max' + IPD_Diag(idx)%level_type = '100 m above grnd' + IPD_Diag(idx)%level = 1 + IPD_Diag(idx)%cnvfac = cn_one + IPD_Diag(idx)%zhour = Model%zhour + IPD_Diag(idx)%fcst_hour = Model%fhour + do nb = 1,nblks + IPD_Diag(idx)%data(nb)%var2p => Diag(nb)%gustmax100m enddo !---rain diff --git a/GFS_layer/GFS_physics_driver.F90 b/GFS_layer/GFS_physics_driver.F90 index 5a184c0e..1779bd32 100644 --- a/GFS_layer/GFS_physics_driver.F90 +++ b/GFS_layer/GFS_physics_driver.F90 @@ -1491,15 +1491,24 @@ subroutine GFS_physics_driver & fm10_neutral, Diag%u10n, Diag%v10n, & Sfcprop%u10n, Sfcprop%v10n, Sfcprop%rhoa) !endif - + ! get 100-m wind components using linear interpolation same as in FV3, then get wind gust at 10 and 100 m + call interpolate_z(im, npz, 100., wz, Statein%ugrs, Diag%u100m) + call interpolate_z(im, npz, 100., wz, Statein%vgrs, Diag%v100m) call compute_gust(im, Diag%u10m, Diag%v10m, Sfcprop%uustar, & Diag%zol, Diag%zlvl, Model%gust_parameter, & - Diag%gust) + Diag%gust10m) + call compute_gust(im, Diag%u100m, Diag%v100m, Sfcprop%uustar, & + Diag%zol, Diag%zlvl, Model%gust_parameter, & + Diag%gust100m) do i=1, im - !find max wind gust - tem = Diag%gust(i) - if (tem > Diag%gustmax(i)) then - Diag%gustmax(i) = tem + !find max wind gusts + tem = Diag%gust10m(i) + if (tem > Diag%gustmax10m(i)) then + Diag%gustmax10m(i) = tem + endif + tem = Diag%gust100m(i) + if (tem > Diag%gustmax100m(i)) then + Diag%gustmax100m(i) = tem endif end do @@ -4335,6 +4344,36 @@ subroutine compute_diagnostics_with_scaled_co2(Model, Statein, Sfcprop, Coupling enddo end subroutine compute_diagnostics_with_scaled_co2 + + subroutine interpolate_z(im, km, zl, hght, a3, a2) + implicit none + integer, intent(in):: im, km + real, intent(in):: hght(im,km+1) ! hght(k) > hght(k+1) + real, intent(in):: a3(im,km) + real, intent(in):: zl + real, intent(out):: a2(im) + ! local: + real zm(km) + integer i,k + + do i=1,im + do k=1,km + zm(k) = 0.5*(hght(i,k)+hght(i,k+1)) + enddo + if( zl >= zm(1) ) then + a2(i,j) = a3(i,1) + elseif ( zl <= zm(km) ) then + a2(i,j) = a3(i,km) + else + do k=1,km-1 + if( zl <= zm(k) .and. zl >= zm(k+1) ) then + a2(i) = a3(i,k) + (a3(i,k+1)-a3(i,k))*(zm(k)-zl)/(zm(k)-zm(k+1)) + exit + endif + enddo + endif + end do + end subroutine interpolate_z !> @} end module module_physics_driver diff --git a/GFS_layer/GFS_typedefs.F90 b/GFS_layer/GFS_typedefs.F90 index dc4fd150..cfc92108 100644 --- a/GFS_layer/GFS_typedefs.F90 +++ b/GFS_layer/GFS_typedefs.F90 @@ -216,8 +216,10 @@ module GFS_typedefs real (kind=kind_phys), pointer :: zorll (:) => null() !< land surface roughness in cm real (kind=kind_phys), pointer :: charnock (:) => null() !< Charnock parameter ! Sofar added Spring 2023 real (kind=kind_phys), pointer :: rhoa (:) => null() !< surface air density in kg/m3 ! Sofar added 9/22/23 - real (kind=kind_phys), pointer :: u10m (:) => null() !< neutral U10 wind in m/s ! Sofar added 11/17/23 - real (kind=kind_phys), pointer :: v10m (:) => null() !< neutral V10 wind in m/s ! Sofar added 11/17/23 + real (kind=kind_phys), pointer :: u10m (:) => null() !< U10 wind in m/s ! Sofar added 11/17/23 + real (kind=kind_phys), pointer :: v10m (:) => null() !< V10 wind in m/s ! Sofar added 11/17/23 + real (kind=kind_phys), pointer :: u100m (:) => null() !< U100 wind in m/s ! Sofar added 11/17/23 + real (kind=kind_phys), pointer :: v100m (:) => null() !< V100 wind in m/s ! Sofar added 11/17/23 real (kind=kind_phys), pointer :: u10n (:) => null() !< neutral U10 wind in m/s ! Sofar added 9/22/23 real (kind=kind_phys), pointer :: v10n (:) => null() !< neutral V10 wind in m/s ! Sofar added 9/22/23 real (kind=kind_phys), pointer :: ztrl (:) => null() !< surface roughness for t and q in cm @@ -681,7 +683,7 @@ module GFS_typedefs 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 :: add_w_freeconv !< add free convective velocity to wind speed in GFDL surface layer scheme - real(kind=kind_phys) :: gust_parameter !< constant in the parameterization of the 10-m wind gust + real(kind=kind_phys) :: gust_parameter !< constant in the parameterization of the gust 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) @@ -1306,8 +1308,10 @@ module GFS_typedefs real (kind=kind_phys), pointer :: runoff (:) => null() !< total water runoff real (kind=kind_phys), pointer :: ep (:) => null() !< potential evaporation real (kind=kind_phys), pointer :: cldwrk (:) => null() !< cloud workfunction (valid only with sas) - real (kind=kind_phys), pointer :: gust (:) => null() !< 10-m wind gust m/s - real (kind=kind_phys), pointer :: gustmax(:) => null() !< 10-m maximum wind gust m/s + real (kind=kind_phys), pointer :: gust10m (:) => null() !< 10-m wind gust m/s + real (kind=kind_phys), pointer :: gustmax10m(:) => null() !< 10-m maximum wind gust m/s + real (kind=kind_phys), pointer :: gust100m (:) => null() !< 100-m wind gust m/s + real (kind=kind_phys), pointer :: gustmax100m(:)=> null() !< 100-m maximum wind gust m/s real (kind=kind_phys), pointer :: dugwd (:) => null() !< vertically integrated u change by OGWD real (kind=kind_phys), pointer :: dvgwd (:) => null() !< vertically integrated v change by OGWD real (kind=kind_phys), pointer :: psmean (:) => null() !< surface pressure (kPa) @@ -1333,6 +1337,8 @@ module GFS_typedefs ! Output - only in physics real (kind=kind_phys), pointer :: u10m (:) => null() !< 10 meter u/v wind speed real (kind=kind_phys), pointer :: v10m (:) => null() !< 10 meter u/v wind speed + real (kind=kind_phys), pointer :: u100m (:) => null() !< 100 meter u/v wind speed + real (kind=kind_phys), pointer :: v100m (:) => null() !< 100 meter u/v wind speed real (kind=kind_phys), pointer :: u10n (:) => null() !< 10 meter u/v neutral wind speed !Sofar added 10/19/23 real (kind=kind_phys), pointer :: v10n (:) => null() !< 10 meter u/v neutral wind speed !Sofar added 10/19/23 real (kind=kind_phys), pointer :: zol (:) => null() !< Dimensionless surface layer stability @@ -4087,9 +4093,11 @@ subroutine diag_create (Diag, IM, Model) allocate (Diag%fm10 (IM)) allocate (Diag%cd (IM)) allocate (Diag%v10mmax (IM)) - allocate (Diag%wind10mmax (IM)) - allocate (Diag%gust (IM)) - allocate (Diag%gustmax (IM)) + allocate (Diag%wind10mmax (IM)) + allocate (Diag%gust10m (IM)) + allocate (Diag%gustmax10m (IM)) + allocate (Diag%gust100m (IM)) + allocate (Diag%gustmax100m (IM)) allocate (Diag%rain (IM)) allocate (Diag%rainc (IM)) allocate (Diag%ice (IM)) @@ -4103,6 +4111,8 @@ subroutine diag_create (Diag, IM, Model) allocate (Diag%totgrpb (IM)) allocate (Diag%u10m (IM)) allocate (Diag%v10m (IM)) + allocate (Diag%u100m (IM)) + allocate (Diag%v100m (IM)) allocate (Diag%u10n (IM)) ! Sofar added: 10/19/13 allocate (Diag%v10n (IM)) ! Sofar added: 10/19/13 allocate (Diag%zol (IM)) @@ -4393,7 +4403,8 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center) Diag%spfhmin = huge Diag%spfhmax = zero Diag%u10mmax = zero - Diag%gustmax = zero + Diag%gustmax10m = zero + Diag%gustmax100m = zero Diag%fm10 = zero Diag%cd = zero Diag%v10mmax = zero @@ -4405,9 +4416,12 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center) Diag%graupel = zero !--- Out - Diag%gust = zero + Diag%gust10m = zero + Diag%gust100m = zero Diag%u10m = zero Diag%v10m = zero + Diag%u100m = zero + Diag%v100m = zero Diag%u10n = zero ! Sofar added: 10/19/23 Diag%v10n = zero ! Sofar added: 10/19/23 Diag%zol = zero diff --git a/gsmphys/gust_parameterization.f90 b/gsmphys/gust_parameterization.f90 index 1d34c07b..9fccbf67 100644 --- a/gsmphys/gust_parameterization.f90 +++ b/gsmphys/gust_parameterization.f90 @@ -1,9 +1,9 @@ -subroutine compute_gust(im, u_10m, v_10m, ustar, zol, z1, gust_parameter, gust) +subroutine compute_gust(im, u, v, ustar, zol, z1, gust_parameter, gust) use machine, only : kind_phys implicit none integer, intent(in) :: im - real(kind=kind_phys), dimension(im), intent(in) :: u_10m, v_10m, ustar, & + real(kind=kind_phys), dimension(im), intent(in) :: u, v, ustar, & zol, z1 real(kind=kind_phys), intent(in) :: gust_parameter real(kind=kind_phys), dimension(im), intent(out) :: gust @@ -23,7 +23,7 @@ subroutine compute_gust(im, u_10m, v_10m, ustar, zol, z1, gust_parameter, gust) if (zol(i) > 0.0_kind_phys) then gust(i) = gust(i) * (max(0.0_kind_phys, 1.0_kind_phys - 1.0_kind_phys / 24.0_kind_phys * H_o_z(i) * zol(i)))**(1.0_kind_phys/3.0_kind_phys) end if - gust(i) = gust(i) + sqrt(u_10m(i)**2 + v_10m(i)**2) + gust(i) = gust(i) + sqrt(u(i)**2 + v(i)**2) end do From fc00b48d679198c6209d849ad534e54d63179fd4 Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Fri, 5 Sep 2025 12:13:42 -0700 Subject: [PATCH 17/20] Implemented the interpolation in height --- GFS_layer/GFS_physics_driver.F90 | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/GFS_layer/GFS_physics_driver.F90 b/GFS_layer/GFS_physics_driver.F90 index 1779bd32..dc75d31a 100644 --- a/GFS_layer/GFS_physics_driver.F90 +++ b/GFS_layer/GFS_physics_driver.F90 @@ -1492,8 +1492,8 @@ subroutine GFS_physics_driver & Sfcprop%u10n, Sfcprop%v10n, Sfcprop%rhoa) !endif ! get 100-m wind components using linear interpolation same as in FV3, then get wind gust at 10 and 100 m - call interpolate_z(im, npz, 100., wz, Statein%ugrs, Diag%u100m) - call interpolate_z(im, npz, 100., wz, Statein%vgrs, Diag%v100m) + call interpolate_z(im, npz, 100., Statein%phii, Statein%ugrs, Diag%u100m) + call interpolate_z(im, npz, 100., Statein%phii, Statein%vgrs, Diag%v100m) call compute_gust(im, Diag%u10m, Diag%v10m, Sfcprop%uustar, & Diag%zol, Diag%zlvl, Model%gust_parameter, & Diag%gust10m) @@ -4345,17 +4345,27 @@ subroutine compute_diagnostics_with_scaled_co2(Model, Statein, Sfcprop, Coupling end subroutine compute_diagnostics_with_scaled_co2 - subroutine interpolate_z(im, km, zl, hght, a3, a2) + subroutine interpolate_z(im, km, zl, phi, a3_in, a2) + use physcons, only: con_g implicit none integer, intent(in):: im, km - real, intent(in):: hght(im,km+1) ! hght(k) > hght(k+1) - real, intent(in):: a3(im,km) + real, intent(in):: phi(im,km+1) ! phi(k+1) > phi(k) + real, intent(in):: a3_in(im,km) real, intent(in):: zl real, intent(out):: a2(im) ! local: + real hght(im,km+1), a3(im, km) ! hght(k) > hght(k+1) real zm(km) integer i,k + ! flip z and a3 + do k=1,km+1 + do i=1,im + hght(i,k) = phi(i,km+2-k) / con_g + if (k <= km) a3(i,k) = a3_in(i,km+1-k) + enddo + enddo + do i=1,im do k=1,km zm(k) = 0.5*(hght(i,k)+hght(i,k+1)) From c8335b105121585d751cab51d8852696e0156114 Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Fri, 5 Sep 2025 12:43:17 -0700 Subject: [PATCH 18/20] fixed some issues --- GFS_layer/GFS_physics_driver.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/GFS_layer/GFS_physics_driver.F90 b/GFS_layer/GFS_physics_driver.F90 index dc75d31a..31de03f8 100644 --- a/GFS_layer/GFS_physics_driver.F90 +++ b/GFS_layer/GFS_physics_driver.F90 @@ -1492,8 +1492,8 @@ subroutine GFS_physics_driver & Sfcprop%u10n, Sfcprop%v10n, Sfcprop%rhoa) !endif ! get 100-m wind components using linear interpolation same as in FV3, then get wind gust at 10 and 100 m - call interpolate_z(im, npz, 100., Statein%phii, Statein%ugrs, Diag%u100m) - call interpolate_z(im, npz, 100., Statein%phii, Statein%vgrs, Diag%v100m) + call interpolate_z(im, levs, 100., Statein%phii, Statein%ugrs, Diag%u100m) + call interpolate_z(im, levs, 100., Statein%phii, Statein%vgrs, Diag%v100m) call compute_gust(im, Diag%u10m, Diag%v10m, Sfcprop%uustar, & Diag%zol, Diag%zlvl, Model%gust_parameter, & Diag%gust10m) @@ -4371,9 +4371,9 @@ subroutine interpolate_z(im, km, zl, phi, a3_in, a2) zm(k) = 0.5*(hght(i,k)+hght(i,k+1)) enddo if( zl >= zm(1) ) then - a2(i,j) = a3(i,1) + a2(i) = a3(i,1) elseif ( zl <= zm(km) ) then - a2(i,j) = a3(i,km) + a2(i) = a3(i,km) else do k=1,km-1 if( zl <= zm(k) .and. zl >= zm(k+1) ) then From 7afadefc6fca4a27b904442488c59b0999b13fae Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Mon, 8 Sep 2025 10:18:55 -0700 Subject: [PATCH 19/20] added instantaeous precip rates --- FV3GFS/FV3GFS_io.F90 | 27 +++++++++++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) diff --git a/FV3GFS/FV3GFS_io.F90 b/FV3GFS/FV3GFS_io.F90 index 36e36bda..4d64b113 100644 --- a/FV3GFS/FV3GFS_io.F90 +++ b/FV3GFS/FV3GFS_io.F90 @@ -4585,7 +4585,7 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block Diag(idx)%axes = 2 Diag(idx)%name = 'rain' Diag(idx)%desc = 'total rain at this time step' - Diag(idx)%unit = 'XXX' + Diag(idx)%unit = 'm' Diag(idx)%mod_name = 'gfs_phys' allocate (Diag(idx)%data(nblks)) do nb = 1,nblks @@ -4596,13 +4596,36 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block Diag(idx)%axes = 2 Diag(idx)%name = 'rainc' Diag(idx)%desc = 'convective rain at this time step' - Diag(idx)%unit = 'XXX' + Diag(idx)%unit = 'm' Diag(idx)%mod_name = 'gfs_phys' allocate (Diag(idx)%data(nblks)) do nb = 1,nblks Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%rainc(:) enddo + + idx = idx + 1 + Diag(idx)%axes = 2 + Diag(idx)%name = 'precip_rate' + Diag(idx)%desc = 'instantaneous precipitation rate at this time step' + Diag(idx)%unit = 'kg/m2/s' + Diag(idx)%mod_name = 'gfs_phys' + allocate (Diag(idx)%data(nblks)) + do nb = 1,nblks + Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%rain(:)/Model%dtf * 1000. + enddo + + idx = idx + 1 + Diag(idx)%axes = 2 + Diag(idx)%name = 'precipc_rate' + Diag(idx)%desc = 'instantaneous convective precipitation rate at this time step' + Diag(idx)%unit = 'kg/m2/s' + Diag(idx)%mod_name = 'gfs_phys' + allocate (Diag(idx)%data(nblks)) + do nb = 1,nblks + Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%rainc(:)/Model%dtf * 1000. + enddo + idx = idx + 1 Diag(idx)%axes = 2 Diag(idx)%name = 'ice' From 741e4777aab996ee34ffb78200631c6855ebc584 Mon Sep 17 00:00:00 2001 From: Wolfgang Langhans Date: Mon, 8 Sep 2025 14:27:40 -0700 Subject: [PATCH 20/20] created new pointers that compute the converted rates --- FV3GFS/FV3GFS_io.F90 | 5 ++--- GFS_layer/GFS_physics_driver.F90 | 5 +++++ GFS_layer/GFS_typedefs.F90 | 6 ++++++ 3 files changed, 13 insertions(+), 3 deletions(-) diff --git a/FV3GFS/FV3GFS_io.F90 b/FV3GFS/FV3GFS_io.F90 index 4d64b113..ab0fdad2 100644 --- a/FV3GFS/FV3GFS_io.F90 +++ b/FV3GFS/FV3GFS_io.F90 @@ -4603,7 +4603,6 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%rainc(:) enddo - idx = idx + 1 Diag(idx)%axes = 2 Diag(idx)%name = 'precip_rate' @@ -4612,7 +4611,7 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block Diag(idx)%mod_name = 'gfs_phys' allocate (Diag(idx)%data(nblks)) do nb = 1,nblks - Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%rain(:)/Model%dtf * 1000. + Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%instant_precip_rate(:) enddo idx = idx + 1 @@ -4623,7 +4622,7 @@ subroutine gfdl_diag_register(Time, Sfcprop, Gfs_diag, Model, Cldprop, Atm_block Diag(idx)%mod_name = 'gfs_phys' allocate (Diag(idx)%data(nblks)) do nb = 1,nblks - Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%rainc(:)/Model%dtf * 1000. + Diag(idx)%data(nb)%var2 => Gfs_diag(nb)%instant_conv_precip_rate(:) enddo idx = idx + 1 diff --git a/GFS_layer/GFS_physics_driver.F90 b/GFS_layer/GFS_physics_driver.F90 index 31de03f8..a6039595 100644 --- a/GFS_layer/GFS_physics_driver.F90 +++ b/GFS_layer/GFS_physics_driver.F90 @@ -3971,6 +3971,11 @@ subroutine GFS_physics_driver & enddo endif + do i = 1, im + Diag%instant_precip_rate(i) = Diag%rain(i) / (dtf*con_p001) + Diag%instant_conv_precip_rate(i) = Diag%rainc(i) / (dtf*con_p001) + end do + deallocate (clw) deallocate (clw_trac_idx) if (Model%do_shoc) then diff --git a/GFS_layer/GFS_typedefs.F90 b/GFS_layer/GFS_typedefs.F90 index cfc92108..58cbf984 100644 --- a/GFS_layer/GFS_typedefs.F90 +++ b/GFS_layer/GFS_typedefs.F90 @@ -1324,6 +1324,8 @@ module GFS_typedefs real (kind=kind_phys), pointer :: wind10mmax(:) => null() !< maximum wind speed real (kind=kind_phys), pointer :: rain (:) => null() !< total rain at this time step real (kind=kind_phys), pointer :: rainc (:) => null() !< convective rain at this time step + real (kind=kind_phys), pointer :: instant_precip_rate(:) => null() !< total instant. precip rate (kg/m2/s) at this time step + real (kind=kind_phys), pointer :: instant_conv_precip_rate(:) => null() !< convective instant. precip rate (kg/m2/s) at this time step real (kind=kind_phys), pointer :: ice (:) => null() !< ice fall at this time step real (kind=kind_phys), pointer :: snow (:) => null() !< snow fall at this time step real (kind=kind_phys), pointer :: graupel(:) => null() !< graupel fall at this time step @@ -4100,6 +4102,8 @@ subroutine diag_create (Diag, IM, Model) allocate (Diag%gustmax100m (IM)) allocate (Diag%rain (IM)) allocate (Diag%rainc (IM)) + allocate (Diag%instant_precip_rate (IM)) + allocate (Diag%instant_conv_precip_rate (IM)) allocate (Diag%ice (IM)) allocate (Diag%snow (IM)) allocate (Diag%graupel (IM)) @@ -4411,6 +4415,8 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center) Diag%wind10mmax = zero Diag%rain = zero Diag%rainc = zero + Diag%instant_precip_rate = zero + Diag%instant_conv_precip_rate = zero Diag%ice = zero Diag%snow = zero Diag%graupel = zero