Skip to content

Commit

Permalink
Merge branch 'main' of github.com:RonnyMeier/CTSM into main
Browse files Browse the repository at this point in the history 
 Conflicts:
	cime_config/testdefs/testmods_dirs/clm/Meier2022_surf_rough/user_nl_clm
	src/biogeophys/FrictionVelocityMod.F90
  • Loading branch information
@ekluzek
ekluzek committed Oct 19, 2022
2 parents 357fea5 + 6be9820 commit a924de8
Show file tree
Hide file tree
Showing 12 changed files with 130 additions and 88 deletions.
5 changes: 4 additions & 1 deletion bld/CLMBuildNamelist.pm
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3885,7 +3885,10 @@ sub setup_logic_friction_vel {
#
my ($opts, $nl_flags, $definition, $defaults, $nl) = @_;

add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'zetamaxstable' );
my $z0param_method = remove_leading_and_trailing_quotes($nl->get_value('z0param_method' ));

add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'zetamaxstable',
'z0param_method'=>$z0param_method );
}

#-------------------------------------------------------------------------------
Expand Down
3 changes: 2 additions & 1 deletion bld/namelist_files/namelist_defaults_ctsm.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -208,7 +208,8 @@ attributes from the config_cache.xml file (with keys converted to upper-case).
<!-- Friction velocity -->
<zetamaxstable phys="clm4_5" >2.0d00</zetamaxstable>
<zetamaxstable phys="clm5_0" >0.5d00</zetamaxstable>
<zetamaxstable phys="clm5_1" >0.5d00</zetamaxstable>
<zetamaxstable phys="clm5_1" z0param_method="ZengWang2007" >0.5d00</zetamaxstable>
<zetamaxstable phys="clm5_1" z0param_method="Meier2022" >2.0d00</zetamaxstable>

<!-- atm2lnd defaults -->
<repartition_rain_snow phys="clm5_1" >.true.</repartition_rain_snow>
Expand Down
2 changes: 1 addition & 1 deletion cime_config/testdefs/testmods_dirs/clm/Meier2022_surf_rough/user_nl_clm
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
z0param_method = 'Meier2022'
use_z0m_snowmelt = .true.
calc_human_stress_indices = 'NONE' ! Currently dies when turned on because of a negative humidity (about -31) in Wet Bulb calculation
use_z0mg_2d = .false.
paramfile = '$DIN_LOC_ROOT/lnd/clm2/paramdata/ctsm51_params.RMz0.c220304.nc'

23 changes: 16 additions & 7 deletions src/biogeophys/BareGroundFluxesMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,7 +117,6 @@ subroutine BareGroundFluxes(bounds, num_noexposedvegp, filter_noexposedvegp, &
integer :: p,c,g,f,j,l ! indices
integer :: iter ! iteration index
real(r8) :: zldis(bounds%begp:bounds%endp) ! reference height "minus" zero displacement height [m]
real(r8) :: zeta ! dimensionless height used in Monin-Obukhov theory
real(r8) :: wc ! convective velocity [m/s]
real(r8) :: dth(bounds%begp:bounds%endp) ! diff of virtual temp. between ref. height and surface
real(r8) :: dthv ! diff of vir. poten. temp. between ref. height and surface
Expand Down Expand Up @@ -234,9 +233,12 @@ 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
zeta => frictionvel_inst%zeta_patch , & ! Output: [real(r8) (:) ] dimensionless stability parameter
z0mg_col => frictionvel_inst%z0mg_col , & ! Output: [real(r8) (:) ] roughness length, momentum [m]
z0hg_col => frictionvel_inst%z0hg_col , & ! Output: [real(r8) (:) ] roughness length, sensible heat [m]
z0qg_col => frictionvel_inst%z0qg_col , & ! Output: [real(r8) (:) ] roughness length, latent heat [m]
Expand Down Expand Up @@ -359,18 +361,25 @@ 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))
zeta(p) = zldis(p)*vkc*grav*thvstar/(ustar(p)**2*thv(c))

if (zeta >= 0._r8) then !stable
zeta = min(zetamax,max(zeta,0.01_r8))
if (zeta(p) >= 0._r8) then !stable
zeta(p) = min(zetamax,max(zeta(p),0.01_r8))
um(p) = max(ur(p),0.1_r8)
else !unstable
zeta = max(-100._r8,min(zeta,-0.01_r8))
zeta(p) = max(-100._r8,min(zeta(p),-0.01_r8))
wc = beta(c)*(-grav*ustar(p)*thvstar*zii(c)/thv(c))**0.333_r8
um(p) = sqrt(ur(p)*ur(p) + wc*wc)
end if
obu(p) = zldis(p)/zeta
obu(p) = zldis(p)/zeta(p)

num_iter(p) = iter
end do
Expand Down
22 changes: 15 additions & 7 deletions src/biogeophys/CanopyFluxesMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -918,14 +918,17 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp,
call endrun(msg = 'unknown z0param_method', additional_msg = errMsg(sourcefile, __LINE__))
end select

! Update the forcing heights
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)

z0hv(p) = z0mv(p)
z0qv(p) = z0mv(p)

! Update the forcing heights
! 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) + 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

end do

found = .false.
Expand Down Expand Up @@ -1390,15 +1393,20 @@ 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(zetamax,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
um(p) = max(ur(p),0.1_r8)
else !unstable
zeta(p) = max(-100._r8,min(zeta(p),-0.01_r8))
if ( ustar(p)*thvstar > 0.0d00 )then
write(iulog,*) 'ustar*thvstart is positive and has to be negative'
write(iulog,*) 'ustar*thvstar is positive and has to be negative'
write(iulog,*) 'p = ', p
write(iulog,*) '-grav*ustar(p)*thvstar*zii/thv(c) = ', -grav*ustar(p)*thvstar*zii/thv(c)
write(iulog,*) 'ustar = ', ustar(p)
Expand Down
28 changes: 14 additions & 14 deletions src/biogeophys/FrictionVelocityMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -555,7 +555,7 @@ subroutine SetRoughnessLengthsAndForcHeightsNonLake(this, bounds, &

frac_veg_nosno => canopystate_inst%frac_veg_nosno_patch , & ! Input: [integer (:) ] fraction of vegetation not covered by snow (0 OR 1) [-]
frac_sno => waterdiagnosticbulk_inst%frac_sno_col , & ! Input: [real(r8) (:) ] fraction of ground covered by snow (0 to 1)
snomelt_accum => waterfluxbulk_inst%qflx_snomelt_accum_col , & ! Input: [real(r8) (:) ] accumulated col snow melt for z0m calculation (m H2O)
snomelt_accum => waterdiagnosticbulk_inst%snomelt_accum_col , & ! Input: [real(r8) (:) ] accumulated col snow melt for z0m calculation (m H2O)
urbpoi => lun%urbpoi , & ! Input: [logical (:) ] true => landunit is an urban point
z_0_town => lun%z_0_town , & ! Input: [real(r8) (:) ] momentum roughness length of urban landunit (m)
z_d_town => lun%z_d_town , & ! Input: [real(r8) (:) ] displacement height of urban landunit (m)
Expand All @@ -576,7 +576,11 @@ subroutine SetRoughnessLengthsAndForcHeightsNonLake(this, bounds, &
case ('ZengWang2007')
if (frac_sno(c) > 0._r8) then
if(use_z0m_snowmelt) then
z0mg(c) = exp(1.4_r8 * (atan((log10(snomelt_accum(c))+0.23_r8)/0.08_r8))-0.31_r8) / 1000._r8
if ( snomelt_accum(c) < 1.e-5_r8 )then
z0mg(c) = exp(1.4_r8 * -rpi/2.0_r8 -0.31_r8) / 1000._r8
else
z0mg(c) = exp(1.4_r8 * (atan((log10(snomelt_accum(c))+0.23_r8)/0.08_r8))-0.31_r8) / 1000._r8
end if
else
z0mg(c) = this%zsno
end if
Expand All @@ -595,13 +599,9 @@ subroutine SetRoughnessLengthsAndForcHeightsNonLake(this, bounds, &
end if
else
z0mg(c) = this%zsno


end if
else if (lun%itype(l) == istice) then
z0mg(c) = this%zglc


else
z0mg(c) = this%zlnd
end if
Expand Down Expand Up @@ -639,17 +639,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)
Expand Down
43 changes: 31 additions & 12 deletions src/biogeophys/LakeFluxesMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -170,7 +170,6 @@ subroutine LakeFluxes(bounds, num_lakec, filter_lakec, num_lakep, filter_lakep,
real(r8) :: ur(bounds%begp:bounds%endp) ! wind speed at reference height [m/s]
real(r8) :: ustar(bounds%begp:bounds%endp) ! friction velocity [m/s]
real(r8) :: wc ! convective velocity [m/s]
real(r8) :: zeta ! dimensionless height used in Monin-Obukhov theory
real(r8) :: zldis(bounds%begp:bounds%endp) ! reference height "minus" zero displacement height [m]
real(r8) :: displa(bounds%begp:bounds%endp) ! displacement (always zero) [m]
real(r8) :: u2m ! 2 m wind speed (m/s)
Expand Down Expand Up @@ -228,8 +227,7 @@ subroutine LakeFluxes(bounds, num_lakec, filter_lakec, num_lakep, filter_lakep,

h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] liquid water (kg/m2)
h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:) ] ice lens (kg/m2)
snomelt_accum => waterfluxbulk_inst%qflx_snomelt_accum_col , & ! Input: [real(r8) (:) ] accumulated col snow melt for z0m calculation (m H2O)

snomelt_accum => waterdiagnosticbulk_inst%snomelt_accum_col , & ! Input: [real(r8) (:) ] accumulated col snow melt for z0m calculation (m H2O)
t_skin_patch => temperature_inst%t_skin_patch , & ! Output: [real(r8) (:) ] patch skin temperature (K)

t_lake => temperature_inst%t_lake_col , & ! Input: [real(r8) (:,:) ] lake temperature (Kelvin)
Expand All @@ -240,6 +238,7 @@ subroutine LakeFluxes(bounds, num_lakec, filter_lakec, num_lakep, filter_lakep,
forc_hgt_t_patch => frictionvel_inst%forc_hgt_t_patch , & ! Input: [real(r8) (:) ] observational height of temperature at pft level [m]
forc_hgt_q_patch => frictionvel_inst%forc_hgt_q_patch , & ! Input: [real(r8) (:) ] observational height of specific humidity at pft level [m]
zetamax => frictionvel_inst%zetamaxstable , & ! Input: [real(r8) ] max zeta value under stable conditions
zeta => frictionvel_inst%zeta_patch , & ! Output: [real(r8) (:) ] dimensionless stability parameter
ram1 => frictionvel_inst%ram1_patch , & ! Output: [real(r8) (:) ] aerodynamical resistance (s/m)

q_ref2m => waterdiagnosticbulk_inst%q_ref2m_patch , & ! Output: [real(r8) (:) ] 2 m height surface specific humidity (kg/kg)
Expand Down Expand Up @@ -366,7 +365,11 @@ subroutine LakeFluxes(bounds, num_lakec, filter_lakec, num_lakep, filter_lakep,

case ('ZengWang2007')
if(use_z0m_snowmelt) then
z0mg(p) = exp(1.4_r8 * (atan((log10(snomelt_accum(c))+0.23_r8)/0.08_r8))-0.31_r8) / 1000._r8
if ( snomelt_accum(c) < 1.e-5_r8 ) then
z0mg(p) = exp(1.4_r8 * -rpi/2.0_r8 -0.31_r8) / 1000._r8
else
z0mg(p) = exp(1.4_r8 * (atan((log10(snomelt_accum(c))+0.23_r8)/0.08_r8))-0.31_r8) / 1000._r8
end if
else
z0mg(p) = params_inst%zsno
end if
Expand All @@ -377,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
Expand Down Expand Up @@ -543,17 +554,17 @@ subroutine LakeFluxes(bounds, num_lakec, filter_lakec, num_lakep, filter_lakep,
qstar = temp2(p)*dqh(p)

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))
zeta(p)=zldis(p)*vkc * grav*thvstar/(ustar(p)**2*thv(c))

if (zeta >= 0._r8) then !stable
zeta = min(zetamax,max(zeta,0.01_r8))
if (zeta(p) >= 0._r8) then !stable
zeta(p) = min(zetamax,max(zeta(p),0.01_r8))
um(p) = max(ur(p),0.1_r8)
else !unstable
zeta = max(-100._r8,min(zeta,-0.01_r8))
zeta(p) = max(-100._r8,min(zeta(p),-0.01_r8))
wc = beta1*(-grav*ustar(p)*thvstar*zii/thv(c))**0.333_r8
um(p) = sqrt(ur(p)*ur(p)+wc*wc)
end if
obu(p) = zldis(p)/zeta
obu(p) = zldis(p)/zeta(p)

if (obuold(p)*obu(p) < 0._r8) nmozsgn(p) = nmozsgn(p)+1

Expand Down Expand Up @@ -619,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
Expand Down
Loading

0 comments on commit a924de8

Please sign in to comment.