Allow nlevgrnd < nlevurb (mainly from Keith Oleson)

[billsacks]

Description of changes

Changes needed to allow nlevgrnd < nlevurb. Most changes are from Keith Oleson, but I did a bit of work on top of his and now am opening this PR.

This branch is about a year out of date (my fault for leaving this partially completed for so long). So the next step is to merge the latest master into this and resolve conflicts.

Specific notes

Contributors other than yourself, if any: @olyson

CTSM Issues Fixed (include github issue #):
Resolves #674

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

Any User Interface Changes (namelist or namelist defaults changes)? No

Testing performed, if any:

1. Ran a short test case that triggers init_interp (mimicking @olyson 's original testing) with code mods to have nlevgrnd = 4; verified that this ran to completion:

./create_newcase --case $scratch/I2000Ctsm50NwpSpGswpGs_nlevurbgrnd_1022a --res f09_g17 --compset I2000Ctsm50NwpSpGswpGs --run-unsupported

./xmlchange NTASKS=144
./xmlchange DEBUG=TRUE

with these diffs:

diff --git a/src/main/clm_varpar.F90 b/src/main/clm_varpar.F90
index 682491039..a314c8e66 100644
--- a/src/main/clm_varpar.F90
+++ b/src/main/clm_varpar.F90
@@ -157,8 +157,12 @@ contains
       nlevsoi     =  20 
       nlevgrnd    =  nlevsoi+5
     else if ( soil_layerstruct == '5SL_3m' ) then
-       nlevsoi     =  5
-       nlevgrnd    =  5
+!KO       nlevsoi     =  5
+!KO       nlevgrnd    =  5
+!KO
+       nlevsoi     =  4
+       nlevgrnd    =  4
+!KO
     else
        write(iulog,*) subname//' ERROR: Unrecognized soil layer structure: ', trim(soil_layerstruct)
        call shr_sys_abort(subname//' ERROR: Unrecognized soil layer structure')
diff --git a/src/main/initVerticalMod.F90 b/src/main/initVerticalMod.F90
index 8ae1df226..50056c675 100644
--- a/src/main/initVerticalMod.F90
+++ b/src/main/initVerticalMod.F90
@@ -285,7 +285,7 @@ contains
        dzsoi(2)= 0.3_r8
        dzsoi(3)= 0.6_r8
        dzsoi(4)= 1.0_r8
-       dzsoi(5)= 1.0_r8
+!KO       dzsoi(5)= 1.0_r8
        
        zisoi(0) = 0._r8
        do j = 1,nlevgrnd

2. LII_D_Ld3_P144x1.f09_g17.I2000Ctsm50NwpSpGswpGs.cheyenne_intel.clm-default with the above code diffs (so nlevgrnd = 4, nlevurb = 5) and the following in user_nl_clm:

finidat = '/glade/scratch/sacks/I2000Ctsm50NwpSpGswpGs_nlevurbgrnd_1022a/run/finidat_interp_dest_saved.nc'

where that is pointing to a copy of the finidat_interp_dest.nc file created by running test case (1).

[billsacks]

@olyson - Please start from the version here (branch nlevurbgrnd on my fork): this adds a few commits on top of your work, mainly to get init_interp working.

As we discussed, the next step is probably to update this to the latest master - unless you want to do any more testing first, now that init_interp works.

I did a test merge of master (ctsm5.1.dev010) into this branch. One of the files with conflicts is initInterpMultilevelContainer.F90, which I did a lot of work on. So I made an attempt at resolving the conflicts in that file; I have attached the result. I haven't done any testing to make sure this runs or even compiles, but I have double-checked the result and it looks at least close to right. So: When you do the conflict resolution from merging ctsm5.1.dev010 into your branch, you can copy this file into src/init_interp, overwriting the version in there and so forcing the conflict resolution to look like my version.

initInterpMultilevelContainer.F90.gz

I also noticed that there are changes in this branch in CanopyTemperatureMod, which no longer exists. The relevant code in that module has been moved into BiogeophysPreFluxCalcsMod. I checked, and the version there is exactly what was in the code in CanopyTemperature prior to your changes – so I think we can simply bring your changes from CanopyTemperature into that module as is – i.e., copy the relevant block of code from your version of CanopyTemperature and replace the corresponding block of code in BiogeophysPreFluxCalcsMod with that.

I haven't looked at other merge conflicts.

[billsacks]

• A todo for myself: adding a test to the test suite that exercises nlevgrnd < nlevurb, maybe adding this to an existing test.

[billsacks]

@olyson as I mentioned by email a few days ago, I had an idea for some testing that I feel could be good to do. This was motivated by looking at the differences in SoilTemperatureMod, and not feeling sure whether the solution of the matrix will do the right thing when nlevgrnd < nlevurb (so I think there are some rows or columns of the matrix that end up not getting filled for soil columns; I wasn't sure if this would end up violating some assumption in the matrix solution). If you're fairly confident that this code (and others) is working correctly in the nlevgrnd < nlevurb case, then we could skip this testing, but if you are as uncertain as I am, then read on.

I like to come up with invariants that we can test for: start with a scenario that we trust, then make some change that exercises the new code but should give bit-for-bit results with the trusted scenario if things are working correctly. In this case, the test I came up with is:

• Set up a global, low resolution (e.g., f10) case with no urban anywhere (by tweaking the PCT_URBAN surface dataset field, and adding that area to PCT_NATVEG; or maybe easier, by setting toosmall_urban to 100, which should hopefully result in 0 urban anywhere). Other than urban being 0% everywhere, the control case should be out-of-the-box. Run for (say) 1 year. The next step may be easiest if this case has relatively small nlevgrnd (e.g., using the 5SL_3m soil layer structure... but other than that, it may be good to use a standard Clm50 model setup rather than Nwp, so that we have full subgrid heterogeneity).
  • Update: I realized that, ideally, we'd include pervious and impervious roads, since they use nlevgrnd rather than nlevurb. So maybe keep urban at its normal cover, but change the urban parameters so urban areas are entirely roads, if that's possible?
• Then run a similar case, but with nlevurb set to nlevgrnd + 2 (for whatever nlevgrnd was in the original case). I'm not sure what this will take, but hopefully, given that there is no urban anywhere, it should be fairly easy to just put in some dummy values where they're needed. If I'm thinking things through correctly, this should be bit-for-bit with the control, since the addition of some urban levels shouldn't have any impact on the solution if urban is 0% area everywhere. This will ensure that simply expanding the number of levels in various arrays doesn't change answers for non-urban landunits.

If it's not too hard, I think it would be good to run this with both Clm45 and Clm50 because I think they use different solution methods in some places (no need for Clm51).

Finally, if it's not too hard, it would be good to confirm that urban areas are also working correctly when nlevgrnd < nlevurb... though this seems less likely to be a problem, so seems lower priority. I think this could be done similarly to the above, but setting urban to 100% everywhere; if it's hard to do that globally, then this one could be a single-point all-urban case (in which case, probably run for longer, like 20 years). Other than setting urban to 100% everywhere, the control case would have out-of-the-box settings. Then run a similar case, but change nlevgrnd to nlevurb - 2 (by changing the soil layer structure); again, I think this should be bit-for-bit.

• Update: Actually, I just realized that roads will cause a problem with this test. If it's easy, this could be handled by changing urban parameters so there are no roads anywhere. But if this is a pain, then this test is probably less important than the first.

• It would be good to do some manual invariant testing like this, unless @olyson feels confident enough in the operation of the code that he feels it isn't necessary.

[billsacks]

@olyson - Okay, I have merged ALMOST all of the files with the latest master; these were mostly pretty easy. Unfortunately, the only one that I haven't merged is going to be by far the hardest to merge, and I'm hoping you'll still be willing to help with this (sorry!): SoilTemperatureMod.F90.

The current status is that I have basically reverted SoilTemperatureMod.F90 to how it was on master; I have just applied the first few diffs from this nlevurbgrnd branch so that the code compiles. There were 101 merge conflicts in this file, so trying to do an automated merge felt kind of pointless. Instead, I think it will be easiest to look at the changes you made originally, and reapply those to the new version, either directly or conceptually. Here's a link to see your original diffs: https://github.com/ESCOMP/CTSM/pull/1195/files/45a668077dcf078be24d31d5ed3f2116db63070f#diff-a3ab86cf37e1bbb63513066dbe487e1d355b4ed31ce3623cabde340895dcead6. That said, it's possible that many of the early diffs in this file will merge cleanly, so if you'd like, we can try to let git do some of the work... let me know if you'd like my help with this. (I would do it by checking out the version of this branch just before I did the merge, then doing 'git merge ctsm5.1.dev010' and seeing what the resulting SoilTemperatureMod.F90 looks like.)

From looking through the diffs, I noticed that many (but not all) changes involve places where nlevgrnd is referenced. I think the main change in SoilTemperatureMod was that Negin (working with Sean) consolidated many of the little subroutines into a much smaller number of generic subroutines. However, there may also be new references to nlevgrnd in the latest version of SoilTemperatureMod, so it could be good (after applying your mods) to search for that variable and ensure that all remaining uses of nlevgrnd are truly supposed to be nlevgrnd.

[olyson]

I think it will be easiest to look at the changes I made originally, and reapply those to the new version manually. As you saw, there were a bunch of little subroutines where I had to make changes. Presumably, with the refactor, this will be not quite as time-consuming. I assume I can start by updating my version of your branch with your latest commit to that branch.

[billsacks]

I think it will be easiest to look at the changes I made originally, and reapply those to the new version manually. As you saw, there were a bunch of little subroutines where I had to make changes. Presumably, with the refactor, this will be not quite as time-consuming. I assume I can start by updating my version of your branch with your latest commit to that branch.

Sounds good, Keith. Yes, that update is the best first thing to do.

[billsacks]

@olyson - I just updated the branch with some infrastructure changes (a few hours of a peeling onion effect while I was restoring your check of the levmaxurbgrnd dimension on the restart file). These shouldn't overlap with anything you're doing with SoilTemperatureMod. So if you haven't started that work yet, you might as well grab this update now, but if you're in the middle of it, I think you can safely finish that before updating to grab my latest commits.

I have run these two tests and confirmed they pass and are bit-for-bit with ctsm5.1.dev010:

ERP_D_P36x2_Ld3.f10_f10_musgs.I1850Clm50BgcCrop.cheyenne_intel.clm-default
ERP_P36x2_D_Ld5.f10_f10_musgs.I2000Ctsm50NwpBgcCropGswpGs.cheyenne_intel.clm-default

[olyson]

Thanks. I'm deep into SoilTemperatureMod now, so I'll wait to update.

[olyson]

I've restored the changes I made to SoilTemperatureMod.F90. I've also updated some new code in SnowHydrologyMod.F90 and SurfaceWaterMod.F90 that was causing failures in debug mode.
I've updated to your last commit and push that to the branch. Although I'm not sure why that update shows as a separate commit with changed files, so perhaps I've done that wrong.
I reran the two tests you ran and they pass.
I also ran a couple of my own tests and they were bfb with ctsm5.1.dev010.
I ran a test with nlevgrnd=nlevsoi=4 (and nlevurb=5) and that ran to completion.
So I think we could proceed with the extra testing you proposed. However, I don't think the code will handle urban areas as roads only. The code in, for example, UrbanFluxes, expects and urban canyon to solve for air temperature.

[billsacks]

Awesome, @olyson - thanks so much! I checked the merge, and it looks like you did that right.

As for the testing I suggested: First off, if you can think of alternative testing that would accomplish something similar, I'm not tied to the testing I suggested: it's just what came to mind initially. If you go with my suggestion, then it seems okay to me to leave out roads if that makes things easier - i.e., ignoring this:

Update: I realized that, ideally, we'd include pervious and impervious roads, since they use nlevgrnd rather than nlevurb. So maybe keep urban at its normal cover, but change the urban parameters so urban areas are entirely roads, if that's possible?

It means that the test wouldn't cover quite all cases, but I think it would still be useful.

But since my testing idea was mostly motivated by the extensive changes in SoilTemperatureMod, and because you just revisited those changes, I am happy to defer to you in determining what is and isn't necessary to test, based on what you do / don't feel confident about. That said, I'm happy to chat with you if some joint brainstorming would be helpful.

[olyson]

I've done testing with 0% urban per your first suggestion (details below). The CLM4.5 and CLM5 tests were successful. There were differences associated with the 100% urban tested as expected because of roads.

Global cases with 0% urban

1. ./create_newcase --case clm50_nlevurbgrnd_f10_GSWP3V1_nourb_2000 --compset I2000Clm50BgcCropGs --res f10_f10_musgs --project P93300641 --run-unsupported
   Set the following in user_nl_clm:
   toosmall_urban = 100
   soil_layerstruct_predefined = '4SL_2m'
   Then run for 1 year – SUCCESSFULL
2. ./create_clone --case clm50_nlevurbgrnd_f10_GSWP3V1_nourb_nlevurb10_2000 --clone clm50_nlevurbgrnd_f10_GSWP3V1_nourb_2000
   Change nlevurb=10 in clm_varpar.F90
   Create a new surface dataset using the new raw urban dataset (which has nlevurb=10). Only differences between this dataset and the default dataset are urban parameters. Specify this new dataset in user_nl_clm:
   /glade/work/oleson/ctsm_sacks_nlevurbgrnd/tools/mksurfdata_map/surfdata_10x15_hist_78pfts_CMIP6_simyr2000_c201028.nc
   Then run for 1 year – SUCCESSFULL and BFB with first simulation.
3. Repeat for CLM4.5 – SUCCESSFULL and BFB
   Single point cases with 100% urban
4. As you said, roads will cause a problem with this proposed test (non-BFB). You can set to pervious road to zero but the impervious road will still exist and will be affected by a change in nlevgrnd. I ran this test anyway and got differences as expected. However, they were fairly small, as one might expect because road is but one of five “surfaces” in the urban canyon.
5. A similar test I’ve done is a 100% urban single-point simulation with nlevgrnd < nlevurb. This compiled and ran fine in DEBUG mode.

[billsacks]

Awesome, @olyson - thanks so much for doing that additional testing! I think I can take this branch from here in the next few days – adding a test and doing a bit of minor cleanup that I noticed when I reviewed it a few days ago.