Soil layer definition clean-up and user-defined option

[slevis-lmwg]

I currently envision the user-defined option getting folded into the clean-up work because I do not expect extensive code changes for either; of course I'm open to different outcomes.

Description of changes

1. Clean-up following @swensosc and @billsacks suggestions in issue Clean up soil vertical layers definition #279

2. Initial proposal for a user-defined option of soil layer structure:

• Use existing namelist variable soil_layerstruct
• Instead of setting to an existing option (e.g. '5SL_3m', etc.),
  set soil_layerstruct = 'user:aa,bbbb,cccc,dddd,eeee,ffff,gggg,...'
  where the string is declared of len=256 and
  'user:' tells the code to expect a user defined structure,
  aa gets assigned to nlevsoi and nlevgrnd in the code,
  bbbb gets assigned to dzsoi(1) in the code,
  cccc gets assigned to dzsoi(2) in the code,
  and so on to dzsoi(nlevgrnd).

This user-defined option should likely be tested with a new (or modified existing) test in the test-suite.

Specific notes

Contributors other than yourself, if any:
@billsacks @dlawrenncar @swensosc @ekluzek

CTSM Issues Fixed (include github issue #):
fixes #279
fixes #728

Are answers expected to change (and if so in what way)?
For NWP cases, because had to change nlevsoi from 5 to 4 to work around the current limitation that nlevsoi must be less than nlevgrnd for the model to not abort.

Any User Interface Changes (namelist or namelist defaults changes)?
New user-defined option to specify soil layer structure.

Testing performed, if any:
This far I have only tried one test (it PASSed) to make sure that a pre-existing soil layer structure works as it did before the clean-up:
./create_test ERS_D_Ld10.f10_f10_musgs.I2000Clm50BgcCropGs.cheyenne_intel.clm-rm_indiv_lunits_and_collapse_to_dom -c /glade/p/cgd/tss/ctsm_baselines/ctsm1.0.dev049

[slevis-lmwg]

I need to change this to read 4 characters per dzsoi(j) since the units are meters and 4 characters will allow
min(dzsoi(j)) = .001 m

[billsacks]

@slevisconsulting thanks a lot for starting to work on this - this is going to be a big help!

As for the user interface: I feel it would be better to use a new vector-valued namelist variable (or multiple variables, if necessary) for this, rather than trying to do this in the existing string variable. I feel this will be more intuitive for users, make error-checking easier, and it will be easier for us to write and maintain the code.

I think the first thing to decide is: for our out-of-the-box cases, do we want to continue to use a string that gives some shorthand name for a hard-coded layer structure, or do we want to switch to always having a list of dzsoi values in the namelist? The latter will be easier (and argues even more strongly for doing this with a new variable that is a vector of reals rather than maintaining the old string value), but my suggestion should be doable even if we go with the former. I'll elaborate on how I envision the former working:

We have two namelist variables; a string (say, soil_layerstruct_predefined) and a vector of reals (say, soil_layerstruct_userdefined). For a given run, exactly one of them would appear in lnd_in. The hard-coded defaults are 'UNSET' and a vector of spvals (or some other value that would never appear in practice, like -999). If both are specified in user_nl_clm, that's an error (caught by build-namelist). The Fortran code would then look like this:

  if (soil_layerstruct_userdefined(1) == spval .and. soil_layerstruct_predefined == 'UNSET') then
     call endrun("Must set one")
  else if (soil_layerstruct_userdefined(1) /= spval .and. soil_layerstruct_predefined /= 'UNSET') then
     call endrun("Can't set both")
  end if

  if (soil_layerstruct_userdefined(1) /= spval) then
     nlevgrnd = ...  ! set based on last value that is NOT spval
     ! Here, set various things based on the userdefined values
  else
     ! Here, set various things based on the predefined string
  end if

  ! Here, set things common to both methods

Does that make sense, or do you see problems with that approach?

[slevis-lmwg]

@billsacks this makes sense to me.

[slevis-lmwg]

Early tests have revealed a couple of preexisting bugs. I will outline them here. Let me know your preferred way to proceed.

This test originally PASSed:
ERS_D_Ld10.f10_f10_musgs.I2000Clm50BgcCropGs.cheyenne_intel.clm-rm_indiv_lunits_and_collapse_to_dom -c /glade/p/cgd/tss/ctsm_baselines/ctsm1.0.dev049

But with soil_layerstruct = '5SL_3m' the test FAILed with this error
Subscript 2 of the array DZ has value 6 which is greater than the upper bound of 5
in LakeTemperatureMod.F90 in a loop with use_lch4 = .true.
I replaced dz with dz_lake in line 960 and got past this error.

Next, the same test FAILed with
Subscript 2 of the array SOM_ADV_COEF has value 6 which is > than the upper bound of 5
in SoilBiogeochemLittVertTranspMod.F90
I included "+1" in the next four statements in SoilBiogeochemStateType.F90
allocate(this%som_adv_coef_col (begc:endc,1:nlevdecomp_full+1))
allocate(this%som_diffus_coef_col (begc:endc,1:nlevdecomp_full+1))
this%som_adv_coef_col(c,1:nlevdecomp_full+1) = 0._r8
this%som_diffus_coef_col(c,1:nlevdecomp_full+1) = 0._r8
and got past this error.

Next, the same test FAILed with
Subscript 1 of the array ZSOI has value 6 which is greater than the upper bound of 5
in SoilBiogeochemLittVertTranspMod.F90
The indexing problem is pervasive here and manifests when nlevgrnd = nlevsoi because the code here assumes nlevgrnd > nlevsoi. This is more complex to resolve, and we can install error checks that inform users of this requirement.

[billsacks]

So it sounds like the 5-layer soil can't be used with BGC currently? I can't tell from your description if the fixes you put in for (1) and (2) are correct fixes (as far as you can tell) or just temporary workarounds to get the test to pass; also, would those fixes change answers for standard cases?

Maybe we can discuss this further at the meeting tomorrow?

[slevis-lmwg]

I will respond. We can also discuss further at tomorrow's meeting:

1. Yes, the 5-layer soil cannot be used with BGC currently.
2. My first fix (dz --> dz_lake) is the correct, as far as I can tell, fix of the first bug; it changes answers from dev049.
3. My second fix is a partial fix (correct, I think) to the larger problem that I described with SoilBiogeochemLittVertTranspMod.F90. Due to the complexity of the larger problem, I backed out the second fix and re-tested the 5-layer soil with a 2-line change that worked:
   nlevgrnd = 6
dzsoi(6) = 1
   ...and a 1-line change that also worked; here I also confirmed that (as I expected) answers change from dev049:
   nlevsoi = 4

[billsacks]

From today's CTSM software meeting:

In the short-term, to fix the problem @slevisconsulting ran into: Set nlevsoi = 4 for the NWP case. We should then rename 5SL_3m to 4SL_2m.

What do we want to do for out-of-the-box cases? Use string or have the new vector of values? For now, maybe stick with string; we may decide to change that later.

The user interface for setting a user-defined soil layer structure should be:

• Specify a vector of real values giving dz for each layer; the length of this vector gives nlevgrnd
• If a user specifies that, then they also need to set a second namelist value that gives the max number of soil layers (i.e., nlevsoi). For now, we require that this be <= (length of dz vector) - 1. (Eventually the requirement will just be that it be <= (length of dz vector).)

[slevis-lmwg]

Testing to date:
./create_test ERS_D_Ld10.f10_f10_musgs.I2000Clm50BgcCropGs.cheyenne_intel.clm-rm_indiv_lunits_and_collapse_to_dom -c /glade/p/cgd/tss/ctsm_baselines/ctsm1.0.dev049
PASS when I remove the dz → dz_lake fix in LakeTemperatureMod.F90, else not bfb w dev049

Rerun w/o comparing to baseline, with the dz → dz_lake fix, and with soil_layerstruct_predefined = '5SL_3m' PASS

Rerun with soil_layerstruct_userdefined = 0.01,0.02,0.04,0.08,0.16,0.32,0.64,1.28,2.56,5.12 and with nlevsoinl = 9 PASS

[billsacks]

@slevisconsulting - Thanks for your updates! I haven't looked closely yet, but at a glance this looks good.

I was thinking about what testing we could / should have for this. I can think of three possibilities:

(1) A few unit tests of the initialization logic. This would be cleanest if we could combine the logic for initializating nlevsoi and nlevgrnd with the other logic in initVertical, but I'm not sure if that's feasible. Otherwise, the unit tests could call clm_varpar_init followed by clm_varcon_init (to allocate variables) followed by initVertical. Ideally, to test all of the logic, we'd also do the namelist reads, but I think that's awkward right now in a unit test, so instead the unit tests could start by setting the appropriate variables that are typically read from namelists, and then call the above routines.

(2) Introduce a new test type that runs two cases: one with a predefined soil layer structure and one with a user-defined soil layer structure that exactly matches the predefined one; ensure the two are bit-for-bit. This should be quite easy to do. (Caveat: the namelist settings probably need to look like 0.1d0 in order to give double precision that matches the user-defined settings.)

(3) Make sure this new option is covered by some system test, but don't actually use the new test type of (2).

I feel like (1) would be ideal, both to give confidence that the settings are being made correctly now and in the future and to reduce our current over-reliance on system tests. However, it would also take the most work to put in place. With (1), I personally wouldn't feel a need for an additional system test. If we do have a system test, (2) would give the most confidence, but it would also contribute to our system test bloat. Overall, I don't have strong feelings on what should be done, though it seems we should probably have some testing of this.

[slevis-lmwg]

Thank you for these suggestions @billsacks. I'm drawn to option (3), partly because it's easiest for me to implement and partly because I suspect that we could go from (3) to (2) rather easily (with your guidance because right now I would not know how).

Regarding the "system test bloat" I wonder whether one of the NWP tests could merge with my "rm_indiv_lunits_and_collapse_to_dom" test at some point. I don't think they need to be 2 tests.

When I'm ready to run the test-suite, I will need to update to the next tag to get the #760 #761 bug-fix.

[billsacks]

Regarding the "system test bloat" I wonder whether one of the NWP tests could merge with my "rm_indiv_lunits_and_collapse_to_dom" test at some point. I don't think they need to be 2 tests.

Good point. Do you want to go ahead and make this change, folded in with this tag or some other? I see that there are currently 3 NWP tests:

aux_clm: ERP_P36x2_D_Ld5.f10_f10_musgs.I2000Ctsm50NwpSpGswpGs.cheyenne_intel.clm-default                      # Include a debug ERP test of the NWP configuration.
aux_clm: SMS_Ld1.nldas2_rnldas2_mnldas2.I2000Ctsm50NwpSpNldasGs.cheyenne_gnu.clm-default                      # Include a short smoke test covering the nldas2 grid and the I2000Ctsm50NwpSpNldasGs compset, which uses NLDAS datm forcing.
aux_clm: SMS_Ld1.nldas2_rnldas2_mnldas2.I2000Ctsm50NwpSpNldasRsGs.cheyenne_gnu.clm-default                    # Include a short smoke test covering the nldas2 grid and the I2000Ctsm50NwpSpNldasRsGs compset, which uses NLDAS datm forcing.

I'd like to keep at least one that uses the default, out-of-the-box NWP configuration, but one of the above could be changed to include rm_indiv_lunits_and_collapse_to_dom. NWP already includes some aspects of rm_indiv_lunits_and_collapse_to_dom, but not all. I don't have a good sense of whether the aspects of rm_indiv_lunits_and_collapse_to_dom that are not included by NWP out-of-the-box might already be sufficiently covered by unit test; if so, we could remove rm_indiv_lunits_and_collapse_to_dom entirely. I'll leave that up to you.

[slevis-lmwg]

@billsacks I'm revisiting this question:

What do we want to do for out-of-the-box cases? Use string or have the new vector of values? For now, maybe stick with string; we may decide to change that later.

...because I've encountered the following conundrum in the test suite about keeping answers bfb:

The code used to say...
<soil_layerstruct structure="fast" >5SL_3m</soil_layerstruct>
<soil_layerstruct structure="standard" phys="clm5_0">20SL_8.5m</soil_layerstruct>
<soil_layerstruct structure="standard" phys="clm4_5">10SL_3.5m</soil_layerstruct>
and now says...
<soil_layerstruct_predefined structure="fast">4SL_2m</soil_layerstruct_predefined>
<soil_layerstruct_predefined structure="standard" phys="clm5_0">UNSET</soil_layerstruct_predefined>
<soil_layerstruct_predefined structure="standard" phys="clm4_5">10SL_3.5m</soil_layerstruct_predefined>

I went to .../cime_config/testdefs/testmods_dirs/clm to add soil_layerstruct_predefined definitions in all the user_nl_clm files and remembered that many (most?) clm4.5 tests share the same user_nl_clm as their clm5.0 counterparts. So adding soil_layerstruct_predefined = '20SL_8.5m' would change clm4.5 answers. Replacing these strings with vectors of values will allow different defaults for clm5 and clm4.5 using the code quoted above. We may still lose bfb because replacing (e.g.) this formula
dzsoi(j) = dzsoi(nlevsoi) + (((j - nlevsoi) * 25._r8)**1.5_r8) / 100._r8
with a vector of values may not be exact.

I'm interested in your (and others') opinion.

[billsacks]

I may be confused here: why did you need to change the namelist_defaults for structure="standard"? The feeling from last week's meeting was that we should stick with a string for the out-of-the-box cases.

[slevis-lmwg]

I changed it to UNSET to ensure that the user was aware that they had to set at least one, either the predefined or the userdefined option.

[billsacks]

I don't think that's what we want: we want the out-of-the-box behavior to be the same as before - otherwise users will always have to explicitly set this, which will be a pain.

[billsacks]

Ah, but I think I see how my earlier suggestions may have been wrong, then. Maybe we need some logic in build-namelist that sets the predefined option to UNSET if the user-defined one is set???

Then the logic in the Fortran could stay as I suggested earlier - i.e., ensuring that only one of these is set.

[slevis-lmwg]

Ok, I'll do that.

[slevis-lmwg]

To combine
ERP_P36x2_D_Ld5.f10_f10_musgs.I2000Ctsm50NwpSpGswpGs.cheyenne_intel.clm-default
with
ERS_D_Ld10.f10_f10_musgs.I2000Clm50BgcCropGs.cheyenne_intel.clm-rm_indiv_lunits_and_collapse_to_dom
I propose replacing the two with this test
ERP_P36x2_D_Ld5.f10_f10_musgs.I2000Ctsm50NwpBgcCropGs.cheyenne_intel.clm-default
What does Gswp do? Should I keep it?

A diff between lnd_in files from the NWP and the rm_indiv... tests tells me that we'd be testing everything from the latter in the new test.

Oh, and I might as well fold this into this PR, as well.

[billsacks]

Thanks, your proposal for combining the tests mostly seems mostly good. I like the idea of introducing a test of NWP with BgcCrop. However:

• It seems like your proposed test won't test the toosmall options; do you feel that's okay?

• The one issue I see with your suggestion is that it leaves compset I2000Ctsm50NwpSpGswpGs untested. I think we like to have at least one test of each out-of-the-box compset. So as much as I like the idea of combining tests: I think for now you should actually leave the current NWP tests as they are and add a new NWP test like the one you suggest.

The Gswp says to use the GSWP3 datm forcing data. I think we decided to be explicit about the forcing data in the new NWP compsets. So yes, you should leave that in the compset alias, and make sure to use the GSWP3 datm forcing in the compset long name.

[slevis-lmwg]

Status of testing...

I used this single test to breeze through the most recent commits:
./create_test ERP_D_Ld5.f10_f10_musgs.I2000Clm50Vic.cheyenne_intel.clm-vrtlay -p P93300041 -c /glade/p/cgd/tss/ctsm_baselines/ctsm1.0.dev052 & PASS

After the last commit (f62c5d0) I ran test-suites on hobart and cheyenne and got DIFF from dev052 predominantly in the clm4.5 tests (likely all of them).

I did git checkout 938f125 and ran both test-suites again.

• Hobart PASSed!
• Cheyenne had a short list of FAILs:
  SMS_Ld1.nldas2_rnldas2_mnldas2.I2000Ctsm50NwpSpNldasGs.cheyenne_gnu.clm-default BASELINE ctsm1.0.dev052: DIFF change nlevsoi from 4 to 5 PASS
  SMS_Ld1.nldas2_rnldas2_mnldas2.I2000Ctsm50NwpSpNldasRsGs.cheyenne_gnu.clm-default RUN change nlevsoi from 4 to 5 PASS
  ERP_P72x2_Lm36.f10_f10_musgs.I2000Clm50BgcCrop.cheyenne_intel.clm-clm50cropIrrigMonth_interp RUN could not launch executable Reran and got: PASS
  ERS_Ly3.f10_f10_musgs.I1850Clm50BgcCropCmip6.cheyenne_intel.clm-basic RUN could not launch executable Reran and got: PASS
  ERS_Ly5_P72x1.f10_f10_musgs.IHistClm45BgcCrop.cheyenne_intel.clm-cropMonthOutput RUN no indication of a problem Reran and got: PASS

Also on Cheyenne I got a long list of FAILures while building (it's happened twice now). I suspect that the job that submits them is running out of time, not the individual tests. How do I modify the requested time? I received this job failure (and the previous one) by email:
Post job file processing error; job 7332387.chadmin1.ib0.cheyenne.ucar.edu on host r14i7n12

[slevis-lmwg]

I resubmitted the cheyenne test suite in git checkout 938f125 and only one test remained PENDing. I will resubmit it. If it PASSES, then we should be good!

[slevis-lmwg]

@billsacks this is ready to merge, and I'm open to further review if you prefer.

[billsacks]

Also on Cheyenne I got a long list of FAILures while building (it's happened twice now). I suspect that the job that submits them is running out of time, not the individual tests. How do I modify the requested time? I received this job failure (and the previous one) by email:
Post job file processing error; job 7332387.chadmin1.ib0.cheyenne.ucar.edu on host r14i7n12

Your explanation is likely correct. The default wallclock time for the job that builds the tests is 6 hours. This is set in python/ctsm/machine_defaults.py, line 46. Why don't you go ahead and change that to 12 hours (the queue maximum); feel free to commit that change to this branch.

[billsacks]

Why don't you go ahead and change that to 12 hours (the queue maximum); feel free to commit that change to this branch.

Actually, seeing that the changes on this branch are basically complete, and may not need additional testing (I'm doing a final review now): we can just make this change in a future tag.

[billsacks]

Your latest changes basically look great to me; I really appreciate the cleanup you did in various places. However, the rework of organic_frac_squared doesn't look quite right: see inline comment below.

[slevis-lmwg]

I have performed the following tests:

1. @billsacks recommended that I submit this test from my branch and again from my master with soil_layerstruct_predefined = '49SL_10m' to manually confirm bfb:
   ERP_D_Ld5.f10_f10_musgs.I2000Clm50Vic.cheyenne_intel.clm-vrtlay
   PASS because my master in test-id d2vv1z is bfb with my branch in test-id taa26m
2. I followed @billsacks instructions for creating a new test type that automatically confirms bfb between a case using soil_layerstruct_predefined and a case using soil_layerstruct_userdefined with dzsoi values identical to those of the former case:
   ./create_test SOILSTRUCTUD_Ld5.f10_f10_musgs.I2000Clm50BgcCropGs.cheyenne_intel.clm-soil_layerstruct PASS
3. om_frac clean-up following Bill’s recommendation unfortunately FAILed when compared to dev052:
   ERP_P36x2_D_Ld5.f10_f10_musgs.IHistClm45BgcCruGs.cheyenne_intel.clm-decStart DIFF
   ERP_D_Ld5.f10_f10_musgs.I2000Clm50Vic.cheyenne_intel.clm-vrtlay PASS
   I will revert to changes that I hope will pass in the next round of testing.

@billsacks I found it: Point 3 in the above post is where I said that your recommendation gave different answers relative to dev052. In any case, let me know what you find.

[billsacks]

@slevisconsulting I think I have fixed the above issues over in #771 . However, I'd be okay with your bringing this to master before #771 is fully reviewed; if so, I would just ask that you change the if (.not. organic_frac_squared) conditional back to if ( soil_layerstruct /= '10SL_3.5m' )then ! apply soil texture from 10 layer input dataset

[slevis-lmwg]

@billsacks I think your suggestion makes sense, hence my most recent commit(s). Do you feel this change warrants running the test suites again?

[billsacks]

@slevisconsulting - If you feel confident in these changes, I don't feel you need to rerun the full test suite, but please at least run one clm50 and one clm45 test with the final version to ensure that they pass and are bit-for-bit.

[slevis-lmwg]

@slevisconsulting - If you feel confident in these changes, I don't feel you need to rerun the full test suite, but please at least run one clm50 and one clm45 test with the final version to ensure that they pass and are bit-for-bit.

Will do.

[slevis-lmwg]

git pull ! to get latest push on hobart; then ran these tests

ERP_Ld5_P48x1.f10_f10_musgs.I2000Clm50BgcCruGs.hobart_nag.clm-reduceOutput -c ctsm1.0.dev052 PASS
ERS_D_Mmpi-serial_Ld5.1x1_brazil.I2000Clm45FatesGs.hobart_nag.clm-FatesColdDef -c ctsm1.0.dev052 PASS