Interpolation of clay, sand and om_frac can use the wrong layer

[billsacks]

Brief summary of bug

It looks like the vertical interpolation of clay, sand and om_frac from the surface dataset to the model's soil layer structure will sometimes use a different soil layer (off by 1) from what should ideally be used.

General bug information

CTSM version you are using: ctsm1.0.dev052

Note that I have made some fixes to this interpolation in #771 , but these don't address this fundamental issue.

Does this bug cause significantly incorrect results in the model's science? Incorrect, yes, but probably not significantly incorrect.

Configurations affected: Potentially all configurations that use a different soil layer structure than 10SL_3.5m (so, by default, all CLM5 configurations). I haven't checked what is actually affected in practice.

Details of bug

This refers to the logic here:

https://github.com/ESCOMP/ctsm/blob/842185f8ce5da306c7ad1a7955cffca827239c48/src/biogeophys/SoilStateInitTimeConstMod.F90#L398-L404

I think it's incorrect to compare zisoi in the model with zisoi on the file. Instead, I think we should be comparing zsoi in the model with zisoi on the file - asking if the node center of a given layer in the model falls in between two interface depths on the file.

A simple example illustrating the problem with the current code is if the zisoi values on the file are shifted by roundoff from the zisoi values in the model; e.g.:

zisoi = 0.1, 0.2, 0.3, 0.4
zisoifl = 0.099999999, 0.199999999, 0.299999999, 0.399999999

In this case, I'm pretty sure that we would put the file's values from level 3 into the model's level 2, etc. - whereas I think we'd want to put the file's values from level 2 into the model's level 2, etc. in this case: In addition to generally feeling that we should be using the model's node center (zsoi rather than zisoi) in this comparison, it also feels wrong to have an algorithm that is so sensitive to roundoff-level differences like this.

[billsacks]

@swensosc I'm tentatively assigning this to you, at least to look at and tell me if this sounds right.

[slevis-lmwg]

@billsacks I decided to look at this, too. I think I agree with you. Assuming you're changing to:
do j = 1,nlevsoifl-1
if (zsoi(lev) >= zisoifl(j) .AND. zsoi(lev) < zisoifl(j+1)) then
then I also think that the next three lines should point to j rather than j+1.

[slevis-lmwg]

Oh, and
do j = 1,nlevsoifl-1
should change to
do j = 1,nlevsoifl

Unless I'm confused about how things line up:
I'm assuming node of layer 2 sits between interfaces 2 and 3.

[slevis-lmwg]

I see in #771 that zisoi ranges 0 to nlevsoi, so my last two comments must be wrong.

[billsacks]

@swensosc doesn't have strong feelings on this. At the next CLM software meeting where @dlawrenncar is present, let's talk about whether this is, indeed, a bug, and the priority of fixing it (and possibly the related issue that the surface dataset should really contain its own definition of the layer structure, rather than having some assumed values in the code).

[billsacks]

Feeling from yesterday's meeting: @dlawrenncar doesn't have strong feelings on what's right here; he's happy to defer to @slevisconsulting and me. Fixing this isn't very high priority, but we might go ahead and do it while it's relatively fresh in our minds.