From 054413acda906cf88772ab84d0512508d0af439d Mon Sep 17 00:00:00 2001 From: barlage Date: Fri, 28 Apr 2023 10:30:45 -0600 Subject: [PATCH 01/16] document thermalz0 scheme --- physics/module_sf_noahmplsm.F90 | 273 ++++++++++++++++++-------------- 1 file changed, 158 insertions(+), 115 deletions(-) diff --git a/physics/module_sf_noahmplsm.F90 b/physics/module_sf_noahmplsm.F90 index ef6f99c44..ae143587f 100644 --- a/physics/module_sf_noahmplsm.F90 +++ b/physics/module_sf_noahmplsm.F90 @@ -5487,143 +5487,186 @@ end subroutine sfcdif3 !>\ingroup NoahMP_LSM ! compute thermal roughness length based on option opt_trs. - subroutine thermalz0(parameters,fveg,z0m,z0mg,zlvl,zpd,ezpd,ustarx, & !in - vegtyp,vaie,ur,csigmaf0,csigmaf1,aone,cdmnv,cdmng,icom, & !in - z0mt,z0ht) !out + + subroutine thermalz0(parameters, fveg, z0m, z0mg, zlvl, zpd, ezpd, & !in + ustarx, vegtyp, vaie, ur, c_sigma_f0, c_sigma_f1, a1, & !in + cdmn_v, cdmn_g, surface_flag, & !in + z0m_out, z0h_out ) !out + ! compute thermal roughness length based on option opt_trs. ! ------------------------------------------------------------------------------------------------- implicit none ! ------------------------------------------------------------------------------------------------- ! inputs - type (noahmp_parameters), intent(in) :: parameters !< - integer , intent(in ) :: vegtyp !< vegetation type - integer , intent(in ) :: icom !< 0=bared 1=vege 2=composition - real (kind=kind_phys), intent(in ) :: fveg !< green vegetation fraction [0.0-1.0] - real (kind=kind_phys), intent(in ) :: z0m !< z0 momentum (m) - real (kind=kind_phys), intent(in ) :: z0mg !< z0 momentum, ground (m) - real (kind=kind_phys), intent(in ) :: zlvl !< reference height [m] - real (kind=kind_phys), intent(in ) :: zpd !< zero plane displacement (m) - real (kind=kind_phys), intent(in ) :: ezpd !< zero plane displacement (m) - real (kind=kind_phys), intent(in ) :: ustarx !< friction velocity (m/s) - real (kind=kind_phys), intent(in ) :: vaie !< reference height [m] - real (kind=kind_phys), intent(in ) :: ur !< wind speed [m/s] - real (kind=kind_phys), intent(inout) :: csigmaf0 !< - real (kind=kind_phys), intent(inout) :: csigmaf1 !< - real (kind=kind_phys), intent(in ) :: aone !< - real (kind=kind_phys), intent(in ) :: cdmnv !< - real (kind=kind_phys), intent(in ) :: cdmng !< - real (kind=kind_phys), intent(out ) :: z0mt !< composited z0 momentum (m) - real (kind=kind_phys), intent(out ) :: z0ht !< composited z0 momentum (m) + type (noahmp_parameters),intent(in ) :: parameters ! parameters data structure + integer , intent(in ) :: vegtyp ! vegetation type + integer , intent(in ) :: surface_flag ! 0=bare 1=vegetation 2=composite + real (kind=kind_phys), intent(in ) :: fveg ! vegetation fraction [0.0-1.0] + real (kind=kind_phys), intent(in ) :: z0m ! z0 momentum [m] + real (kind=kind_phys), intent(in ) :: z0mg ! z0 momentum, ground [m] + real (kind=kind_phys), intent(in ) :: zlvl ! reference height [m] + real (kind=kind_phys), intent(in ) :: zpd ! zero plane displacement [m] + real (kind=kind_phys), intent(in ) :: ezpd ! grid zero plane displacement [m] + real (kind=kind_phys), intent(in ) :: ustarx ! friction velocity [m/s] + real (kind=kind_phys), intent(in ) :: vaie ! exposed LAI + SAI [m2/m2] + real (kind=kind_phys), intent(in ) :: ur ! wind speed [m/s] + real (kind=kind_phys), intent(in ) :: a1 ! Blumel 99 eqn 43 + real (kind=kind_phys), intent(in ) :: cdmn_v ! neutral momentum drag coefficient for vegetation + real (kind=kind_phys), intent(in ) :: cdmn_g ! neutral momentum drag coefficient for bare ground + real (kind=kind_phys), intent(inout) :: c_sigma_f0 ! C factor for no vegetation Blumel99 eqn 35 + real (kind=kind_phys), intent(inout) :: c_sigma_f1 ! C factor for full vegetation Blumel99 eqn 39 + real (kind=kind_phys), intent(out ) :: z0m_out ! output z0 momentum [m] + real (kind=kind_phys), intent(out ) :: z0h_out ! output z0 heat [m] ! local - real (kind=kind_phys) :: czil1 ! canopy based czil - real (kind=kind_phys) :: coeffa - real (kind=kind_phys) :: coeffb - real (kind=kind_phys) :: csigmafveg - real (kind=kind_phys) :: gsigma - real (kind=kind_phys) :: sigmaa - real (kind=kind_phys) :: cdmn - real (kind=kind_phys) :: kbsigmafveg - real (kind=kind_phys) :: reyn - real (kind=kind_phys) :: kbsigmaf0 - real (kind=kind_phys) :: kbsigmaf1 + real (kind=kind_phys) :: czil ! Zilitinkevich factor + real (kind=kind_phys) :: coeff_a ! slope of Blumel99 eqn 40 Blumel99 eqn 41 + real (kind=kind_phys) :: coeff_b ! intercept of Blumel99 eqn 40 Blumel99 eqn 42 + real (kind=kind_phys) :: c_sigma_fveg ! estimated C factor Blumel99 eqn 40 + real (kind=kind_phys) :: g_sigma ! weighting function Blumel99 eqn 22 + real (kind=kind_phys) :: sigma_a ! momentum partition factor Blumel99 eqn 8 + real (kind=kind_phys) :: cdmn ! grid neutral momentum drag coefficient Blumel99 eqn 21 + real (kind=kind_phys) :: reyn ! roughness Reynolds number Blumel99 eqn 36c + real (kind=kind_phys) :: kb_sigma_f0 ! bare ground kb^-1 Blumel99 eqn 36ab + real (kind=kind_phys) :: kb_sigma_f1 ! vegetated kb^-1 Blumel99 eqn 38 + real (kind=kind_phys) :: kb_sigma_fveg! grid estimated kb^-1 Blumel99 eqn 34 + + integer, parameter :: bare_flag = 0, vegetated_flag = 1, composite_flag = 2 + integer, parameter :: z0heqz0m = 1, & + chen09 = 2, & + tessel = 3, & + blumel99 = 4 + real (kind=kind_phys), parameter :: blumel_gamma = 0.5, & + blumel_zeta = 1.0, & + viscosity = 1.5e-5 ! ------------------------------------------------------------------------------------------------- - czil1 = 0.5 - coeffa = 0.0 - coeffb = 0.0 - csigmafveg= 0.0 - gsigma = 0.0 - cdmn = 0.0 - reyn = 0.0 - sigmaa = 0.0 - kbsigmafveg = 0.0 - kbsigmaf0 = 0.0 - kbsigmaf1 = 0.0 - if( icom == 2 )then - if (opt_trs == 1) then - z0mt = fveg * z0m + (1.0 - fveg) * z0mg - z0ht = z0mt - elseif (opt_trs == 2) then - z0mt = fveg * z0m + (1.0 - fveg) * z0mg - czil1=10.0 ** (- 0.4 * parameters%hvt) - z0ht = fveg * z0m*exp(-czil1*0.4*258.2*sqrt(ustarx*z0m)) & - +(1.0 - fveg) * z0mg*exp(-czil1*0.4*258.2*sqrt(ustarx*z0mg)) - elseif (opt_trs == 3) then - z0mt = fveg * z0m + (1.0 - fveg) * z0mg - if (vegtyp.le.5) then - z0ht = fveg * z0m + (1.0 - fveg) * z0mg*0.1 + czil = 0.5 + coeff_a = 0.0 + coeff_b = 0.0 + c_sigma_fveg = 0.0 + g_sigma = 0.0 + cdmn = 0.0 + reyn = 0.0 + sigma_a = 0.0 + kb_sigma_fveg = 0.0 + kb_sigma_f0 = 0.0 + kb_sigma_f1 = 0.0 + + surface_flag_select : select case(surface_flag) + + case (composite_flag) ! calculate grid based z0m and z0h + + if (opt_trs == z0heqz0m) then + + z0m_out = fveg * z0m + (1.0 - fveg) * z0mg ! probably should be log + z0h_out = z0m_out + + elseif (opt_trs == chen09) then + + z0m_out = fveg * z0m + (1.0 - fveg) * z0mg ! probably should be log + czil = 10.0 ** (- 0.4 * parameters%hvt) + z0h_out = fveg * z0m * exp(-czil*0.4*258.2*sqrt(ustarx*z0m )) & + + (1.0 - fveg) * z0mg * exp(-czil*0.4*258.2*sqrt(ustarx*z0mg)) + + elseif (opt_trs == tessel) then + + z0m_out = fveg * z0m + (1.0 - fveg) * z0mg ! probably should be log + if (vegtyp <= 5) then + z0h_out = fveg * z0m + (1.0 - fveg) * z0mg * 0.1 else - z0ht = fveg * z0m*0.01 + (1.0 - fveg) * z0mg*0.1 + z0h_out = fveg * z0m * 0.01 + (1.0 - fveg) * z0mg * 0.1 endif - elseif (opt_trs == 4) then - coeffa = (csigmaf0 - csigmaf1)/(1.0 - exp(-1.0*aone)) - coeffb = csigmaf0 - coeffa - csigmafveg = coeffa * exp(-1.0*aone*fveg) + coeffb - gsigma = fveg**0.5 + fveg*(1.0-fveg)*1.0 -! -! 0.5 ~ 1.0 for the 0.5 place; 0 ~ 1.0 for the 1.0 place, adjustable empirical + elseif (opt_trs == blumel99) then + + coeff_a = (c_sigma_f0 - c_sigma_f1)/(1.0 - exp(-1.0*a1)) ! Blumel99 eqn 41 + coeff_b = c_sigma_f0 - coeff_a ! Blumel99 eqn 42 + c_sigma_fveg = coeff_a * exp(-1.0*a1*fveg) + coeff_b ! Blumel99 eqn 40 + +! blumel_gamma = 0.5 ~ 1.0 and blumel_zeta = 0 ~ 1.0, adjustable empirical ! canopy roughness geometry parameter; currently fveg = 0.78 has the largest ! momentum flux; can test the fveg-based average by setting 0.5 to 1.0 and 1.0 -! to 0.0 ! see Blumel; JAM,1998 -! +! to 0.0 ! see Blumel; JAM,1999 - cdmn = gsigma*cdmnv + (1.0-gsigma)*cdmng - z0mt = (zlvl - ezpd)*exp(-0.4/sqrt(cdmn)) + g_sigma = fveg**blumel_gamma + fveg*(1.0-fveg)*blumel_zeta ! Blumel99 eqn 22 + cdmn = g_sigma*cdmn_v + (1.0-g_sigma)*cdmn_g ! Blumel99 eqn 21 + z0m_out = (zlvl - ezpd)*exp(-0.4/sqrt(cdmn)) ! Blumel99 eqn 24 + kb_sigma_fveg = c_sigma_fveg/log((zlvl-ezpd)/z0m_out) - & + log((zlvl-ezpd)/z0m_out) ! Blumel99 eqn 34 + z0h_out = z0m_out/exp(kb_sigma_fveg) - kbsigmafveg = csigmafveg/log((zlvl-ezpd)/z0mt) - log((zlvl-ezpd)/z0mt) - z0ht = z0mt/exp(kbsigmafveg) - endif - - elseif( icom == 0 )then - - z0mt = z0mg - if (opt_trs == 1) then - z0ht = z0mt - elseif (opt_trs == 2) then - czil1=10.0 ** (- 0.4 * parameters%hvt) - z0ht =z0mt*exp(-czil1*0.4*258.2*sqrt(ustarx*z0mt)) - elseif (opt_trs == 3) then - if (vegtyp.le.5) then - z0ht = z0mt - else - z0ht = z0mt*0.01 endif - elseif (opt_trs == 4) then - reyn = ustarx*z0mt/(1.5e-05) - if (reyn .gt. 2.0) then - kbsigmaf0 = 2.46*reyn**0.25 - log(7.4) - else - kbsigmaf0 = - log(0.397) + + case (bare_flag) ! calculate z0m and z0h over bare tile + + z0m_out = z0mg + + if (opt_trs == z0heqz0m) then + + z0h_out = z0m_out + + elseif (opt_trs == chen09) then + + czil = 10.0 ** (- 0.4 * parameters%hvt) + z0h_out = z0m_out * exp(-czil*0.4*258.2*sqrt(ustarx*z0m_out)) + + elseif (opt_trs == tessel) then + + if (vegtyp <= 5) then + z0h_out = z0m_out + else + z0h_out = z0m_out * 0.01 + endif + + elseif (opt_trs == blumel99) then + + reyn = ustarx*z0m_out/viscosity ! Blumel99 eqn 36c + if (reyn > 2.0) then + kb_sigma_f0 = 2.46*reyn**0.25 - log(7.4) ! Blumel99 eqn 36a + else + kb_sigma_f0 = - log(0.397) ! Blumel99 eqn 36b + endif + + z0h_out = max(z0m_out/exp(kb_sigma_f0),1.0e-6) + c_sigma_f0 = log((zlvl-zpd)/z0m_out) * & + (log((zlvl-zpd)/z0m_out) + kb_sigma_f0) ! Blumel99 eqn 35 + endif - z0ht = max(z0mt/exp(kbsigmaf0),1.0e-6) - csigmaf0 = log((zlvl-zpd)/z0mt)*(log((zlvl-zpd)/z0mt) + kbsigmaf0) - endif + case (vegetated_flag) ! calculate z0m and z0h over vegetated tile - elseif( icom == 1 )then - - z0mt = z0m - if (opt_trs == 1) then - z0ht = z0mt - elseif (opt_trs == 2) then - czil1= 10.0 ** (- 0.4 * parameters%hvt) - z0ht = z0mt*exp(-czil1*0.4*258.2*sqrt(ustarx*z0mt)) - elseif (opt_trs == 3) then - if (vegtyp.le.5) then - z0ht = z0mt - else - z0ht = z0mt*0.01 - endif - elseif (opt_trs == 4) then - sigmaa = 1.0 - (0.5/(0.5+vaie))*exp(-vaie**2/8.0) - kbsigmaf1 = 16.4*(sigmaa*vaie**3)**(-0.25)*sqrt(parameters%dleaf*ur/log((zlvl-zpd)/z0mt)) - z0ht = z0mt/exp(kbsigmaf1) - csigmaf1 = log((zlvl-zpd)/z0mt)*(log((zlvl-zpd)/z0mt)+kbsigmaf1) ! for output for interpolation + z0m_out = z0m + + if (opt_trs == z0heqz0m) then + + z0h_out = z0m_out + + elseif (opt_trs == chen09) then + + czil = 10.0 ** (- 0.4 * parameters%hvt) + z0h_out = z0m_out * exp(-czil*0.4*258.2*sqrt(ustarx*z0m_out)) + + elseif (opt_trs == tessel) then + + if (vegtyp <= 5) then + z0h_out = z0m_out + else + z0h_out = z0m_out*0.01 endif - endif + + elseif (opt_trs == blumel99) then + + sigma_a = 1.0 - (0.5/(0.5+vaie)) * exp(-vaie**2/8.0) ! Blumel99 eqn 8 + kb_sigma_f1 = 16.4 * (sigma_a*vaie**3)**(-0.25) * & ! Blumel99 eqn 38 + sqrt(parameters%dleaf*ur/log((zlvl-zpd)/z0m_out)) + z0h_out = z0m_out/exp(kb_sigma_f1) + c_sigma_f1 = log((zlvl-zpd)/z0m_out)*(log((zlvl-zpd)/z0m_out)+kb_sigma_f1) ! Blumel99 eqn 39 + + endif + + end select surface_flag_select end subroutine thermalz0 From c39a6a50cb8791d0678a7e6bd3aa111f2ec75c32 Mon Sep 17 00:00:00 2001 From: barlage Date: Mon, 1 May 2023 09:09:04 -0600 Subject: [PATCH 02/16] add z0m option to account for phenology --- physics/module_sf_noahmplsm.F90 | 57 +++++++++++++++++++++++++++------ physics/noahmpdrv.F90 | 4 ++- 2 files changed, 50 insertions(+), 11 deletions(-) diff --git a/physics/module_sf_noahmplsm.F90 b/physics/module_sf_noahmplsm.F90 index ae143587f..afea081f3 100644 --- a/physics/module_sf_noahmplsm.F90 +++ b/physics/module_sf_noahmplsm.F90 @@ -163,14 +163,19 @@ module module_sf_noahmplsm ! 1 -> liu, et al. 2016 integer :: opt_trs !< options for thermal roughness scheme - ! **1 -> z0h=z0 - ! 2 -> czil - ! 3 -> ec style - ! 4 -> kb inversed + ! **1 -> z0h=z0m + ! 2 -> czil = f(canopy height) from Chen09 + ! 3 -> ec style from TESSEL + ! 4 -> kb inverse from Blumel99 integer :: opt_diag !< options for surface 2m/q diagnostic approach ! 1 -> external GFS sfc_diag ! **2 -> original NoahMP 2-title ! 3 -> NoahMP 2-title + internal GFS sfc_diag + + integer :: opt_z0m !< options for momentum roughness length + ! **1 -> use z0m from MPTABLE + ! 2 -> z0m = f(canopy height, LAI/SAI) + !------------------------------------------------------------------------------------------! ! physical constants: ! !------------------------------------------------------------------------------------------! @@ -1974,6 +1979,10 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in real (kind=kind_phys) :: ezpd real (kind=kind_phys) :: aone + real (kind=kind_phys) :: canopy_density_factor + real (kind=kind_phys) :: vai_limited + real (kind=kind_phys) :: z0m_hvt_ratio(20) + !jref:end real (kind=kind_phys), parameter :: mpe = 1.e-6 @@ -2012,6 +2021,11 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in csigmaf1 = 0.0 csigmaf0 = 0.0 aone = 0.0 + + canopy_density_factor = 1.0 + vai_limited = 2.0 + z0m_hvt_ratio = (/ 0.545,0.055,0.047,0.050,0.050,0.182,0.545,0.046,0.050,0.120, & + 0.060,0.075,0.067,0.093,0.000,0.000,0.000,0.075,0.100,0.060 /) ! @@ -2054,12 +2068,32 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in zpdg = snowh if(veg) then - z0m = parameters%z0mvt - zpd = 0.65 * parameters%hvt - if(snowh.gt.zpd) zpd = snowh + + if(opt_z0m == 1) then + + z0m = parameters%z0mvt + zpd = 0.65 * parameters%hvt + + elseif(opt_z0m == 2) then + + z0m = z0m_hvt_ratio(vegtyp) * parameters%hvt + zpd = 0.65 * parameters%hvt + if(vegtyp /= 13) then + vai_limited = min(vai, 2.0) + canopy_density_factor = (1.0 - exp(-vai_limited)) / (1.0 - exp(-2.0)) + z0m = exp(canopy_density_factor * log(z0m) + (1.0 - canopy_density_factor) * log(z0mg)) + zpd = canopy_density_factor * zpd + end if + + end if + + if(snowh.gt.zpd) zpd = snowh + else + z0m = z0mg zpd = zpdg + end if ! special case for urban @@ -10154,9 +10188,10 @@ end subroutine psn_crop !>\ingroup NoahMP_LSM !! - subroutine noahmp_options(idveg ,iopt_crs ,iopt_btr ,iopt_run ,iopt_sfc ,iopt_frz , & - iopt_inf ,iopt_rad ,iopt_alb ,iopt_snf ,iopt_tbot, iopt_stc, & - iopt_rsf , iopt_soil, iopt_pedo, iopt_crop ,iopt_trs, iopt_diag ) + subroutine noahmp_options(idveg , iopt_crs , iopt_btr , iopt_run , iopt_sfc , iopt_frz , & + iopt_inf, iopt_rad , iopt_alb , iopt_snf , iopt_tbot, iopt_stc , & + iopt_rsf, iopt_soil, iopt_pedo, iopt_crop, iopt_trs , iopt_diag, & + iopt_z0m ) implicit none @@ -10180,6 +10215,7 @@ subroutine noahmp_options(idveg ,iopt_crs ,iopt_btr ,iopt_run ,iopt_sfc integer, intent(in) :: iopt_crop !< crop model option (0->none; 1->liu et al.) integer, intent(in) :: iopt_trs !< thermal roughness scheme option (1->z0h=z0; 2->rb reversed) integer, intent(in) :: iopt_diag !< surface 2m t/q diagnostic approach + integer, intent(in) :: iopt_z0m !< momentum roughness length option ! ------------------------------------------------------------------------------------------------- @@ -10202,6 +10238,7 @@ subroutine noahmp_options(idveg ,iopt_crs ,iopt_btr ,iopt_run ,iopt_sfc opt_crop = iopt_crop opt_trs = iopt_trs opt_diag = iopt_diag + opt_z0m = iopt_z0m end subroutine noahmp_options diff --git a/physics/noahmpdrv.F90 b/physics/noahmpdrv.F90 index 88e75637b..c1200e829 100644 --- a/physics/noahmpdrv.F90 +++ b/physics/noahmpdrv.F90 @@ -451,6 +451,7 @@ subroutine noahmpdrv_run & integer :: iopt_pedo = 1 ! option for pedotransfer function integer :: iopt_crop = 0 ! option for crop model integer :: iopt_gla = 2 ! option for glacier treatment + integer :: iopt_z0m = 2 ! option for z0m treatment ! ! --- local inputs to noah-mp and glacier subroutines; listed in order in noah-mp call @@ -848,7 +849,8 @@ subroutine noahmpdrv_run & call noahmp_options(idveg ,iopt_crs, iopt_btr , iopt_run, iopt_sfc, & iopt_frz, iopt_inf , iopt_rad, iopt_alb, & iopt_snf, iopt_tbot, iopt_stc, iopt_rsf, & - iopt_soil,iopt_pedo, iopt_crop,iopt_trs,iopt_diag) + iopt_soil,iopt_pedo, iopt_crop,iopt_trs, & + iopt_diag,iopt_z0m) if ( vegetation_category == isice_table ) then From 3ca1cb482ff14b5d85ecb9c931ad21907380e192 Mon Sep 17 00:00:00 2001 From: Michael Barlage Date: Mon, 1 May 2023 11:46:36 -0400 Subject: [PATCH 03/16] revert z0m option to original for testing --- physics/noahmpdrv.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/noahmpdrv.F90 b/physics/noahmpdrv.F90 index c1200e829..fb1859cc9 100644 --- a/physics/noahmpdrv.F90 +++ b/physics/noahmpdrv.F90 @@ -451,7 +451,7 @@ subroutine noahmpdrv_run & integer :: iopt_pedo = 1 ! option for pedotransfer function integer :: iopt_crop = 0 ! option for crop model integer :: iopt_gla = 2 ! option for glacier treatment - integer :: iopt_z0m = 2 ! option for z0m treatment + integer :: iopt_z0m = 1 ! option for z0m treatment ! ! --- local inputs to noah-mp and glacier subroutines; listed in order in noah-mp call From c92f164bad27bd6678675981ff569bf563334456 Mon Sep 17 00:00:00 2001 From: Michael Barlage Date: Mon, 1 May 2023 16:37:35 -0400 Subject: [PATCH 04/16] add gfs stability inside noahmp; remove dependence on sfc_diff --- physics/module_sf_noahmplsm.F90 | 226 +++++++++++++++++++++++++++++++- physics/noahmpdrv.F90 | 5 +- 2 files changed, 226 insertions(+), 5 deletions(-) diff --git a/physics/module_sf_noahmplsm.F90 b/physics/module_sf_noahmplsm.F90 index afea081f3..33150852a 100644 --- a/physics/module_sf_noahmplsm.F90 +++ b/physics/module_sf_noahmplsm.F90 @@ -8,7 +8,6 @@ module module_sf_noahmplsm use module_wrf_utl #endif use machine , only : kind_phys -use sfc_diff, only : stability implicit none @@ -5511,11 +5510,234 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur tvs = tgb/prsik1x * virtfac endif - call stability (zlvlb, zvfun1, gdx, tv1, thv1, ur, z0m, z0h, tvs, grav, thsfc_loc, & + call gfs_stability (zlvlb, zvfun1, gdx, tv1, thv1, ur, z0m, z0h, tvs, grav, thsfc_loc, & rb1, fm,fh,fm10,fh2,cm,ch,stress1,fv) end subroutine sfcdif3 +!== begin gfs_stability ================================================================================== + +subroutine gfs_stability & +! --- inputs: + ( z1, zvfun, gdx, tv1, thv1, wind, z0max, ztmax, tvs, grav, & + thsfc_loc, & +! --- outputs: + rb, fm, fh, fm10, fh2, cm, ch, stress, ustar) + +! Documentation below refers to UTN and STN which are: +! UTN (Unstable Tech Note) : NCEP Office Note 356 +! STN (Stable Tech Note) : NCEP Office Note 321 + +integer, parameter :: kp = kind_phys +real (kind=kind_phys), parameter :: ca=0.4_kind_phys ! ca - von karman constant + +real(kind=kind_phys), intent(in) :: z1 ! height model level +real(kind=kind_phys), intent(in) :: zvfun ! vegetation adjustment factor +real(kind=kind_phys), intent(in) :: gdx ! grid spatial dimension +real(kind=kind_phys), intent(in) :: tv1 ! virtual temperature at model level +real(kind=kind_phys), intent(in) :: thv1 ! virtual potential temperature at model level +real(kind=kind_phys), intent(in) :: wind ! wind speed at model level +real(kind=kind_phys), intent(in) :: z0max ! momentum roughness length +real(kind=kind_phys), intent(in) :: ztmax ! thermal roughness length +real(kind=kind_phys), intent(in) :: tvs ! surface virtual temperature +real(kind=kind_phys), intent(in) :: grav ! local gravity +logical, intent(in) :: thsfc_loc ! use local theta reference flag + +real(kind=kind_phys), intent(out) :: rb ! bulk richardson number [-] +real(kind=kind_phys), intent(out) :: fm ! phi momentum function (UTN 1.1) [-] +real(kind=kind_phys), intent(out) :: fh ! phi heat function (UTN 1.2) [-] +real(kind=kind_phys), intent(out) :: fm10 ! 10-meter phi momentum function [-] +real(kind=kind_phys), intent(out) :: fh2 ! 2-meter phi heat function [-] +real(kind=kind_phys), intent(out) :: cm ! momentum exchange coeficient [-] +real(kind=kind_phys), intent(out) :: ch ! heat exchange coeficient [-] +real(kind=kind_phys), intent(out) :: stress ! surface stress [m2/s2] +real(kind=kind_phys), intent(out) :: ustar ! friction velocity [m/s] + +! --- locals: +real(kind=kind_phys), parameter :: a0 = -3.975 ! UTN 2.37 +real(kind=kind_phys), parameter :: a1 = 12.32 ! UTN 2.37 +real(kind=kind_phys), parameter :: b1 = -7.755 ! UTN 2.37 +real(kind=kind_phys), parameter :: b2 = 6.041 ! UTN 2.37 +real(kind=kind_phys), parameter :: a0p = -7.941 ! UTN 2.38 +real(kind=kind_phys), parameter :: a1p = 24.75 ! UTN 2.38 +real(kind=kind_phys), parameter :: b1p = -8.705 ! UTN 2.38 +real(kind=kind_phys), parameter :: b2p = 7.899 ! UTN 2.38 + +real(kind=kind_phys), parameter :: alpha = 5.0 ! alpha in e.g., STN 1.10 +real(kind=kind_phys), parameter :: alpha4 = 4.0 * alpha ! term in aa +real(kind=kind_phys), parameter :: xkrefsqr = 0.3 ! baseline maximum z/L +real(kind=kind_phys), parameter :: xkmin = 0.05 ! min multiplier for grid size and vegetation +real(kind=kind_phys), parameter :: xkgdx = 3000.0 ! critical grid scale for diffusivity[m^0.5] +real(kind=kind_phys), parameter :: zolmin = -10.0 ! minimum z/L +real(kind=kind_phys), parameter :: zero = 0.0 +real(kind=kind_phys), parameter :: one = 1.0 + +real(kind=kind_phys) :: aa +real(kind=kind_phys) :: aa0 +real(kind=kind_phys) :: bb +real(kind=kind_phys) :: bb0 +real(kind=kind_phys) :: dtv +real(kind=kind_phys) :: adtv +real(kind=kind_phys) :: hl1 +real(kind=kind_phys) :: hl12 +real(kind=kind_phys) :: pm +real(kind=kind_phys) :: ph +real(kind=kind_phys) :: pm10 +real(kind=kind_phys) :: ph2 +real(kind=kind_phys) :: z1i +real(kind=kind_phys) :: fms +real(kind=kind_phys) :: fhs +real(kind=kind_phys) :: hl0 +real(kind=kind_phys) :: hl0inf +real(kind=kind_phys) :: hlinf +real(kind=kind_phys) :: hl110 +real(kind=kind_phys) :: hlt +real(kind=kind_phys) :: hltinf +real(kind=kind_phys) :: olinf +real(kind=kind_phys) :: tem1 +real(kind=kind_phys) :: tem2 +real(kind=kind_phys) :: zolmax + +real(kind=kind_phys) xkzo + +z1i = one / z1 ! inverse of model height + +! +! set background diffusivities with one for gdx >= xkgdx and +! as a function of horizontal grid size for gdx < xkgdx +! (i.e., gdx/xkgdx for gdx < xkgdx) +! + +if(gdx >= xkgdx) then + xkzo = one +else + xkzo = gdx / xkgdx +endif + +tem1 = tv1 - tvs +if(tem1 > zero) then ! for stable case, adjust for vegetation cover + tem2 = xkzo * zvfun + xkzo = min(max(tem2, xkmin), xkzo) +endif + +zolmax = xkrefsqr / sqrt(xkzo) ! maximum z/L + +! compute stability indices (rb and hlinf) + + dtv = thv1 - tvs + adtv = max(abs(dtv),0.001_kp) + dtv = sign(1.0_kp,dtv) * adtv + + if(thsfc_loc) then ! Use local potential temperature + rb = max(-5000.0_kp, (grav+grav) * dtv * z1 & + / ((thv1 + tvs) * wind * wind)) + else ! Use potential temperature referenced to 1000 hPa + rb = max(-5000.0_kp, grav * dtv * z1 & + / (tv1 * wind * wind)) + endif + + tem1 = one / z0max ! 1/z0m + tem2 = one / ztmax ! 1/z0t + fm = log((z0max+z1) * tem1) ! neutral phi_m + fh = log((ztmax+z1) * tem2) ! neutral phi_h + fm10 = log((z0max+10.0_kp) * tem1) ! neutral phi_m at 10 meters + fh2 = log((ztmax+2.0_kp) * tem2) ! neutral phi_h at 2 meters + hlinf = rb * fm * fm / fh ! z/L STN 2.7 + hlinf = min(max(hlinf,zolmin),zolmax) ! z/L, xi in STN/UTN +! +! stable case +! + if (dtv >= zero) then + hl1 = hlinf ! z/L, xi in STN + if(hlinf > 0.25_kp) then ! z/L > 0.25, do two iterations + tem1 = hlinf * z1i ! 1/L + hl0inf = z0max * tem1 ! z0m/z1, zi_0 in STN + hltinf = ztmax * tem1 ! z0t/z1, zi_0 in STN + aa = sqrt(one + alpha4 * hlinf) ! sqrt term of STN 2.16 with z + aa0 = sqrt(one + alpha4 * hl0inf) ! sqrt term of STN 2.16 with z0m + bb = aa ! sqrt term of STN 2.16 with z + bb0 = sqrt(one + alpha4 * hltinf) ! sqrt term of STN 2.16 with z0t + pm = aa0 - aa + log( (aa + one)/(aa0 + one) ) ! psi_m STN 3.11 + ph = bb0 - bb + log( (bb + one)/(bb0 + one) ) ! psi_h STN 3.11 + fms = fm - pm ! phi_m STN 3.10 + fhs = fh - ph ! phi_h STN 3.10 + hl1 = fms * fms * rb / fhs ! z/L iteration STN 3.8 + hl1 = min(hl1, zolmax) ! z/L iteration + endif +! +! second iteration +! + tem1 = hl1 * z1i ! 1/L + hl0 = z0max * tem1 ! z0m/z1 + hlt = ztmax * tem1 ! z0t/z1 + aa = sqrt(one + alpha4 * hl1) ! sqrt term of STN 2.16 with z + aa0 = sqrt(one + alpha4 * hl0) ! sqrt term of STN 2.16 with z0m + bb = aa ! sqrt term of STN 2.16 with z + bb0 = sqrt(one + alpha4 * hlt) ! sqrt term of STN 2.16 with z0t + pm = aa0 - aa + log( (one+aa)/(one+aa0) ) ! psi_m STN 3.11 + ph = bb0 - bb + log( (one+bb)/(one+bb0) ) ! psi_h STN 3.11 + hl110 = hl1 * 10.0_kp * z1i ! 10/L + aa = sqrt(one + alpha4 * hl110) ! sqrt term of STN 2.16 with z=10m + pm10 = aa0 - aa + log( (one+aa)/(one+aa0) ) ! psi_m STN 3.11 with z=10m + hl12 = (hl1+hl1) * z1i ! 2/L +! aa = sqrt(one + alpha4 * hl12) + bb = sqrt(one + alpha4 * hl12) ! sqrt term of STN 2.16 with z=2m + ph2 = bb0 - bb + log( (one+bb)/(one+bb0) ) ! psi_m STN 3.11 with z=2m +! +! unstable case - check for unphysical obukhov length +! see steps in UTN Sec. D +! + else ! dtv < 0 case + + olinf = z1 / hlinf ! z/L, xi in UTN + tem1 = 50.0_kp * z0max ! 50 * z0m, z/L limit for calc methods, see UTN Sec. E + if(abs(olinf) <= tem1) then ! + hlinf = -z1 / tem1 ! + hlinf = max(hlinf, zolmin) + endif +! +! get pm and ph +! + if (hlinf >= -0.5_kp) then + hl1 = hlinf + pm = (a0 + a1*hl1) * hl1 / (one+ (b1+b2*hl1) *hl1) ! psi_m UTN 2.37 + ph = (a0p + a1p*hl1) * hl1 / (one+ (b1p+b2p*hl1)*hl1) ! psi_h UTN 2.38 + hl110 = hl1 * 10.0_kp * z1i ! 10/L + pm10 = (a0 + a1*hl110) * hl110/(one+(b1+b2*hl110)*hl110) ! psi_m UTN 2.37 with z=10m + hl12 = (hl1+hl1) * z1i ! 2/L + ph2 = (a0p + a1p*hl12) * hl12/(one+(b1p+b2p*hl12)*hl12) ! psi_h UTN 2.38 with z=2m + else ! z/L < -0.5 + hl1 = -hlinf ! -z/L + tem1 = one / sqrt(hl1) ! sqrt(-z/L) + pm = log(hl1) + 2.0_kp * sqrt(tem1) - 0.8776_kp ! UTN 2.64, first three terms + ph = log(hl1) + 0.5_kp * tem1 + 1.386_kp ! UTN 2.65, first three terms + hl110 = hl1 * 10.0_kp * z1i ! 10/L + pm10 = log(hl110) + 2.0_kp/sqrt(sqrt(hl110)) - 0.8776_kp ! psi_m UTN 2.64 with z=10m + hl12 = (hl1+hl1) * z1i ! 2/L + ph2 = log(hl12) + 0.5_kp / sqrt(hl12) + 1.386_kp ! psi_h UTN 2.65 with z=2m + endif + + endif ! end of if (dtv >= 0 ) then loop +! +! finish the exchange coefficient computation to provide fm and fh +! + fm = fm - pm ! phi_m + fh = fh - ph ! phi_h + fm10 = fm10 - pm10 ! phi_m at 10m + fh2 = fh2 - ph2 ! phi_h at 2m + cm = ca * ca / (fm * fm) ! momentum exchange coef = k^2/phi_m^2 + ch = ca * ca / (fm * fh) ! heat exchange coef = k^2/phi_m/phi_h + tem1 = 0.00001_kp/z1 ! minimum exhange coef (?) + cm = max(cm, tem1) + ch = max(ch, tem1) + stress = cm * wind * wind ! surface stress = Cm*U*U + ustar = sqrt(stress) ! friction velocity + + return +!................................. + end subroutine gfs_stability +!--------------------------------- + !== begin thermalz0 !================================================================================== diff --git a/physics/noahmpdrv.F90 b/physics/noahmpdrv.F90 index fb1859cc9..d232b759f 100644 --- a/physics/noahmpdrv.F90 +++ b/physics/noahmpdrv.F90 @@ -203,8 +203,7 @@ subroutine noahmpdrv_run & use machine , only : kind_phys use funcphys, only : fpvs - use sfc_diff, only : stability -! use module_sf_noahmplsm + use module_sf_noahmplsm, only : gfs_stability use module_sf_noahmp_glacier use noahmp_tables @@ -1189,7 +1188,7 @@ subroutine noahmpdrv_run & ! if ( .not. do_mynnsfclay) then !GFS sfcdiff if ( iopt_sfc .ne. 4 ) then !GFS sfcdiff - call stability & + call gfs_stability & (zf(i), zvfun(i), gdx, virtual_temperature, vptemp,wind(i), z0_total, z0h_total, & tvs1, con_g, thsfc_loc, & rb1(i), fm1(i), fh1(i), fm101(i), fh21(i), cm(i), ch(i), stress1(i), ustar1(i)) From 5c438fd3622b42a5f621f89497d8bc054264a402 Mon Sep 17 00:00:00 2001 From: Michael Barlage Date: Mon, 1 May 2023 17:22:11 -0400 Subject: [PATCH 05/16] add z0m to hvt ratio to mptable --- physics/module_sf_noahmplsm.F90 | 6 ++---- physics/noahmp_tables.f90 | 12 +++++++++--- physics/noahmpdrv.F90 | 3 ++- physics/noahmptable.tbl | 2 ++ 4 files changed, 15 insertions(+), 8 deletions(-) diff --git a/physics/module_sf_noahmplsm.F90 b/physics/module_sf_noahmplsm.F90 index 33150852a..b6e751e31 100644 --- a/physics/module_sf_noahmplsm.F90 +++ b/physics/module_sf_noahmplsm.F90 @@ -219,6 +219,7 @@ module module_sf_noahmplsm real (kind=kind_phys) :: z0mvt !< momentum roughness length (m) real (kind=kind_phys) :: hvt !< top of canopy (m) real (kind=kind_phys) :: hvb !< bottom of canopy (m) + real (kind=kind_phys) :: z0mhvt !< ratio of z0m to hvt real (kind=kind_phys) :: den !< tree density (no. of trunks per m2) real (kind=kind_phys) :: rc !< tree crown radius (m) real (kind=kind_phys) :: mfsno !< snowmelt m parameter () @@ -1980,7 +1981,6 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in real (kind=kind_phys) :: canopy_density_factor real (kind=kind_phys) :: vai_limited - real (kind=kind_phys) :: z0m_hvt_ratio(20) !jref:end @@ -2023,8 +2023,6 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in canopy_density_factor = 1.0 vai_limited = 2.0 - z0m_hvt_ratio = (/ 0.545,0.055,0.047,0.050,0.050,0.182,0.545,0.046,0.050,0.120, & - 0.060,0.075,0.067,0.093,0.000,0.000,0.000,0.075,0.100,0.060 /) ! @@ -2075,7 +2073,7 @@ subroutine energy (parameters,ice ,vegtyp ,ist ,nsnow ,nsoil , & !in elseif(opt_z0m == 2) then - z0m = z0m_hvt_ratio(vegtyp) * parameters%hvt + z0m = parameters%z0mhvt * parameters%hvt zpd = 0.65 * parameters%hvt if(vegtyp /= 13) then vai_limited = min(vai, 2.0) diff --git a/physics/noahmp_tables.f90 b/physics/noahmp_tables.f90 index 0e44b3cfc..cc7856af6 100644 --- a/physics/noahmp_tables.f90 +++ b/physics/noahmp_tables.f90 @@ -44,6 +44,7 @@ module noahmp_tables real :: z0mvt_table(mvt) !< momentum roughness length (m) real :: hvt_table(mvt) !< top of canopy (m) real :: hvb_table(mvt) !< bottom of canopy (m) + real :: z0mhvt_table(mvt) !< ratio of z0m to hvt real :: den_table(mvt) !< tree density (no. of trunks per m2) real :: rc_table(mvt) !< tree crown radius (m) real :: mfsno_table(mvt) !< snowmelt curve parameter () @@ -323,7 +324,8 @@ subroutine read_mp_table_parameters sai_sep, sai_oct, sai_nov, sai_dec, lai_jan, lai_feb, lai_mar, lai_apr, & lai_may, lai_jun, lai_jul, lai_aug, lai_sep, lai_oct, lai_nov, lai_dec, & rhol_vis, rhol_nir, rhos_vis, rhos_nir, taul_vis, taul_nir, taus_vis, taus_nir,& - ch2op, dleaf, z0mvt, hvt, hvb, den, rc, mfsno, scffac, xl, cwpvt, c3psn, kc25, & + ch2op, dleaf, z0mvt, hvt, hvb, z0mhvt, & + den, rc, mfsno, scffac, xl, cwpvt, c3psn, kc25, & akc, ko25, ako, avcmx, aqe, ltovrc, dilefc, dilefw, rmf25, sla, fragr, tmin, & vcmx25, tdlef, bp, mp, qe25, rms25, rmr25, arm, folnmx, wdpool, wrrat, mrp, & nroot, rgl, rs, hs, topt, rsmax, rtovrc, rswoodc, bf, wstrc, laimin, & @@ -331,7 +333,8 @@ subroutine read_mp_table_parameters namelist / noahmp_usgs_veg_categories / veg_dataset_description, nveg namelist / noahmp_usgs_parameters / isurban, iswater, isbarren, isice, iscrop, eblforest, natural, & lcz_1, lcz_2, lcz_3, lcz_4, lcz_5, lcz_6, lcz_7, lcz_8, lcz_9, lcz_10, lcz_11, & - ch2op, dleaf, z0mvt, hvt, hvb, den, rc, mfsno, scffac, xl, cwpvt, c3psn, kc25, & + ch2op, dleaf, z0mvt, hvt, hvb, z0mhvt, & + den, rc, mfsno, scffac, xl, cwpvt, c3psn, kc25, & akc, ko25, ako, avcmx, aqe, ltovrc, dilefc, dilefw, rmf25, sla, fragr, tmin, & vcmx25, tdlef, bp, mp, qe25, rms25, rmr25, arm, folnmx, wdpool, wrrat, mrp, & nroot, rgl, rs, hs, topt, rsmax, rtovrc, rswoodc, bf, wstrc, laimin, & @@ -343,7 +346,8 @@ subroutine read_mp_table_parameters namelist / noahmp_modis_veg_categories / veg_dataset_description, nveg namelist / noahmp_modis_parameters / isurban, iswater, isbarren, isice, iscrop, eblforest, natural, & lcz_1, lcz_2, lcz_3, lcz_4, lcz_5, lcz_6, lcz_7, lcz_8, lcz_9, lcz_10, lcz_11, & - ch2op, dleaf, z0mvt, hvt, hvb, den, rc, mfsno, scffac, xl, cwpvt, c3psn, kc25, & + ch2op, dleaf, z0mvt, hvt, hvb, z0mhvt, & + den, rc, mfsno, scffac, xl, cwpvt, c3psn, kc25, & akc, ko25, ako, avcmx, aqe, ltovrc, dilefc, dilefw, rmf25, sla, fragr, tmin, & vcmx25, tdlef, bp, mp, qe25, rms25, rmr25, arm, folnmx, wdpool, wrrat, mrp, & nroot, rgl, rs, hs, topt, rsmax, rtovrc, rswoodc, bf, wstrc, laimin, & @@ -502,6 +506,7 @@ subroutine read_mp_table_parameters z0mvt_table = -1.0e36 hvt_table = -1.0e36 hvb_table = -1.0e36 + z0mhvt_table = -1.0e36 den_table = -1.0e36 rc_table = -1.0e36 mfsno_table = -1.0e36 @@ -814,6 +819,7 @@ subroutine read_mp_table_parameters z0mvt_table (1:nveg) = z0mvt (1:nveg) hvt_table (1:nveg) = hvt (1:nveg) hvb_table (1:nveg) = hvb (1:nveg) + z0mhvt_table (1:nveg) = z0mhvt (1:nveg) den_table (1:nveg) = den (1:nveg) rc_table (1:nveg) = rc (1:nveg) mfsno_table (1:nveg) = mfsno (1:nveg) diff --git a/physics/noahmpdrv.F90 b/physics/noahmpdrv.F90 index d232b759f..aad0d1ca5 100644 --- a/physics/noahmpdrv.F90 +++ b/physics/noahmpdrv.F90 @@ -450,7 +450,7 @@ subroutine noahmpdrv_run & integer :: iopt_pedo = 1 ! option for pedotransfer function integer :: iopt_crop = 0 ! option for crop model integer :: iopt_gla = 2 ! option for glacier treatment - integer :: iopt_z0m = 1 ! option for z0m treatment + integer :: iopt_z0m = 2 ! option for z0m treatment ! ! --- local inputs to noah-mp and glacier subroutines; listed in order in noah-mp call @@ -1335,6 +1335,7 @@ subroutine transfer_mp_parameters (vegtype,soiltype,slopetype, & parameters%z0mvt = z0mvt_table(vegtype) !momentum roughness length (m) parameters%hvt = hvt_table(vegtype) !top of canopy (m) parameters%hvb = hvb_table(vegtype) !bottom of canopy (m) + parameters%z0mhvt = z0mhvt_table(vegtype) !momentum roughness length (m) parameters%den = den_table(vegtype) !tree density (no. of trunks per m2) parameters%rc = rc_table(vegtype) !tree crown radius (m) parameters%mfsno = mfsno_table(vegtype) !snowmelt m parameter () diff --git a/physics/noahmptable.tbl b/physics/noahmptable.tbl index 02e59b37a..add1737be 100644 --- a/physics/noahmptable.tbl +++ b/physics/noahmptable.tbl @@ -59,6 +59,7 @@ z0mvt = 1.00, 0.15, 0.15, 0.15, 0.14, 0.50, 0.12, 0.06, 0.09, 0.50, 0.80, 0.85, 1.10, 1.09, 0.80, 0.00, 0.12, 0.50, 0.00, 0.10, 0.30, 0.20, 0.03, 0.00, 0.01, 0.00, 0.00, hvt = 15.0, 2.00, 2.00, 2.00, 1.50, 8.00, 1.00, 1.10, 1.10, 10.0, 16.0, 18.0, 20.0, 20.0, 16.0, 0.00, 0.50, 10.0, 0.00, 0.50, 4.00, 2.00, 0.50, 0.00, 0.10, 0.00, 0.00, hvb = 1.00, 0.10, 0.10, 0.10, 0.10, 0.15, 0.05, 0.10, 0.10, 0.10, 11.5, 7.00, 8.00, 8.50, 10.0, 0.00, 0.05, 0.10, 0.00, 0.10, 0.10, 0.10, 0.10, 0.00, 0.10, 0.00, 0.00, + z0mhvt= 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.00, 0.05, 0.05, 0.05, 0.05, 0.00, 0.05, 0.00, 0.00, den = 0.01, 25.0, 25.0, 25.0, 25.0, 25.0, 100., 10.0, 10.0, 0.02, 0.10, 0.28, 0.02, 0.28, 0.10, 0.01, 10.0, 0.10, 0.01, 1.00, 1.00, 1.00, 1.00, 0.00, 0.01, 0.01, 0.01, rc = 1.00, 0.08, 0.08, 0.08, 0.08, 0.08, 0.03, 0.12, 0.12, 3.00, 1.40, 1.20, 3.60, 1.20, 1.40, 0.01, 0.10, 1.40, 0.01, 0.30, 0.30, 0.30, 0.30, 0.00, 0.01, 0.01, 0.01, !mfsno = 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, @@ -218,6 +219,7 @@ z0mvt = 1.09, 1.10, 0.85, 0.80, 0.80, 0.20, 0.06, 0.60, 0.50, 0.12, 0.30, 0.15, 1.00, 0.14, 0.00, 0.00, 0.00, 0.30, 0.20, 0.03, hvt = 20.0, 20.0, 18.0, 16.0, 16.0, 1.10, 1.10, 13.0, 10.0, 1.00, 5.00, 2.00, 15.0, 1.50, 0.00, 0.00, 0.00, 4.00, 2.00, 0.50, hvb = 8.50, 8.00, 7.00, 11.5, 10.0, 0.10, 0.10, 0.10, 0.10, 0.05, 0.10, 0.10, 1.00, 0.10, 0.00, 0.00, 0.00, 0.30, 0.20, 0.10, + z0mhvt= 0.0545, 0.055, 0.047, 0.050, 0.050, 0.182, 0.0545, 0.046, 0.050, 0.120, 0.060, 0.075, 0.067, 0.093, 0.000, 0.000, 0.000, 0.075, 0.100, 0.060, den = 0.28, 0.02, 0.28, 0.10, 0.10, 10.0, 10.0, 10.0, 0.02, 100., 5.05, 25.0, 0.01, 25.0, 0.00, 0.01, 0.01, 1.00, 1.00, 1.00, rc = 1.20, 3.60, 1.20, 1.40, 1.40, 0.12, 0.12, 0.12, 3.00, 0.03, 0.75, 0.08, 1.00, 0.08, 0.00, 0.01, 0.01, 0.30, 0.30, 0.30, !mfsno = 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, 2.50, From 54b406ef06b3be00e5002187b507d1c19072d132 Mon Sep 17 00:00:00 2001 From: Michael Barlage Date: Mon, 1 May 2023 17:38:15 -0400 Subject: [PATCH 06/16] use blumel99 approach for bare soil for chen09 trs option --- physics/module_sf_noahmplsm.F90 | 19 +++++++------------ 1 file changed, 7 insertions(+), 12 deletions(-) diff --git a/physics/module_sf_noahmplsm.F90 b/physics/module_sf_noahmplsm.F90 index b6e751e31..51e7fe3f9 100644 --- a/physics/module_sf_noahmplsm.F90 +++ b/physics/module_sf_noahmplsm.F90 @@ -5846,10 +5846,10 @@ subroutine thermalz0(parameters, fveg, z0m, z0mg, zlvl, g_sigma = fveg**blumel_gamma + fveg*(1.0-fveg)*blumel_zeta ! Blumel99 eqn 22 cdmn = g_sigma*cdmn_v + (1.0-g_sigma)*cdmn_g ! Blumel99 eqn 21 - z0m_out = (zlvl - ezpd)*exp(-0.4/sqrt(cdmn)) ! Blumel99 eqn 24 + z0m_out = (zlvl - ezpd)*exp(-0.4/sqrt(cdmn)) ! Blumel99 eqn 24 kb_sigma_fveg = c_sigma_fveg/log((zlvl-ezpd)/z0m_out) - & log((zlvl-ezpd)/z0m_out) ! Blumel99 eqn 34 - z0h_out = z0m_out/exp(kb_sigma_fveg) + z0h_out = z0m_out/exp(kb_sigma_fveg) endif @@ -5861,11 +5861,6 @@ subroutine thermalz0(parameters, fveg, z0m, z0mg, zlvl, z0h_out = z0m_out - elseif (opt_trs == chen09) then - - czil = 10.0 ** (- 0.4 * parameters%hvt) - z0h_out = z0m_out * exp(-czil*0.4*258.2*sqrt(ustarx*z0m_out)) - elseif (opt_trs == tessel) then if (vegtyp <= 5) then @@ -5874,13 +5869,13 @@ subroutine thermalz0(parameters, fveg, z0m, z0mg, zlvl, z0h_out = z0m_out * 0.01 endif - elseif (opt_trs == blumel99) then + elseif (opt_trs == blumel99 .or. opt_trs == chen09) then reyn = ustarx*z0m_out/viscosity ! Blumel99 eqn 36c if (reyn > 2.0) then - kb_sigma_f0 = 2.46*reyn**0.25 - log(7.4) ! Blumel99 eqn 36a + kb_sigma_f0 = 2.46*reyn**0.25 - log(7.4) ! Blumel99 eqn 36a else - kb_sigma_f0 = - log(0.397) ! Blumel99 eqn 36b + kb_sigma_f0 = - log(0.397) ! Blumel99 eqn 36b endif z0h_out = max(z0m_out/exp(kb_sigma_f0),1.0e-6) @@ -5912,8 +5907,8 @@ subroutine thermalz0(parameters, fveg, z0m, z0mg, zlvl, elseif (opt_trs == blumel99) then - sigma_a = 1.0 - (0.5/(0.5+vaie)) * exp(-vaie**2/8.0) ! Blumel99 eqn 8 - kb_sigma_f1 = 16.4 * (sigma_a*vaie**3)**(-0.25) * & ! Blumel99 eqn 38 + sigma_a = 1.0 - (0.5/(0.5+vaie)) * exp(-vaie**2/8.0) ! Blumel99 eqn 8 + kb_sigma_f1 = 16.4 * (sigma_a*vaie**3)**(-0.25) * & ! Blumel99 eqn 38 sqrt(parameters%dleaf*ur/log((zlvl-zpd)/z0m_out)) z0h_out = z0m_out/exp(kb_sigma_f1) c_sigma_f1 = log((zlvl-zpd)/z0m_out)*(log((zlvl-zpd)/z0m_out)+kb_sigma_f1) ! Blumel99 eqn 39 From 591791cb1515b4c9087bacbbba9a1f0b6319a770 Mon Sep 17 00:00:00 2001 From: Michael Barlage Date: Mon, 1 May 2023 17:43:52 -0400 Subject: [PATCH 07/16] use log averaging for z0m --- physics/module_sf_noahmplsm.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/module_sf_noahmplsm.F90 b/physics/module_sf_noahmplsm.F90 index 51e7fe3f9..28e6460a2 100644 --- a/physics/module_sf_noahmplsm.F90 +++ b/physics/module_sf_noahmplsm.F90 @@ -5814,19 +5814,19 @@ subroutine thermalz0(parameters, fveg, z0m, z0mg, zlvl, if (opt_trs == z0heqz0m) then - z0m_out = fveg * z0m + (1.0 - fveg) * z0mg ! probably should be log + z0m_out = exp(fveg * log(z0m) + (1.0 - fveg) * log(z0mg)) z0h_out = z0m_out elseif (opt_trs == chen09) then - z0m_out = fveg * z0m + (1.0 - fveg) * z0mg ! probably should be log + z0m_out = exp(fveg * log(z0m) + (1.0 - fveg) * log(z0mg)) czil = 10.0 ** (- 0.4 * parameters%hvt) z0h_out = fveg * z0m * exp(-czil*0.4*258.2*sqrt(ustarx*z0m )) & + (1.0 - fveg) * z0mg * exp(-czil*0.4*258.2*sqrt(ustarx*z0mg)) elseif (opt_trs == tessel) then - z0m_out = fveg * z0m + (1.0 - fveg) * z0mg ! probably should be log + z0m_out = exp(fveg * log(z0m) + (1.0 - fveg) * log(z0mg)) if (vegtyp <= 5) then z0h_out = fveg * z0m + (1.0 - fveg) * z0mg * 0.1 else From e0460595b81b80c2c5a5c6a25df62cdb54500827 Mon Sep 17 00:00:00 2001 From: Michael Barlage Date: Mon, 1 May 2023 17:49:32 -0400 Subject: [PATCH 08/16] modify compositing for trs = chen09 and tessel --- physics/module_sf_noahmplsm.F90 | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/physics/module_sf_noahmplsm.F90 b/physics/module_sf_noahmplsm.F90 index 28e6460a2..090e120dd 100644 --- a/physics/module_sf_noahmplsm.F90 +++ b/physics/module_sf_noahmplsm.F90 @@ -5820,17 +5820,25 @@ subroutine thermalz0(parameters, fveg, z0m, z0mg, zlvl, elseif (opt_trs == chen09) then z0m_out = exp(fveg * log(z0m) + (1.0 - fveg) * log(z0mg)) - czil = 10.0 ** (- 0.4 * parameters%hvt) - z0h_out = fveg * z0m * exp(-czil*0.4*258.2*sqrt(ustarx*z0m )) & - + (1.0 - fveg) * z0mg * exp(-czil*0.4*258.2*sqrt(ustarx*z0mg)) + czil = 10.0 ** (- 0.4 * parameters%hvt) + + reyn = ustarx*z0m_out/viscosity ! Blumel99 eqn 36c + if (reyn > 2.0) then + kb_sigma_f0 = 2.46*reyn**0.25 - log(7.4) ! Blumel99 eqn 36a + else + kb_sigma_f0 = - log(0.397) ! Blumel99 eqn 36b + endif + + z0h_out = exp( fveg * log(z0m * exp(-czil*0.4*258.2*sqrt(ustarx*z0m))) + & + (1.0 - fveg) * log(max(z0m/exp(kb_sigma_f0),1.0e-6)) ) elseif (opt_trs == tessel) then z0m_out = exp(fveg * log(z0m) + (1.0 - fveg) * log(z0mg)) if (vegtyp <= 5) then - z0h_out = fveg * z0m + (1.0 - fveg) * z0mg * 0.1 + z0h_out = fveg * log(z0m) + (1.0 - fveg) * log(z0mg * 0.1) else - z0h_out = fveg * z0m * 0.01 + (1.0 - fveg) * z0mg * 0.1 + z0h_out = fveg * log(z0m * 0.01) + (1.0 - fveg) * log(z0mg * 0.1) endif elseif (opt_trs == blumel99) then From 8ccdf650e0035ac367cce0443e5c2e6f4a09eaa1 Mon Sep 17 00:00:00 2001 From: Michael Barlage Date: Mon, 1 May 2023 17:51:55 -0400 Subject: [PATCH 09/16] add limits to rb --- physics/module_sf_noahmplsm.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/physics/module_sf_noahmplsm.F90 b/physics/module_sf_noahmplsm.F90 index 090e120dd..48d5024cf 100644 --- a/physics/module_sf_noahmplsm.F90 +++ b/physics/module_sf_noahmplsm.F90 @@ -5028,8 +5028,7 @@ subroutine ragrb(parameters,iter ,vai ,rhoair ,hg ,tah , & !in tmprb = cwpc*50. / (1. - exp(-cwpc/2.)) rb = tmprb * sqrt(parameters%dleaf/uc) - rb = max(rb,20.0) -! rb = 200 + rb = min(max(rb, 5.0),50.0) ! limit rb to 5-50, typically rb<50 end subroutine ragrb From 60fb0e81846a54d5acb3e6fce6ba74f9ae34fb85 Mon Sep 17 00:00:00 2001 From: Michael Barlage Date: Mon, 1 May 2023 17:57:28 -0400 Subject: [PATCH 10/16] lower glacier snow limit to 100mm --- physics/module_sf_noahmp_glacier.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/physics/module_sf_noahmp_glacier.F90 b/physics/module_sf_noahmp_glacier.F90 index bd6b016f1..7b7512fa4 100644 --- a/physics/module_sf_noahmp_glacier.F90 +++ b/physics/module_sf_noahmp_glacier.F90 @@ -2646,9 +2646,9 @@ subroutine snowwater_glacier (nsnow ,nsoil ,imelt ,dt ,sfctmp , & !in !to obtain equilibrium state of snow in glacier region - if(sneqv > 2000.) then ! 2000 mm -> maximum water depth + if(sneqv > 100.) then ! 100 mm -> maximum water depth bdsnow = snice(0) / dzsnso(0) - snoflow = (sneqv - 2000.) + snoflow = (sneqv - 100.) snice(0) = snice(0) - snoflow dzsnso(0) = dzsnso(0) - snoflow/bdsnow snoflow = snoflow / dt From 63aed7ca255f8b420ce50523c73a79c51294dd21 Mon Sep 17 00:00:00 2001 From: Michael Barlage Date: Mon, 1 May 2023 17:59:44 -0400 Subject: [PATCH 11/16] bug fix for mixing ratio calculation in canres --- physics/module_sf_noahmplsm.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/module_sf_noahmplsm.F90 b/physics/module_sf_noahmplsm.F90 index 48d5024cf..38f7313e6 100644 --- a/physics/module_sf_noahmplsm.F90 +++ b/physics/module_sf_noahmplsm.F90 @@ -6179,7 +6179,7 @@ subroutine canres (parameters,par ,sfctmp,rcsoil ,eah ,sfcprs , & !in ! compute q2 and q2sat q2 = 0.622 * eah / (sfcprs - 0.378 * eah) !specific humidity [kg/kg] - q2 = q2 / (1.0 + q2) !mixing ratio [kg/kg] + q2 = q2 / (1.0 - q2) !mixing ratio [kg/kg] call calhum(parameters,sfctmp, sfcprs, q2sat, dqsdt2) From 12d8ab54d5aa33f22bab21cbf0186891acc2ab2b Mon Sep 17 00:00:00 2001 From: Michael Barlage Date: Mon, 1 May 2023 18:06:16 -0400 Subject: [PATCH 12/16] fix bug in infil for infiltration limit --- physics/module_sf_noahmplsm.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/module_sf_noahmplsm.F90 b/physics/module_sf_noahmplsm.F90 index 38f7313e6..22e177912 100644 --- a/physics/module_sf_noahmplsm.F90 +++ b/physics/module_sf_noahmplsm.F90 @@ -8635,7 +8635,7 @@ subroutine infil (parameters,nsoil ,dt ,zsoil ,sh2o ,sice , & !in call wdfcnd2 (parameters,wdf,wcnd,sh2o(1),sicemax,1) infmax = max (infmax,wcnd) - infmax = min (infmax,px) + infmax = min (infmax,px/dt) runsrf= max(0., qinsur - infmax) pddum = qinsur - runsrf From 9525a1650b804e09b95a8f1059ea881dc80a0b20 Mon Sep 17 00:00:00 2001 From: Michael Barlage Date: Mon, 1 May 2023 18:15:40 -0400 Subject: [PATCH 13/16] add canopy heat to mptable for tuning --- physics/module_sf_noahmplsm.F90 | 3 ++- physics/noahmp_tables.f90 | 9 ++++++--- physics/noahmpdrv.F90 | 1 + physics/noahmptable.tbl | 3 +++ 4 files changed, 12 insertions(+), 4 deletions(-) diff --git a/physics/module_sf_noahmplsm.F90 b/physics/module_sf_noahmplsm.F90 index 22e177912..336ae5a26 100644 --- a/physics/module_sf_noahmplsm.F90 +++ b/physics/module_sf_noahmplsm.F90 @@ -224,6 +224,7 @@ module module_sf_noahmplsm real (kind=kind_phys) :: rc !< tree crown radius (m) real (kind=kind_phys) :: mfsno !< snowmelt m parameter () real (kind=kind_phys) :: scffac !< snow cover factor (m) + real (kind=kind_phys) :: cbiom !< canopy biomass heat capacity parameter (m) real (kind=kind_phys) :: saim(12) !< monthly stem area index, one-sided real (kind=kind_phys) :: laim(12) !< monthly leaf area index, one-sided real (kind=kind_phys) :: sla !< single-side leaf area per kg [m2/kg] @@ -4245,7 +4246,7 @@ subroutine vege_flux(parameters,nsnow ,nsoil ,isnow ,vegtyp ,veg , & end if ! canopy heat capacity - hcv = 0.02*vaie*cwat + canliq*cwat/denh2o + canice*cice/denice !j/m2/k + hcv = parameters%cbiom*vaie*cwat + canliq*cwat/denh2o + canice*cice/denice !j/m2/k b = sav-irc-shc-evc-tr+pahv !additional w/m2 ! a = fveg*(4.*cir*tv**3 + csh + (cev+ctr)*destv) !volumetric heat capacity diff --git a/physics/noahmp_tables.f90 b/physics/noahmp_tables.f90 index cc7856af6..de207d0cc 100644 --- a/physics/noahmp_tables.f90 +++ b/physics/noahmp_tables.f90 @@ -49,6 +49,7 @@ module noahmp_tables real :: rc_table(mvt) !< tree crown radius (m) real :: mfsno_table(mvt) !< snowmelt curve parameter () real :: scffac_table(mvt) !< snow cover factor (m) + real :: cbiom_table(mvt) !< canopy biomass heat capacity parameter (m) real :: saim_table(mvt,12) !< monthly stem area index, one-sided real :: laim_table(mvt,12) !< monthly leaf area index, one-sided real :: sla_table(mvt) !< single-side leaf area per kg [m2/kg] @@ -325,7 +326,7 @@ subroutine read_mp_table_parameters lai_may, lai_jun, lai_jul, lai_aug, lai_sep, lai_oct, lai_nov, lai_dec, & rhol_vis, rhol_nir, rhos_vis, rhos_nir, taul_vis, taul_nir, taus_vis, taus_nir,& ch2op, dleaf, z0mvt, hvt, hvb, z0mhvt, & - den, rc, mfsno, scffac, xl, cwpvt, c3psn, kc25, & + den, rc, mfsno, scffac, cbiom, xl, cwpvt, c3psn, kc25, & akc, ko25, ako, avcmx, aqe, ltovrc, dilefc, dilefw, rmf25, sla, fragr, tmin, & vcmx25, tdlef, bp, mp, qe25, rms25, rmr25, arm, folnmx, wdpool, wrrat, mrp, & nroot, rgl, rs, hs, topt, rsmax, rtovrc, rswoodc, bf, wstrc, laimin, & @@ -334,7 +335,7 @@ subroutine read_mp_table_parameters namelist / noahmp_usgs_parameters / isurban, iswater, isbarren, isice, iscrop, eblforest, natural, & lcz_1, lcz_2, lcz_3, lcz_4, lcz_5, lcz_6, lcz_7, lcz_8, lcz_9, lcz_10, lcz_11, & ch2op, dleaf, z0mvt, hvt, hvb, z0mhvt, & - den, rc, mfsno, scffac, xl, cwpvt, c3psn, kc25, & + den, rc, mfsno, scffac, cbiom, xl, cwpvt, c3psn, kc25, & akc, ko25, ako, avcmx, aqe, ltovrc, dilefc, dilefw, rmf25, sla, fragr, tmin, & vcmx25, tdlef, bp, mp, qe25, rms25, rmr25, arm, folnmx, wdpool, wrrat, mrp, & nroot, rgl, rs, hs, topt, rsmax, rtovrc, rswoodc, bf, wstrc, laimin, & @@ -347,7 +348,7 @@ subroutine read_mp_table_parameters namelist / noahmp_modis_parameters / isurban, iswater, isbarren, isice, iscrop, eblforest, natural, & lcz_1, lcz_2, lcz_3, lcz_4, lcz_5, lcz_6, lcz_7, lcz_8, lcz_9, lcz_10, lcz_11, & ch2op, dleaf, z0mvt, hvt, hvb, z0mhvt, & - den, rc, mfsno, scffac, xl, cwpvt, c3psn, kc25, & + den, rc, mfsno, scffac, cbiom, xl, cwpvt, c3psn, kc25, & akc, ko25, ako, avcmx, aqe, ltovrc, dilefc, dilefw, rmf25, sla, fragr, tmin, & vcmx25, tdlef, bp, mp, qe25, rms25, rmr25, arm, folnmx, wdpool, wrrat, mrp, & nroot, rgl, rs, hs, topt, rsmax, rtovrc, rswoodc, bf, wstrc, laimin, & @@ -511,6 +512,7 @@ subroutine read_mp_table_parameters rc_table = -1.0e36 mfsno_table = -1.0e36 scffac_table = -1.0e36 + cbiom_table = -1.0e36 rhol_table = -1.0e36 rhos_table = -1.0e36 taul_table = -1.0e36 @@ -824,6 +826,7 @@ subroutine read_mp_table_parameters rc_table (1:nveg) = rc (1:nveg) mfsno_table (1:nveg) = mfsno (1:nveg) scffac_table (1:nveg) = scffac (1:nveg) + cbiom_table (1:nveg) = cbiom (1:nveg) xl_table (1:nveg) = xl (1:nveg) cwpvt_table (1:nveg) = cwpvt (1:nveg) c3psn_table (1:nveg) = c3psn (1:nveg) diff --git a/physics/noahmpdrv.F90 b/physics/noahmpdrv.F90 index aad0d1ca5..6831d17a2 100644 --- a/physics/noahmpdrv.F90 +++ b/physics/noahmpdrv.F90 @@ -1340,6 +1340,7 @@ subroutine transfer_mp_parameters (vegtype,soiltype,slopetype, & parameters%rc = rc_table(vegtype) !tree crown radius (m) parameters%mfsno = mfsno_table(vegtype) !snowmelt m parameter () parameters%scffac = scffac_table(vegtype) !snow cover factor + parameters%cbiom = cbiom_table(vegtype) !canopy biomass heat capacity parameter (m) parameters%saim = saim_table(vegtype,:) !monthly stem area index, one-sided parameters%laim = laim_table(vegtype,:) !monthly leaf area index, one-sided parameters%sla = sla_table(vegtype) !single-side leaf area per kg [m2/kg] diff --git a/physics/noahmptable.tbl b/physics/noahmptable.tbl index add1737be..3ffd5b532 100644 --- a/physics/noahmptable.tbl +++ b/physics/noahmptable.tbl @@ -67,6 +67,7 @@ mfsno = 4.00, 3.00, 3.00, 3.00, 4.00, 4.00, 2.00, 2.00, 2.00, 2.00, 1.00, 1.00, 1.00, 1.00, 1.00, 3.00, 3.00, 3.00, 3.00, 3.50, 3.50, 3.50, 3.50, 2.50, 3.50, 3.50, 3.50, ! c. he 12/17/2020: optimized snow cover factor (m) in scf formulation to replace original constant 2.5*z0,z0=0.002m, based on evaluation with snotel swe and modis scf, surface albedo scffac= 0.042, 0.014, 0.014, 0.014, 0.026, 0.026, 0.020, 0.018, 0.016, 0.020, 0.008, 0.008, 0.008, 0.008, 0.008, 0.030, 0.020, 0.020, 0.016, 0.030, 0.030, 0.030, 0.030, 0.030, 0.030, 0.030, 0.030, + cbiom = 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, ! row 1: vis ! row 2: near ir @@ -228,6 +229,8 @@ ! c. he 12/17/2020: optimized snow cover factor (m) in scf formulation to replace original constant 2.5*z0,z0=0.002m, based on evaluation with snotel swe and modis scf, surface albedo ! scffac = 0.008, 0.008, 0.008, 0.008, 0.008, 0.016, 0.016, 0.020, 0.020, 0.020, 0.020, 0.014, 0.042, 0.026, 0.030, 0.016, 0.030, 0.030, 0.030, 0.030, scffac = 0.005, 0.005, 0.005, 0.005, 0.005, 0.008, 0.008, 0.010, 0.010, 0.010, 0.010, 0.007, 0.021, 0.013, 0.015, 0.008, 0.015, 0.015, 0.015, 0.015, + cbiom = 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, + ! row 1: vis ! row 2: near ir rhol_vis=0.07, 0.10, 0.07, 0.10, 0.10, 0.07, 0.07, 0.07, 0.10, 0.11, 0.105, 0.11, 0.00, 0.11, 0.00, 0.00, 0.00, 0.10, 0.10, 0.10, From 6ab1b969002d265d0a82acd5f2fb261a390d3633 Mon Sep 17 00:00:00 2001 From: Michael Barlage Date: Mon, 1 May 2023 18:19:41 -0400 Subject: [PATCH 14/16] add two changes missed in previous update --- physics/module_sf_noahmplsm.F90 | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/physics/module_sf_noahmplsm.F90 b/physics/module_sf_noahmplsm.F90 index 336ae5a26..67e24091e 100644 --- a/physics/module_sf_noahmplsm.F90 +++ b/physics/module_sf_noahmplsm.F90 @@ -9103,7 +9103,8 @@ subroutine groundwater(parameters,nsnow ,nsoil ,dt ,sice ,zsoil , & !in fff = parameters%bexp(iwt) / 3.0 ! calibratable, c.he changed based on gy niu's update rsbmx = hk(iwt) * 1.0e3 * exp(3.0) ! mm/s, calibratable, c.he changed based on gy niu's update - qdis = (1.0-fcrmax)*rsbmx*exp(-parameters%timean)*exp(-fff*(zwt-2.0)) +! qdis = (1.0-fcrmax)*rsbmx*exp(-parameters%timean)*exp(-fff*(zwt-2.0)) + qdis = (1.0-fcrmax)*rsbmx*exp(-parameters%timean)*exp(-fff*zwt) ! c.he changed based on gy niu's update ! matric potential at the layer above the water table @@ -9114,7 +9115,9 @@ subroutine groundwater(parameters,nsnow ,nsoil ,dt ,sice ,zsoil , & !in ! recharge rate qin to groundwater - ka = hk(iwt) +! ka = hk(iwt) +! harmonic average, c.he changed based on gy niu's update + ka = 2.0*(hk(iwt)*parameters%dksat(iwt)*1.0e3) / (hk(iwt)+parameters%dksat(iwt)*1.0e3) wh_zwt = - zwt * 1.e3 !(mm) wh = smpfz - znode(iwt)*1.e3 !(mm) From 36e389231736c27a322c52cb0f9b81c606210e74 Mon Sep 17 00:00:00 2001 From: Michael Barlage Date: Mon, 1 May 2023 18:25:07 -0400 Subject: [PATCH 15/16] fix inout on ustar in sfcdif 1 and 2 --- physics/module_sf_noahmp_glacier.F90 | 2 +- physics/module_sf_noahmplsm.F90 | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/module_sf_noahmp_glacier.F90 b/physics/module_sf_noahmp_glacier.F90 index 7b7512fa4..492c4a50d 100644 --- a/physics/module_sf_noahmp_glacier.F90 +++ b/physics/module_sf_noahmp_glacier.F90 @@ -1583,7 +1583,7 @@ subroutine sfcdif1_glacier(iter ,zlvl ,zpd ,z0h ,z0m , & !in #endif ! outputs - real (kind=kind_phys), intent(out) :: fv !< friction velocity (m/s) + real (kind=kind_phys), intent(inout) :: fv !< friction velocity (m/s) real (kind=kind_phys), intent(out) :: cm !< drag coefficient for momentum real (kind=kind_phys), intent(out) :: ch !< drag coefficient for heat real (kind=kind_phys), intent(out) :: ch2 !< drag coefficient for heat diff --git a/physics/module_sf_noahmplsm.F90 b/physics/module_sf_noahmplsm.F90 index 67e24091e..6c199531d 100644 --- a/physics/module_sf_noahmplsm.F90 +++ b/physics/module_sf_noahmplsm.F90 @@ -5084,7 +5084,7 @@ subroutine sfcdif1(parameters,iter ,sfctmp ,rhoair ,h ,qair , & !in real (kind=kind_phys), intent(out) :: cm !< drag coefficient for momentum real (kind=kind_phys), intent(out) :: ch !< drag coefficient for heat - real (kind=kind_phys), intent(out) :: fv !< friction velocity (m/s) + real (kind=kind_phys), intent(inout) :: fv !< friction velocity (m/s) real (kind=kind_phys), intent(out) :: ch2 !< drag coefficient for heat ! locals @@ -5239,7 +5239,7 @@ subroutine sfcdif2(parameters,iter ,z0 ,thz0 ,thlm ,sfcspd , & !in real (kind=kind_phys), intent(inout) :: akhs real (kind=kind_phys), intent(inout) :: rlmo real (kind=kind_phys), intent(inout) :: wstar2 - real (kind=kind_phys), intent(out) :: ustar + real (kind=kind_phys), intent(inout) :: ustar real (kind=kind_phys) zz, pslmu, pslms, pslhu, pslhs real (kind=kind_phys) xx, pspmu, yy, pspms, psphu, psphs From 2fc95fef3e35dcd114c323c214fedc014845d54d Mon Sep 17 00:00:00 2001 From: Michael Barlage Date: Mon, 1 May 2023 18:29:46 -0400 Subject: [PATCH 16/16] remove zvfun and gdx effect from sfcdif3 --- physics/module_sf_noahmplsm.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/physics/module_sf_noahmplsm.F90 b/physics/module_sf_noahmplsm.F90 index 6c199531d..1888a26e8 100644 --- a/physics/module_sf_noahmplsm.F90 +++ b/physics/module_sf_noahmplsm.F90 @@ -5501,6 +5501,9 @@ subroutine sfcdif3(parameters,iloc ,jloc ,iter ,sfctmp ,qair ,ur tem2 = max(fveg, 0.1_kind_phys) zvfun1 = sqrt(tem1 * tem2) gdx = sqrt(garea1) + + gdx = 3000.0 ! this will remove gdx effect + zvfun1 = 1.0 ! this will remove zvfun effect if(thsfc_loc) then ! Use local potential temperature tvs = tgb * virtfac