Allow users to specify snow layer structure more flexibly

[billsacks]

The current snow layer structure is hard-coded. We should probably let users specify this more flexibly, as in #728 .

Details of the "user interface" still need to be worked out. For example, there are currently separate min and max depths for each layer; should we allow/force the user to specify these independently, or should we have the user specify a single layer structure from which we derive both min and max depths of each layer?

In addition to the number of layers and layer depths, I have noticed that there is a hard-coded assumption that appears in a few places throughout the code that the minimum depth for having a resolved snow layer is 0.01 m; this should probably be connected to the snow layer structure itself.

I'm not sure what to do about h2osno_max. Probably this needs to be settable independently from the snow layer structure, but how do we ensure that h2osno_max is "reasonable" given the snow layer structure?

[billsacks]

@swensosc here is the issue we were talking about.

[billsacks]

From today's ctsm-software meeting: We'll have the users specify the minimum amount of information needed to generate the full vectors, dzmin, dzmax_l and dzmax_u. Note that these can be generated via recursive formulas, as described in the attached document. (This document was written by @whlipscomb Aug, 2013, and was later implemented by @lvankampenhout .)

A firn model for CLM.pdf

The key equations in that document for these purposes are (1.3), (1.4) and (1.5):

dzmin(n) = dzmax_u(n-1)/2
dzmax_u(n) = 2*dzmax_u(n-1) + 0.01
dzmax_l(n) = dzmax_u(n) + dzmax_l(n-1)

The document states that these equations can be used to generate values for layers 3 and higher. To me it looks like they mostly work for layer 2 as well, but that dzmin(2) uses a different formula, probably to ensure that dzmin(2) > dzmin(1). It may be that the problem of dzmin(n) being <= dzmin(n-1) could exist for higher layers, too, depending on the initial values of dzmin and dzmax_u. So I'd suggest that we have the user specify just dzmin(1), dzmax_u(1) and dzmax_l(1), then generate the rest, but that we come up with some general way of ensuring that dzmin(n) > dzmin(n-1) for all n.

We'll need to ensure that the user-specified values satisfy: dzmin(1) < dzmax_u(1) < dzmax_l(1). I think it's important that this relationship holds for all n (not just 1), but I believe that's guaranteed given the above recursive formulas, as long as it holds for n=1. It also looks like these recursive formulas ensure that dzmax_u and dzmax_l grow monotonically, but as noted above, we'll need to do something more to ensure that dzmin always grows monotonically.

Examples of how we could deal with the dzmin issue are:

if (dzmin(n) <= dzmin(n-1)) then
   dzmin(n) = (dzmin(n-1) + dzmax_u(n-1))/2
end if

or

if (dzmin(n) <= dzmin(n-1)) then
   dzmin(n) = dzmax_l(n-1)/2
end if

though the latter could violate some other constraints sometimes - possibly leading to dzmin(n) > dzmax_u(n).

We may need to give some thought to other constraints: I could imagine there being pathological situations where chosen values of dzmin, dzmax_l and dzmax_u (either in layer 1 or the generated values for n > 1) lead to problems during snow layer combination or division. @slevisconsulting would you be able to think through what additional constraints there should be, if any?

[billsacks]

One other note from today's meeting: Once we generalize this, we'll do away with the current hard-coded vector, allowing even the standard, out-of-the-box cases to be generated with these equations. I expect roundoff-level changes from this. We should print out the old and new values of these three vectors to ensure that differences are no greater than roundoff-level.

[billsacks]

@slevisconsulting see the changes I just made in ctsm1.0.dev059, which consolidate the lake and non-lake references to 0.01 as a threshold for snowpack initiation, in SnowHydrologyMod.F90: subroutine BuildFilter_SnowpackInitialized. I'm not sure if that's the only place where this hard-coded value appears.

[lvankampenhout]

Just as a comment, not sure if this is relevant for this issue, but many of the more detailed snow/firn models (IMAU-FDM, Crocus, SNOWPACK) implement snow layers in a Lagrangian fashion, i.e. snow layers are not statically layered but are ever changing in vertical order and layers can "advect downwards" as it were. What I mean to say with this is, if a more realistic implementation of snow stratigraphy is the overarching goal here, then just redefining layer thicknesses but keeping the Eulerian structure is not necessarily going to improve matters much, and a Lagrangian scheme would make perhaps more sense. However, this is not compatible with the current implementation of snow layers (as being an extension of soil layers) and the rework will be more substantial.

It is not clear to me from the first few posts, however, what the rationale is. (If it is about removing hardcoded numbers, then it makes sense).

[billsacks]

@lvankampenhout thanks for your thoughts - this is helpful to know. For now the goal is much more limited: simply removing some hard-coded numbers and letting the user specify these. But this is good to keep in mind for the future.

[slevis-lmwg]

Results from three tests with different nlevsno values:

nlevsno = 12: ./create_test SMS_P720x1_Ln3.hcru_hcru.I2000Clm50BgcCruGs.cheyenne_intel.clm-coldStart -c /glade/p/cgd/tss/ctsm_baselines/ctsm1.0.dev059 PASS because dzmax_l for layers > 1 not used due to insufficient snow

nlevsno = 5: ./create_test SMS_Ld1.nldas2_rnldas2_mnldas2.I2000Ctsm50NwpSpNldasGs.cheyenne_gnu.clm-default -c /glade/p/cgd/tss/ctsm_baselines/ctsm1.0.dev059 PASS but not bfb because
dzmax_l(2) = 8.000000000000000E-002 instead of 7e-2
dzmax_l(3) = 0.190000000000000 instead of 0.18
dzmax_l(4) = 0.420000000000000 instead of 0.41
dzmax_l(5) = 0.890000000000000 instead of 0.88

nlevsno = 3: ./create_test SMS_Lm1_D.f10_f10_musgs.I2000Clm50BgcCrop.cheyenne_intel.clm-snowlayers_3_monthly -c /glade/p/cgd/tss/ctsm_baselines/ctsm1.0.dev059 PASS but not bfb because
dzmax_l(2) = 8.000000000000000E-002 instead of 7e-2
dzmax_l(3) = 0.190000000000000 instead of 0.18

The equations in the "firn model" pdf (linked to an earlier post) were written for layers 3 and higher. Changing the equation for dzmax_l in the pdf to
dzmax_l(j) = dzmax_u(j) + dzmax_l(j-1) - 0.01_r8 for j == 2
allows the third test above to PASS!

nlevsno = 5 test resubmitted PASS

Longer nlevsno = 12 test
SMS_D_Ly6_Mmpi-serial.1x1_smallvilleIA.IHistClm45BgcCropQianGs.cheyenne_intel.clm-cropMonthOutput PASS

[billsacks]

The equations in the "firn model" pdf (linked to an earlier post) were written for layers 3 and higher. Changing the equation for dzmax_l in the pdf to
dzmax_l(j) = dzmax_u(j) + dzmax_l(j-1) - 0.01_r8 for j == 2
allows the third test above to PASS!

Ah, yes, good catch. Sorry that I missed that before!

In that case, since there are special cases for two out of the three values in level 2, I'd suggest that we have the user explicitly set all three values in levels 1 and 2, then use the level 2 values to derive levels 3+ (along with error checking to ensure that the constraints I mentioned in #729 (comment) are satisfied). I haven't looked at your PR yet, so maybe you already do something like this... I'll look at that next.