From 902f1cb76f6371da5f3215ed6e9cd358b3f1f54e Mon Sep 17 00:00:00 2001 From: Keith Oleson Date: Mon, 21 Mar 2022 14:07:31 -0600 Subject: [PATCH] Update forcing heights and restrict zeta to less than 20 for now --- src/biogeophys/BareGroundFluxesMod.F90 | 11 ++++++- src/biogeophys/CanopyFluxesMod.F90 | 7 +++- src/biogeophys/FrictionVelocityMod.F90 | 44 +++++++++++++++++++++----- src/biogeophys/LakeFluxesMod.F90 | 22 +++++++++++-- 4 files changed, 71 insertions(+), 13 deletions(-) diff --git a/src/biogeophys/BareGroundFluxesMod.F90 b/src/biogeophys/BareGroundFluxesMod.F90 index 11d1790fc9..73c169fa13 100644 --- a/src/biogeophys/BareGroundFluxesMod.F90 +++ b/src/biogeophys/BareGroundFluxesMod.F90 @@ -234,7 +234,9 @@ subroutine BareGroundFluxes(bounds, num_noexposedvegp, filter_noexposedvegp, & rh_ref2m_r => waterdiagnosticbulk_inst%rh_ref2m_r_patch , & ! Output: [real(r8) (:) ] Rural 2 m height surface relative humidity (%) rh_ref2m => waterdiagnosticbulk_inst%rh_ref2m_patch , & ! Output: [real(r8) (:) ] 2 m height surface relative humidity (%) - forc_hgt_u_patch => frictionvel_inst%forc_hgt_u_patch , & ! Input: + forc_hgt_u_patch => frictionvel_inst%forc_hgt_u_patch , & ! Output: [real(r8) (:) ] observational height of wind at patch level [m] + forc_hgt_t_patch => frictionvel_inst%forc_hgt_t_patch , & ! Output: [real(r8) (:) ] observational height of temperature at patch level [m] + forc_hgt_q_patch => frictionvel_inst%forc_hgt_q_patch , & ! Output: [real(r8) (:) ] observational height of specific humidity at patch level [m] u10_clm => frictionvel_inst%u10_clm_patch , & ! Input: [real(r8) (:) ] 10 m height winds (m/s) zetamax => frictionvel_inst%zetamaxstable , & ! Input: [real(r8) ] max zeta value under stable conditions z0mg_col => frictionvel_inst%z0mg_col , & ! Output: [real(r8) (:) ] roughness length, momentum [m] @@ -359,6 +361,13 @@ subroutine BareGroundFluxes(bounds, num_noexposedvegp, filter_noexposedvegp, & end select z0qg_patch(p) = z0hg_patch(p) + ! Update the forcing heights for new roughness lengths + ! TODO(KWO, 2022-03-15) Only for Meier2022 for now to maintain bfb with ZengWang2007 + if (z0param_method == 'Meier2022') then + forc_hgt_u_patch(p) = forc_hgt_u_patch(g) + z0mg_patch(p) + displa(p) + forc_hgt_t_patch(p) = forc_hgt_t_patch(g) + z0hg_patch(p) + displa(p) + forc_hgt_q_patch(p) = forc_hgt_q_patch(g) + z0qg_patch(p) + displa(p) + end if thvstar = tstar*(1._r8+0.61_r8*forc_q(c)) + 0.61_r8*forc_th(c)*qstar zeta = zldis(p)*vkc*grav*thvstar/(ustar(p)**2*thv(c)) diff --git a/src/biogeophys/CanopyFluxesMod.F90 b/src/biogeophys/CanopyFluxesMod.F90 index 0a945f975c..31d73e1b37 100644 --- a/src/biogeophys/CanopyFluxesMod.F90 +++ b/src/biogeophys/CanopyFluxesMod.F90 @@ -1393,7 +1393,12 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, if (zeta(p) >= 0._r8) then !stable ! remove stability cap when biomass heat storage is active if(use_biomass_heat_storage) then - zeta(p) = min(100._r8,max(zeta(p),0.01_r8)) + ! TODO(KWO, 2022-03-15) Only for Meier2022 for now to maintain bfb with ZengWang2007 + if (z0param_method == 'Meier2022') then + zeta(p) = min(20._r8,max(zeta(p),0.01_r8)) + else + zeta(p) = min(100._r8,max(zeta(p),0.01_r8)) + end if else zeta(p) = min(zetamax,max(zeta(p),0.01_r8)) endif diff --git a/src/biogeophys/FrictionVelocityMod.F90 b/src/biogeophys/FrictionVelocityMod.F90 index 585659a981..5d532b1145 100644 --- a/src/biogeophys/FrictionVelocityMod.F90 +++ b/src/biogeophys/FrictionVelocityMod.F90 @@ -703,17 +703,17 @@ subroutine SetRoughnessLengthsAndForcHeightsNonLake(this, bounds, & if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then if (frac_veg_nosno(p) == 0) then forc_hgt_u_patch(p) = forc_hgt_u(g) + z0mg(c) + displa(p) - forc_hgt_t_patch(p) = forc_hgt_t(g) + z0mg(c) + displa(p) - forc_hgt_q_patch(p) = forc_hgt_q(g) + z0mg(c) + displa(p) + forc_hgt_t_patch(p) = forc_hgt_t(g) + z0hg(c) + displa(p) + forc_hgt_q_patch(p) = forc_hgt_q(g) + z0qg(c) + displa(p) else - forc_hgt_u_patch(p) = forc_hgt_u(g) + z0m(p) + displa(p) - forc_hgt_t_patch(p) = forc_hgt_t(g) + z0m(p) + displa(p) - forc_hgt_q_patch(p) = forc_hgt_q(g) + z0m(p) + displa(p) + forc_hgt_u_patch(p) = forc_hgt_u(g) + z0mv(p) + displa(p) + forc_hgt_t_patch(p) = forc_hgt_t(g) + z0hv(p) + displa(p) + forc_hgt_q_patch(p) = forc_hgt_q(g) + z0qv(p) + displa(p) end if else if (lun%itype(l) == istwet .or. lun%itype(l) == istice) then - forc_hgt_u_patch(p) = forc_hgt_u(g) + z0mg(c) - forc_hgt_t_patch(p) = forc_hgt_t(g) + z0mg(c) - forc_hgt_q_patch(p) = forc_hgt_q(g) + z0mg(c) + forc_hgt_u_patch(p) = forc_hgt_u(g) + z0mg(c) + displa(p) + forc_hgt_t_patch(p) = forc_hgt_t(g) + z0hg(c) + displa(p) + forc_hgt_q_patch(p) = forc_hgt_q(g) + z0qg(c) + displa(p) else if (urbpoi(l)) then forc_hgt_u_patch(p) = forc_hgt_u(g) + z_0_town(l) + z_d_town(l) forc_hgt_t_patch(p) = forc_hgt_t(g) + z_0_town(l) + z_d_town(l) @@ -893,6 +893,10 @@ subroutine FrictionVelocity(this, lbn, ubn, fn, filtern, & zldis(n) = forc_hgt_u_patch(n)-displa(n) end if zeta(n) = zldis(n)/obu(n) + ! TODO(KWO, 2022-03-15) Only for Meier2022 for now to maintain bfb with ZengWang2007 + if (z0param_method == 'Meier2022') then + zeta(n) = min(20._r8,zldis(n)/obu(n)) + end if if (zeta(n) < -zetam) then ustar(n) = vkc*um(n)/(log(-zetam*obu(n)/z0m(n))& - this%StabilityFunc1(-zetam) & @@ -992,6 +996,10 @@ subroutine FrictionVelocity(this, lbn, ubn, fn, filtern, & zldis(n) = forc_hgt_t_patch(n)-displa(n) end if zeta(n) = zldis(n)/obu(n) + ! TODO(KWO, 2022-03-15) Only for Meier2022 for now to maintain bfb with ZengWang2007 + if (z0param_method == 'Meier2022') then + zeta(n) = min(20._r8,zldis(n)/obu(n)) + end if if (zeta(n) < -zetat) then temp1(n) = vkc/(log(-zetat*obu(n)/z0h(n))& - this%StabilityFunc2(-zetat) & @@ -1016,6 +1024,10 @@ subroutine FrictionVelocity(this, lbn, ubn, fn, filtern, & else zldis(n) = forc_hgt_q_patch(pfti(n))-displa(n) zeta(n) = zldis(n)/obu(n) + ! TODO(KWO, 2022-03-15) Only for Meier2022 for now to maintain bfb with ZengWang2007 + if (z0param_method == 'Meier2022') then + zeta(n) = min(20._r8,zldis(n)/obu(n)) + end if if (zeta(n) < -zetat) then temp2(n) = vkc/(log(-zetat*obu(n)/z0q(n)) & - this%StabilityFunc2(-zetat) & @@ -1038,6 +1050,10 @@ subroutine FrictionVelocity(this, lbn, ubn, fn, filtern, & else zldis(n) = forc_hgt_q_patch(n)-displa(n) zeta(n) = zldis(n)/obu(n) + ! TODO(KWO, 2022-03-15) Only for Meier2022 for now to maintain bfb with ZengWang2007 + if (z0param_method == 'Meier2022') then + zeta(n) = min(20._r8,zldis(n)/obu(n)) + end if if (zeta(n) < -zetat) then temp2(n) = vkc/(log(-zetat*obu(n)/z0q(n)) & - this%StabilityFunc2(-zetat) & @@ -1060,6 +1076,10 @@ subroutine FrictionVelocity(this, lbn, ubn, fn, filtern, & zldis(n) = 2.0_r8 + z0h(n) zeta(n) = zldis(n)/obu(n) + ! TODO(KWO, 2022-03-15) Only for Meier2022 for now to maintain bfb with ZengWang2007 + if (z0param_method == 'Meier2022') then + zeta(n) = min(20._r8,zldis(n)/obu(n)) + end if if (zeta(n) < -zetat) then temp12m(n) = vkc/(log(-zetat*obu(n)/z0h(n))& - this%StabilityFunc2(-zetat) & @@ -1083,6 +1103,10 @@ subroutine FrictionVelocity(this, lbn, ubn, fn, filtern, & else zldis(n) = 2.0_r8 + z0q(n) zeta(n) = zldis(n)/obu(n) + ! TODO(KWO, 2022-03-15) Only for Meier2022 for now to maintain bfb with ZengWang2007 + if (z0param_method == 'Meier2022') then + zeta(n) = min(20._r8,zldis(n)/obu(n)) + end if if (zeta(n) < -zetat) then temp22m(n) = vkc/(log(-zetat*obu(n)/z0q(n)) - & this%StabilityFunc2(-zetat) + this%StabilityFunc2(z0q(n)/obu(n)) & @@ -1111,6 +1135,10 @@ subroutine FrictionVelocity(this, lbn, ubn, fn, filtern, & zldis(n) = forc_hgt_u_patch(n)-displa(n) end if zeta(n) = zldis(n)/obu(n) + ! TODO(KWO, 2022-03-15) Only for Meier2022 for now to maintain bfb with ZengWang2007 + if (z0param_method == 'Meier2022') then + zeta(n) = min(20._r8,zldis(n)/obu(n)) + end if if (min(zeta(n), 1._r8) < 0._r8) then tmp1 = (1._r8 - 16._r8*min(zeta(n),1._r8))**0.25_r8 tmp2 = log((1._r8+tmp1*tmp1)/2._r8) diff --git a/src/biogeophys/LakeFluxesMod.F90 b/src/biogeophys/LakeFluxesMod.F90 index b53b0fd23b..6dc9d1f0c0 100644 --- a/src/biogeophys/LakeFluxesMod.F90 +++ b/src/biogeophys/LakeFluxesMod.F90 @@ -380,9 +380,17 @@ subroutine LakeFluxes(bounds, num_lakec, filter_lakec, num_lakep, filter_lakep, ! Surface temperature and fluxes - forc_hgt_u_patch(p) = forc_hgt_u(g) + z0mg(p) - forc_hgt_t_patch(p) = forc_hgt_t(g) + z0mg(p) - forc_hgt_q_patch(p) = forc_hgt_q(g) + z0mg(p) + ! Update forcing heights for updated roughness lengths + ! TODO(KWO, 2022-03-15) Only for Meier2022 for now to maintain bfb with ZengWang2007 + if (z0param_method == 'Meier2022') then + forc_hgt_u_patch(p) = forc_hgt_u(g) + z0mg(p) + forc_hgt_t_patch(p) = forc_hgt_t(g) + z0hg(p) + forc_hgt_q_patch(p) = forc_hgt_q(g) + z0qg(p) + else + forc_hgt_u_patch(p) = forc_hgt_u(g) + z0mg(p) + forc_hgt_t_patch(p) = forc_hgt_t(g) + z0mg(p) + forc_hgt_q_patch(p) = forc_hgt_q(g) + z0mg(p) + end if ! Find top layer jtop(c) = snl(c) + 1 @@ -622,6 +630,14 @@ subroutine LakeFluxes(bounds, num_lakec, filter_lakec, num_lakep, filter_lakep, z0qg(p) = z0hg(p) end if + ! Update forcing heights for updated roughness lengths + ! TODO(KWO, 2022-03-15) Only for Meier2022 for now to maintain bfb with ZengWang2007 + if (z0param_method == 'Meier2022') then + forc_hgt_u_patch(p) = forc_hgt_u(g) + z0mg(p) + forc_hgt_t_patch(p) = forc_hgt_t(g) + z0hg(p) + forc_hgt_q_patch(p) = forc_hgt_q(g) + z0qg(p) + end if + end do ! end of filtered pft loop iter = iter + 1