Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Surface roughness modifications (cont'd) #2045

Merged
merged 72 commits into from
Aug 23, 2023

Conversation

slevis-lmwg
Copy link
Contributor

@slevis-lmwg slevis-lmwg commented Jun 28, 2023

Description of changes

This PR continues from where #1596 left off.

Specific notes

Contributors other than yourself, if any:
@ekluzek @RonnyMeier

CTSM Issues Fixed (include github issue #):
Fixes #1596
Fixes #1316

Are answers expected to change (and if so in what way)?
Yes

Any User Interface Changes (namelist or namelist defaults changes)?
Yes, NLCOMP diffs in the test-suites

Testing performed, if any:
Test-suite results reported below

Ronny Meier and others added 30 commits March 25, 2021 08:37
…ime output (the latter is still in development). This version of the code is running stable.
…ghness as htop isn't set yet. Also use a different threshold for snomelt_accum suggested by Ronny Meier
… of a negative number being done, in this case set wc to zero and continue
…ss is used to prevent a numerical error that kills the run in calc of WetBulb temp because of negative relative humidity
Resolved conflicts:
doc/ChangeLog
doc/ChangeSum
@slevis-lmwg

This comment was marked as outdated.

Copy link
Collaborator

@ekluzek ekluzek left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, awesome. I have some questions and things to do, more parameters mostly. Some of this could be made into an issue and done later.

We'll decide on this when we talk it over together.

select case (z0param_method)
case ('ZengWang2007')
if (frac_sno(c) > 0._r8) then
if(use_z0m_snowmelt) then
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There was some discussion on this in #1596. The question I wonder about is this as presented adds the "use_z0m_snowmelt option for ZengWang2007 (and use_z0m_snowmelt is a new option for this new work). So I wonder if we want to enable this option for the old structure? I kind of think we should only let it work on the new Meir2022 option.

The discussion in question is here: #1596 (comment)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For code maintenance unless this option is needed for science we'd like to simplify it and not add this option for ZengWang.

@olyson or @swensosc do you have opinions on how useful this would be for running with ZengWang?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we simplify though we should change the buildnml so that you can NOT use this option with ZengWang.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not aware of any need for this. Might be better to simplify.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed

if ( snomelt_accum(c) < 1.e-5_r8 )then
z0mg(c) = exp(-b1_param * rpi * 0.5_r8 + b4_param) * 1.e-3_r8
else
z0mg(c) = exp(b1_param * (atan((log10(snomelt_accum(c))+0.23_r8)/0.08_r8)) + b4_param) * 1.e-3_r8
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The 0.23_r8 is repeated 4 times here and should be made into a parameter either in the code or params file.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also the 0.08

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Replaced


case ('Meier2022')
lt = max(0.00001_r8,elai(p)+esai(p))
displa(p) = htop(p) * (1._r8 - (1._r8 - exp(-(cd1_param * lt)**0.5_r8)) / (cd1_param*lt)**0.5_r8)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wonder if these **0.5 terms should be changed to sqrt? Although that will likely change answers. So maybe we make it an issue and do it as an answer changing tag?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's not do this. But, possibly add some comments about it....

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Addressed with comments

*z0v_c(patch%itype(p)) * lt * 0.25_r8
U_ustar = U_ustar_ini

do while (delt > 0.0001_r8)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There's this 0.0001 here and the 0.00001 above that I wonder if they should be parameters? Although they don't seem to be repeated....

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since, they aren't repeated we think we'll leave them alone.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This expression is confusing though. For example the **(-0.5) which means 1/sqrt(...).

We think adding a comment about this would be helpful.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Replaced

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm, I didn't realize that they were different, and made them both 1.e-4_r8. I should probably back this change out...

@@ -214,6 +214,8 @@ attributes from the config_cache.xml file (with keys converted to upper-case).
<zetamaxstable phys="clm4_5" >2.0d00</zetamaxstable>
<zetamaxstable phys="clm5_0" >0.5d00</zetamaxstable>
<zetamaxstable phys="clm5_1" >0.5d00</zetamaxstable>
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we should remove this line, as the user should have to choose between ZengWang and Meier. I'm not sure it matters, but I think it's clearer with it removed...

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should make sure z0param_method HAS to be selected in the buildnml.

case ('Meier2022')

! After Yang et al. (2008)
z0hg_patch(p) = 70._r8 * nu_param / ustar(p) * exp( -beta_param * ustar(p)**(0.5_r8) * (abs(tstar))**(0.25_r8))
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The 70 is repeated and should become a parameter. Actually I think this expression is repeated and maybe should become a function. This come later...

Also the **0.5...

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We think we should leave the **0.5 and **0.25 should stay as they are. The **0.25 should probably stay that way as that's better than sqrt(sqrt(...)), so the other one should as well.

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
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This 0.61 is somewhat outside the scope here, but should become a parameter. We should maybe make an issue for that.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@@ -1014,6 +1014,15 @@ subroutine Wet_BulbS (Tc_6,rh,wbt)
! !LOCAL VARIABLES:
!EOP
!
#ifndef NDEBUG
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we really want this #ifdef here?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This was on purpose. But, I'm not sure if really should be done.

Was it expensive?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This was on purpose. But, I'm not sure if really should be done (the #if...)

Was it expensive? No not more expensive than the expression below which is expensive, and can blow up for bad values.

And this WetBulb stuff is already known to be expensive, hence we don't turn it on that often.

So final bottom-line remove the #ifndef and the #endif

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Gone!


! After Yang et al. (2008)
z0hg_patch(p) = 70._r8 * nu_param / ustar(p) * exp( -beta_param * ustar(p)**(0.5_r8) * (abs(tstar))**(0.25_r8))

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's remove this commented out code. We can't maintain it when commented out. And without scientific effort we can't add it, as we don't know if it really offers scientific value.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed

@@ -1007,7 +1066,12 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp,
! changed by K.Sakaguchi from here
! transfer coefficient over bare soil is changed to a local variable
! just for readability of the code (from line 680)
csoilb = vkc / (params_inst%a_coef * (z0mg(c) * uaf(p) / 1.5e-5_r8)**params_inst%a_exp)
! RM: Does this need to be updated if Ya08 is used too? Proposed formulation (definitely double-check!)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove this commented out alternative, that we can't maintain.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed

! Surface roughness parameterization
!----------------------------------------------------------

character(len=64), public :: z0param_method
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should have a comment...

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Added

@ekluzek
Copy link
Collaborator

ekluzek commented Aug 17, 2023

This has more work to be done on it, so won't be done immediately. But, we'll schedule it for next week....

@ekluzek
Copy link
Collaborator

ekluzek commented Aug 17, 2023

When I do a diff I see a lot of red for white-space at the ends of lines. It would be good to remove that for the new code coming in. The red makes me hate that white-space at the end.

@slevis-lmwg
Copy link
Contributor Author

slevis-lmwg commented Aug 18, 2023

Addressed most suggestions and now running two preliminary tests before submitting test-suites:
PASS ./create_test ERS_Ly3.f10_f10_mg37.I2000Clm51BgcCrop.cheyenne_intel.clm-Meier2022_surf_rough
PASS ./create_test ERP_P72x2_Lm25.f10_f10_mg37.I2000Clm51BgcCrop.cheyenne_intel.clm-monthly -c /glade/p/cgd/tss/ctsm_baselines/ctsm5.1.dev134

@slevis-lmwg
Copy link
Contributor Author

slevis-lmwg commented Aug 21, 2023

Test-suites comparing against dev134: OK

@slevis-lmwg slevis-lmwg merged commit 0e7457f into ESCOMP:master Aug 23, 2023
2 checks passed
@slevis-lmwg slevis-lmwg deleted the meier_main_pr1596 branch August 23, 2023 17:59
@apcraig
Copy link

apcraig commented Aug 24, 2023

I'd like to ask a question about this PR.

I am adding CTSM to RASM (CESM1.2). I'm currently doing so with release-clm5.0.35, which I think is the last release from several months ago. Anyway, we need to couple the (logz0) to WRF. Several years I did the same with an earlier version of CTSM (for the NNA project, still ongoing), but I'm trying to shift to a more recent and robust release. In the older version, there was a subroutine called SetActualRoughnessLengths which computed z0's from a call in clm_drv. In release-clm5.0.35, that subroutine is gone and it looks like z0's are only computed in Fates, if Fates is on (just from looking at the code). Is that correct?

In our test run, use_fates=.false. and it looks like z0's are never computed. Several questions.

  • Does this PR and recent changes address this shortcoming? Are z0's now computed when fates is off?
  • Should I be running with fates on?
  • More generally, would you recommend another version more recent than release-clm5.0.35 for us to use? I prefer to use a more robust version than some recent hash in master, but there are trade-offs.

Thanks.

@slevis-lmwg
Copy link
Contributor Author

Hi @apcraig ,

My interpretation of this PR's mods:
This PR modified how z0 terms were calculated (eg see description of changes in the ChangeLog ). So z0 terms were already being calculated a different way before this PR (eg see changes in BaregroundFluxes, BiogeophysPreFluxCalcs and elsewhere). Looking at this PR's mods may help you more easily find where z0 is calculated. This PR does not turn on the new z0 calculations by default, so it keeps answers unchanged from the last baseline.

You should have access to z0 calculations with or without Fates. I do not know why Fates does its own calculation of z0. Maybe it was easier (or better) than using the clm calculation.

As you suggest, it may not make sense to work with this most recent version UNLESS you wish to use the new z0 calculations (z0param_method = "Meier2022").

@apcraig
Copy link

apcraig commented Aug 25, 2023

Thanks @slevis-lmwg. OK, I'll have a look at this version. I think an important point is that in release-clm5.0.35, z0 is only calculated by fates. At some point since that version, it sounds like this was updated and z0 is calculated with or without fates. It sounds like I don't need the latest version, but I need some version since release-clm5.0.35. Can you recommend a particular ctsm tags or hash that is relatively robust and tested that might meet my needs, I assume something newer than release-clm5.0.35 and something older than the current master. Thanks!

@slevis-lmwg
Copy link
Contributor Author

When I type git grep z0 in a copy of clm5.0.00 that I have, I see calculations of z0, so I wonder why release-clm5.0.35 is different. I don't feel I have followed the model development closely enough to be able to recommend the right version for your work.

@samsrabin samsrabin added enhancement new capability or improved behavior of existing capability science Enhancement to or bug impacting science labels Aug 8, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement new capability or improved behavior of existing capability science Enhancement to or bug impacting science
Projects
Status: Ready to eat (Done!)
Status: Done (non release/external)
Development

Successfully merging this pull request may close these issues.

New parameterization of vegetation surface roughness for forests
5 participants