Fix interp. of surfdat soil layers so we can use interp. for 10SL case

[billsacks]

Previously, on master, there was special-purpose code for the 10SL case
that avoided doing interpolation from the surface dataset to the soil
variables in the model. (In a recent commit, @slevisconsulting changed
the logic to be based on organic_frac_squared, but that dependence was
not correct.) It would be cleaner - especially now that we allow
user-defined soil layer structures - if we could use the general
interpolation code always, rather than sometimes having special-purpose
code that avoids doing the interpolation.

This commit accomplishes that generality. In order to preserve answers
for clm45 cases, I needed to make three changes:

1. The constant 0.025 needed an _r8 appended to it; otherwise, zisoifl
   could differ by roundoff from zisoi. I'm confident that this is a
   correct fix.

2. I changed which level is used when zisoi(lev) is exactly equal to
   zisoifl(j) for some j: I changed the conditional in the following:

   if (zisoi(lev) >= zisoifl(j) .AND. zisoi(lev) < zisoifl(j+1)) then
   clay = clay3d(g,j+1)
   sand = sand3d(g,j+1)
   om_frac = organic3d(g,j+1)/organic_max
   endif

   to:

   if (zisoi(lev) > zisoifl(j) .AND. zisoi(lev) <= zisoifl(j+1)) then

   I am pretty sure, but not positive, that this is a correct fix in
   general: Previously, when the zisoi values in the model exactly lined
   up with the zisoi values in the file, we would set clay in model
   level j equal to the value from level j+1 from the surface dataset
   (and similarly for sand and om_frac); in the new code, we use level j
   from the surface dataset for model level j in this case. I assume
   this is the right thing to do, though am not 100% positive.

3. I changed the way zisoifl is calculated for the lowest layer, so that
   it matches zisoi(nlevsoi) when running with 10SL_3.5m. I think this
   is the correct thing to do, but I'm not positive of this. Previously,
   we had:

   zisoi(10) = 0.38018819123227207690E+01
   zisoifl(10) = 0.34330930154359391437E+01

   This commit changes zisoifl(10) to match zisoi(10).

With these changes, the following tests are bit-for-bit with master:

SMS_D_Ld1_P4x1.f10_f10_musgs.I2000Clm45BgcCropQianRsGs.bishorn_gnu.clm-default
(with additional changes to define this compset properly)

SMS_D_Ld1_P4x1.f10_f10_musgs.I2000Clm50BgcCropQianRsGs.bishorn_gnu.clm-default

[billsacks]

Diff (1) from my main PR comment: addition of _r8

[slevis-lmwg]

Why do we calculate these depths if they are on the file and could be read? Seems risky to be calculating depths that are assumed to be on a file, when the file may change and we could be calculating the wrong depths at that point.

[billsacks]

Yup, I completely agree that these should be on the file. But they aren't on there right now. Long-term, the best solution would absolutely be to put these levels on the surface datasets, but I felt that was too much to fold in with #759 or this PR.

[billsacks]

Diff (3) from my main PR comment (and related changes around here): changing definition of zisoifl for lowest level

[slevis-lmwg]

Considering zsoifl represents the depths in the soil texture file, we introduce an abstraction here by defining a zsoifl for a soil layer that doesn't really exist, i.e. nlevsoifl + 1. Do you feel comfortable with that?

[billsacks]

I agree that this feels weird, and certainly I'd feel less comfortable with this if zsoifl were used for anything more than defining zisoifl. But since zsoifl is just used to define zisoifl, I figured that the best thing to do here is to define zsoifl so that the zisoifl values end up as desired - which involves setting an assumed value for zsoifl(nlevsoifl+1). Unless you see a cleaner way to accomplish this goal?

[slevis-lmwg]

Sounds good. To help with the long term development and institutional memory of what is going on here, how about a comment that explains how this code might change if we were reading the depths from the file (a file); and which specific file we are referring to; and the danger of switching to a new file given the current implementation. If such a comment exists somewhere else, I suggest copying it here, too.

[billsacks]

Comment added in 40a5881

[billsacks]

Diff (2) from my main PR comment: changing the logic for when zisoi == zisoifl

[slevis-lmwg]

What happens here for, let's say, lev = 20 if zisoi(20) happens to be greater than zisoifl(10). I think that clay, sand, and om_frac do not get valid values for lev = 20 in that case. Even if somehow they do get valid values, maybe we need an explicit statement, so that there is no doubt about their values. I'm thinking something like this:
else if (zisoi(lev) > zisoifl(j+1)) then
clay = clay3d(g,j+1)
sand = sand3d(g,j+1)
om_frac = organic3d(g,j+1)/organic_max

[slevis-lmwg]

Oh, and what if zisoi(lev) happens to be less than or equal to zisoifl(j). That case is not accounted for, either, I think.

[billsacks]

That's a good point. I believe these problems existed before this PR; do you agree? (I'm just trying to distinguish between possible new problems in this PR and pre-existing problems.)

For this

Oh, and what if zisoi(lev) happens to be less than or equal to zisoifl(j). That case is not accounted for, either, I think.

I'm thinking you mean: zisoi(lev) is <= zisoifl(1), for lev = 2 or more.

Now, that said: I believe the code is actually working in a reasonable way right now, it just isn't as explicit as it could be: if these variables (clay, sand, om_frac) aren't explicitly set, they'll remain at their values from the last iteration of the level loop. So, for the case where zisoi(lev) is <= zisoifl(1), I believe we'll end up using values from level 1 of the file (because these will have been set in the lev .eq. 1 block); this seems reasonable. And, for the first case you raise, I believe we'll end up using values equal to the those in the model layer above this one; this may not be the right thing to do (e.g., if the last model layer took values from file layer 8: I think we'd really want to use file layer 10 in this case), but it doesn't seem terrible.

We could fix this to be more explicit; if so, I think this code would become cleaner - and unit testable! - if we extracted this logic for finding the right layer into its own routine, tied in with the parenthetical under "Caveat c" in this comment #759 (comment) .

However, it might make the most sense to wait and do this as part of resolving #772 , once we agree on how we want to solve that one - rather than fixing this implementation just to throw it out and replace it with something more correct.

What do you think?

[slevis-lmwg]

I'm fine with that. Again though, how about adding a comment in the code with your best understanding of how this works now versus how it should work ideally.

[slevis-lmwg]

moving the logic matching soil layers in the model with soil layers on the file to outside the column loop, since this matching is independent of column

Regarding the parenthetical under caveat c in #759, if there's a chance that some day we will end up with varying soil layer depths, it may not be worth changing. I'm just raising as a potential concern; @swensosc and others may have better insight.

[billsacks]

Comment added in 40a5881

[billsacks]

@slevisconsulting and @swensosc - the changes here appear to resolve the issues we were discussing in #759 . However, I'd like a review from both of you to sign off on the changes I made. As noted above, I have verified that this is bfb for two cases, but it's possible that this changes answers in some other cases that I didn't test. I expect answer changes in some situations because I have made some changes to the interpolation from per-layer variables on the surface dataset to those in memory; I'm just not sure if those changes will actually impact any of our current soil layer structures (I think there's a good chance they won't, but they might).

In making these changes, I came to feel that the current interpolation has some more fundamental problems, but I'll open a separate issue for that.

[billsacks]

Also, @slevisconsulting - let me point out what I did here in terms of git/GitHub: In order to open this PR in ESCOMP, I pushed your branch up to ESCOMP. When this PR is merged, it will get merged into your branch on ESCOMP (not on your fork). You can then update your local branch to match this, then push your local branch back up to your fork (or I can do this).

[slevis-lmwg]

@billsacks I'm posting a few things to think about. I'm interested in your responses.

[slevis-lmwg]

Why do we calculate these depths if they are on the file and could be read? Seems risky to be calculating depths that are assumed to be on a file, when the file may change and we could be calculating the wrong depths at that point.

[slevis-lmwg]

Considering zsoifl represents the depths in the soil texture file, we introduce an abstraction here by defining a zsoifl for a soil layer that doesn't really exist, i.e. nlevsoifl + 1. Do you feel comfortable with that?

[slevis-lmwg]

What happens here for, let's say, lev = 20 if zisoi(20) happens to be greater than zisoifl(10). I think that clay, sand, and om_frac do not get valid values for lev = 20 in that case. Even if somehow they do get valid values, maybe we need an explicit statement, so that there is no doubt about their values. I'm thinking something like this:
else if (zisoi(lev) > zisoifl(j+1)) then
clay = clay3d(g,j+1)
sand = sand3d(g,j+1)
om_frac = organic3d(g,j+1)/organic_max

[slevis-lmwg]

Oh, and what if zisoi(lev) happens to be less than or equal to zisoifl(j). That case is not accounted for, either, I think.

[billsacks]

@swensosc doesn't really have any feelings on this, and feels that whatever we came up with is probably better than what was there before. So I'm removing him as a reviewer.