CTSM
CTSM copied to clipboard
ESCOMP
Surface roughness modifications
This branch contains the surface roughness (z0) modifications published in: https://doi.org/10.5194/gmd-15-2365-2022
When changing the namelist input z0param_method from ZengWang2007 (default) to Meier2022 the following modifications are activated:
- A new parameterization of the vegetation surface roughness based on Raupach (1992) with optimized parameters to match the data collected in Hu et al. (2020) for different types of vegetation. This requires several new PFT-specific input parameters in the parameter file.
- A spatially explicit z0m input field for bare soil based on the data of Prigent et al. (2005). This may be activated specifically by the user through the namelist input use_z0mg_2d. This requires a new input variable in the surfdata file.
- The parameterization of z0m for snow based on accumulated snow melt as proposed in Brock et al. (2006). This may be activated specifically by the user through the namelist input use_z0m_snowmelt.
- The parameterization of Yang et al. (2008) for z0h and z0q over bare soil, snow, and glaciers.
Note that the study in GMD also proposes new globally constant values for the z0m of bare soil, snow, and ice. To "activate" those the parameter file needs to be changed at the moment. The original and modified parameter files and surfdata files will be shared by ftp.
Open issues/questions (already discussed with @ekluzek, @dlawrenncar, and @olyson:
- How to incorporate the data of Prigent et al. (2005) in the surfdata generation. I will write an email about this to Catherine Prigent.
- One statement marked in CanopyFluxesMod should probably be changed when using Meier2022 (i.e., Yang et al. (2008) formulation should be used instead of Zeng and Dickinson (1998)).
- At the moment one needs to change the parameter file to switch between the original and proposed globally constant z0m values for bare soil, snow, and ice. This is obviously not very user friendly and prone for mistakes.
- The introduction of Yang et al. (2008) frequently results in z0h and z0q larger than z0m. This is only rarely observed in the field and in contradiction to the theory that z0h and z0q should be smaller because heat and water vapor need to be transported through molecular diffusion in the surface sublayer.
Issues Resolved: #1316
I know the date for your paper isn't completely known right now (hence the XXXX in the namelist name). But, it seems like we'll want a real year when this comes into main-dev of CTSM. So I propose we go with Meir2022? I think in the context of namelist options for CTSM being off a bit in the namelist name is acceptable. Especially if the exact reference is pointed to in the code. And we should put the reference in appropriate places where it's refereed to.
Thanks @ekluzek for reviewing the code. I don't have any strong opinions regarding the tag. Meier2022 should most likely be appropriate. We can also go for MeierDavin2022 or GMD2022. Not sure if CESM has any guidlines on how to choose such tags.
The commented-out code sections marked with "! --> Use this for CLM-YY" were used for the sensitivity tests that are presented in the appendix of the GMD study. Those can all be removed (sorry I forgot to do this). What could be considered to be included as an additional option is the formulation of Owen and Thomson (1963) to compute z0hg and z0qg.
Also, I realized that we actually opened an issue for this development a couple of months ago: New parameterization of vegetation surface roughness for forests #1316
We talked about this a bit this morning.
One aspect that we should keep in mind about this is how this interacts with FATES. Our testing ensures we won't break FATES with it. Two questions are can part of it be turned with FATES, and secondly can/should FATES be updated to bring these kind of changes to the FATES side of things? Before we do that I'd like to make a little more progress so the code is a little more clearer for a FATES scientist to evaluate.
Specifically, the next steps for this are:
- Fix the subscript overflow in DEBUG mode
- Change the MeierXXXX namelist tag to Meier2022
- Tag that as a branch tag
- Remove the Local time changes
- Remove the PFT on separate columns option
- Update to the latest model version
- Make sure all model configurations work with the new option (or disable using them together)
- Run more extensive testing to make sure we agree this is good to come in
- We may remove the 2D option that requires changes to surface datasets
- Some of the parameters should be made into parameters or placed on the paramsfile
- Possibly have it turned on by default for clm5_1 physics
Note that SMS_D_Ly2.f10_f10_mg37.I1850Clm51SpNoAnthro.cheyenne_intel.clm-Meier2022_surf_rough is failing at startup because "Forcing height is below canopy height for patch index".
@olyson and @swensosc I'm running into some science issues with this branch that I'll need you to look into. I think I should go ahead with our plan to bring it up to the latest tag, and then pass it over to one or both of you to look into. Just FYI for now, unless you or others think you should start looking into it now. I'm adding the "next" label so that we'll talk about this Thursday.
I've got most of the updates that needed to happen that were talked about above. It's updated to ctsm5.1.dev077 and other changes removed and the namelist option is Meier2022 now. There are some tests that fail though, including some that don't turn the new changes on.
ERP_D_Ld5_P48x1.f10_f10_mg37.I2000Clm50BgcCru.izumi_nag.clm-noFUN_flexCN ERP_D_Ld5_P48x1.f10_f10_mg37.I2000Clm50BgcCru.izumi_nag.clm-reduceOutput SMS_D_Ld10.f10_f10_mg37.I2000Clm50BgcCrop.izumi_intel.clm-tracer_consistency SMS_D_Ld5.f10_f10_mg37.I2000Clm50BgcCrop.izumi_nag.clm-irrig_alternate ERP_D_Ly3.f09_f09_mg17.I2000Clm51Bgc.cheyenne_intel.clm-Meier2022_surf_rough_all_f09NonCrop ERS_Ly3.f10_f10_mg37.I2000Clm51BgcCrop.cheyenne_intel.clm-Meier2022_surf_rough LWISO_Ld10.f10_f10_mg37.I2000Clm50BgcCrop.cheyenne_gnu.clm-coldStart SMS_D.f10_f10_mg37.I2000Clm51Bgc.cheyenne_gnu.clm-Meier2022_surf_rough SMS_D.f10_f10_mg37.I2000Clm51Bgc.cheyenne_intel.clm-Meier2022_surf_rough SMS_D_Ln1.f09_g17.I1850Clm50BgcNoAnthro.cheyenne_intel.clm-Meier2022_surf_rough
The surf_rough tests fail on izumi, because the files are only on cheyenne right now.
OK, one simple thing I needed to do was update the Ronny Meier paramsfile to the latest, and I did that so that fixes some problems. But, now all of the tests are failing with the following problem. I saw this with the testing in ctsm5.1.dev029, but I didn't see all of the tests failing.
Error: Forcing height is below canopy height for patch index iam = 0: local patch index = 13 iam = 0: global patch index = 1935 iam = 0: global column index = 1366 iam = 0: global landunit index = 457 iam = 0: global gridcell index = 145 iam = 0: gridcell longitude = 15.0000000 iam = 0: gridcell latitude = 30.0000000 iam = 0: pft type = 5 iam = 0: column type = 1 iam = 0: landunit type = 1 ENDRUN: ERROR: ERROR in /glade/work/erik/ctsm_worktrees/surf_rough/src/biogeophys/CanopyFluxesMod.F90 at line 975
I met with @olyson and he is going to verify my changes as well as look into some of the problems I ran into. There could be a mistake in my changes to prevent floating point overflow. Or there could be something with the updated code that's working differently than the version that Ronny ran with.
Regarding the forcing height error, I think this statement in biogeophys/BiogeophysPreFluxCalcsMod.F90:
if ( is_first_step() .or. get_nstep() <= GetBalanceCheckSkipSteps() ) then
should be:
if ( is_first_step() .or. get_nstep() <= GetBalanceCheckSkipSteps()-1 ) then
because at that point htop is available and displacement height can be calculated. Otherwise, the forcing height will be calculated as
forc_hgt_u_patch(p) = forc_hgt_u(g) + z0mg(c) + displa(p) = forc_hgt_u_patch(p) = 30m, since z0mg and displa are zero.
Then displa is recalculated in CanopyFluxes and that calculation can be greater than 30m and hence zldis(p) = forc_hgt_u_patch(p) - displa(p) is negative.
Tried that and the forcing height errors went away.
@olyson awesome. Thanks for figuring that out. I'll add that in and make sure my tests work now.
This check you put in is still triggering, so I was going to look at that next:
ustar*thvstart is positive and has to be negative
More generally, there is inconsistency between the forcing height which is based on this calculation in BiogeophysPreFluxCalcsMod.F90:
displa(p) = htop(p) * (1._r8 - (1._r8 - exp(-(7.5_r8 * (pftcon%z0v_LAImax(patch%itype(p))))*0.5_r8)) & / (7.5_r8(pftcon%z0v_LAImax(patch%itype(p)) ))**0.5_r8)
and the later calculation of displa in CanopyFluxes:
displa(p) = htop(p) * (1._r8 - (1._r8 - exp(-(7.5_r8 * lt)**0.5_r8)) / (7.5_r8*lt)**0.5_r8)
So I think we should be updating the forcing height for this new displa in CanopyFluxes.
Tried this and it fixed the negative ustars in my test that were occurring in nsteps 0,1,2. Error was occurring because zldis, while positive, was less than the roughness length, so ln(zldis/z0) was negative.
More generally, there is inconsistency between the forcing height which is based on this calculation in BiogeophysPreFluxCalcsMod.F90:
displa(p) = htop(p) * (1._r8 - (1._r8 - exp(-(7.5_r8 * (pftcon%z0v_LAImax(patch%itype(p))))*0.5_r8)) & / (7.5_r8(pftcon%z0v_LAImax(patch%itype(p)) ))**0.5_r8)
and the later calculation of displa in CanopyFluxes:
displa(p) = htop(p) * (1._r8 - (1._r8 - exp(-(7.5_r8 * lt)**0.5_r8)) / (7.5_r8*lt)**0.5_r8)
So I think we should be updating the forcing height for this new displa in CanopyFluxes.
Hi,
I think that's a quick one to answer: Those displa and z0mv values are a first guess only. In the default version of the model, they are computed as if the LAI+SAI went towards infinity (i.e., the highest possible values with this vegetation height and type). With the proposed formulations, z0mv has a maximum at LAImax. Therefore, I used this value for a first guess, analogues to the default implementation, as indicated by the comment "Compute as if elai+esai = LAImax - LAIoff in CanopyFluxes". Note that LAIoff is now 0. So, the "- LAIoff" could be removed in the comment.
I hope that makes sense, Ronny
@olyson to check into the branch we would need to get you collaborator access. We can ask @RonnyMeier to do that for you. But, in the meantime if you have proposed changes that you can show me, I can check them in and try it out.
You can try:
/glade/work/oleson/ctsm_MeierZ0/src/biogeophys/BiogeophysPreFluxCalcsMod.F90 /glade/work/oleson/ctsm_MeierZ0/src/biogeophys/CanopyFluxesMod.F90
The cheyenne tests you listed above all pass now except for:
LWISO_Ld10.f10_f10_mg37.I2000Clm50BgcCrop.cheyenne_gnu.clm-coldStart SMS_D_Ln1.f09_g17.I1850Clm50BgcNoAnthro.cheyenne_intel.clm-Meier2022_surf_rough
The first doesn't have Meier2022 on and I haven't found any useful error messages for the second test yet.
Actually, what I have now in /glade/work/oleson/ctsm_MeierZ0/src/biogeophys/CanopyFluxesMod.F90 makes more sense, i.e.,
z0hv(p) = z0mv(p)
z0qv(p) = z0mv(p)
+ ! Update the forcing heights
+ 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 do
As @ekluzek noted, the tests were also failing due to negative humidity in HumanIndexMod, so he turned that off by default. I've turned that on now in my latest version of the code and the tests are still failing. In this particular test:
SMS_D.f10_f10_mg37.I2000Clm51Bgc.cheyenne_intel.clm-Meier2022_surf_rough
there is one pft that is triggering the error. CanopyFluxes is producing a slightly negative 2-m specific humidity which is then converted to relative humidity and then fed into HumanIndexMod. This seems to be occurring under very stable conditions, with tall trees with large elai+esai, a displacement height near the canopy top and a very small ustar (numbers below):
zeta = 100 z0h = z0q = 4.29 htop = 36 displa = 31.3 ustar = 0.0114 elai+esai = 7.1
With the stress indices turned off, the test runs to completion even with the negative humidities. I tried restricting the max zeta to 60 (instead of 100) and that eliminated the negative humidities, at least for that test. Not a real solution though.
Regarding this test that was failing (despite not having Meier2022 on):
LWISO_Ld10.f10_f10_mg37.I2000Clm50BgcCrop.cheyenne_gnu.clm-coldStart
This was failing with:
59:ERROR in CompareBulkToTracer: tracer does not agree with bulk water 59:Called from: after first stage of hydrology 59:Variable: qflx_snomelt_accum_col 59:First difference at index: 15 59:Bulk : 0.94434303138415165E-04 59:Tracer: 0.98722208432612920E-04 59:ratio: 0.10000000000000000E+01 59:Bulk*ratio: 0.94434303138415165E-04
This is a new variable that was added to WaterFluxType.F90 to be used in the new soil roughness calculation. After consultation with @billsacks , I've moved this to WaterDiagnosticBulkType since it is a diagnostic variable (diagnostic in the sense that it isn't a fundamental state variable) and only used in the roughness calculation. We don't need any tracer/isotope quantities associated with this. I've also renamed it from qflx_snomelt_accum to snowmelt_accum since it isn't a flux. That test passes now.
This test:
SMS_D_Ln1.f09_g17.I1850Clm50BgcNoAnthro.cheyenne_intel.clm-Meier2022_surf_rough
is failing because of:
20220310 182753.068 ERROR PET0205 src/addon/NUOPC/src/NUOPC_Base.F90:950 Invalid argument - setClock timeStep=10800s is not a divisor of runDuration=1800s
Found this in the PET log files. This seems to be a 1 time-step run. Seems unrelated to the roughness parameterization...
This test:
SMS_D_Ln1.f09_g17.I1850Clm50BgcNoAnthro.cheyenne_intel.clm-Meier2022_surf_rough
is failing because of:
20220310 182753.068 ERROR PET0205 src/addon/NUOPC/src/NUOPC_Base.F90:950 Invalid argument - setClock timeStep=10800s is not a divisor of runDuration=1800s
If you are using an active runoff model with standard coupling interval of 3 hours, then the test needs to run for at least 3 hours.
Is there a reason why this test needs to be at f09 resolution? I'm in the process of trying to reduce our high resolution testing because it is expensive to run these high-resolution tests with every run of the test suite, and they can occasionally hold up testing by getting stuck in the queue. Making them _Ln1 doesn't really help because most of the time in these short high-res tests is spent in model initialization, so whether it's one time step or a few days doesn't really matter much.
Thanks, that explanation makes sense. I didn't create these tests so I'll leave it up to @ekluzek to respond regarding the question about resolution.
@ekluzek it seems like the ZengWang2007 version of z0m in FrictionVelocity needs some protection for log10 as is done for Meier2022.
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
else
z0mg(c) = this%zsno
end if
Add an if statement:
if ( snomelt_accum(c) < 1.e-5_r8 )then
z0mg(c) = exp(1.4_r8 * -rpi/2.0_r8 -0.31_r8) / 1000._r8
Another question is whether we should even be using that expression for z0mg for ZengWang2007. It used to be just zsno right?
@billsacks the noAnthro test needs to be f09 because it's the only resolution for Potential Vegetation surface dataset. But, for now I'm just adding more tests than are really needed while we are working on the code. I don't think I'll keep that test when it comes to master. So I'll add a note to that effect for some of the tests that should be removed. One thing that I really want to do is to run the whole test suite with the Meier surface roughness on for everything to make sure everything works.
Anyway, @olyson I'll work with this test and add a note to it that we won't keep it long term.
@olyson you probably have more code that could be checked in. I asked @RonnyMeier to add collaborator access, so you should likely checkin your latest changes once you get that.
@ekluzek , is there a reason why the accumulated snow melt parameterization and the roughness length dataset are also being applied to the ZengWang2007 parameterization?
@olyson no, that shouldn't happen. That would be a mistake in the code. We want to preserve answers when the new code isn't being used. I did show that it is for a lot of tests, but if there are places where the new code affects the previous parameterization we should remove that.
I think the one remaining problem now is the negative 2-m humidity. I can solve that by limiting zeta in stable conditions, but not sure if we really want to do that or if there is another way to solve the problem. I think we need to consult with some of the group at this point.
@olyson could you call a meeting with the relevant people that should be involved in discussing that last problem?
You had also asked me to check the expression for z0mg when snomelt_accum = 0. I think your formulation is correct.