Bug in finding layers with given depths in MPAS-Ocean high-frequency output

[xylar]

I believe there is a bug, introduced 9 years ago in this commit:
3c77a2e#diff-904437d562948f083fb74894026c6c469cba2ff0e057a29cc1b33e80534736d9
and propagated many places from there within the high-frequency analysis member:

         ! for now, just get close enough
         iLevelTarget = 1
         do iLevel=2,nVertLevels
           if(refBottomDepth(iLevel) > 100.0_RKIND) then
              iLevelTarget = iLevel-1
              exit
           endif
         enddo

As far as I can tell, the level containing a depth of 100 m has index iLevel not iLevel - 1.

As an example, here are the values of refBottomDepth for the 80-layer grid:

 refBottomDepth = 1.98678646935581, 4.15411939554175, 6.51847301561518, 
    9.09782457480988, 11.9117370565737, 14.9816089343096, 18.3308398810405, 
    21.9848287330396, 25.9713890184824, 30.3208291170111, 35.0661986371497, 
    40.2435334164226, 45.8920987538932, 52.0547139142045, 58.7780736084832, 
    66.1129821226365, 74.1148329705315, 82.8439984303522, 92.3661273493374, 
    102.752599076403, 114.081046662963, 126.435778783603, 139.908344944848, 
    154.597987883586, 170.612224022118, 188.067210086335, 207.088461553922, 
    227.811080298273, 250.38029644816, 274.951640909234, 301.691130033472, 
    330.775262096312, 362.390624492984, 396.733309415583, 434.007926026355, 
    474.426080833997, 518.204206086992, 565.560782817531, 616.712940125286, 
    671.872076779591, 731.238809447229, 794.997370867892, 863.309316527158, 
    936.307184595767, 1014.08826238255, 1096.70883766392, 1184.17979797662, 
    1276.46376106484, 1373.47419118303, 1475.07667395912, 1581.092497774, 
    1691.30425184573, 1805.46295097985, 1923.29619676545, 2044.51666335532, 
    2168.83040784201, 2295.94444415594, 2425.57314654357, 2557.44346644347, 
    2691.29872299424, 2826.90114580615, 2964.03328758224, 3102.49847295099, 
    3242.1205209147, 3382.74290983668, 3524.2275800194, 3666.4535034604, 
    3809.31514413569, 3952.72088681684, 4096.5915105551, 4240.85873802501, 
    4385.46388801116, 4530.35664731102, 4675.493968649, 4820.83908311944, 
    4966.3606305003, 5112.03189866471, 5257.83015784407, 5403.73608416561, 
    5549.73326321215 ;

If you just count the entries, you will see that a depth of 100 m would be the 20th index. This is the layer where refBottomDepth = 102.752599076403. But if you follow the logic of the original code, copied many times, you incorrectly subtract one and get the 19th entry, which is one layer to shallow with refBottomDepth = 92.3661273493374. This is not the layer that contains a depth of 100 m.

[xylar]

I am making this bug report mostly because this error is now also potentially propagating to other code (see, E3SM-Ocean-Discussion#94 (comment)).

[xylar]

@mark-petersen and @vanroekel, if you agree that this is a bug, I'm happy to submit a fix.

[vanroekel]

@xylar thanks for pinging me on this. I agree. this is a bug. It should be iLevel. Quick thought though, do we maybe want to compare refBottomDepth(iLevel-1) and refBottomDepth(iLevel) and pick the one closer to 100m? Maybe not too important though.

[xylar]

I think the code should have simply been:

         ! for now, just get close enough
         do iLevel=1,nVertLevels
           if(refBottomDepth(iLevel) > 100.0_RKIND) then
              iLevelTarget = iLevel
              exit
           endif
         enddo

possibly with some handling of the possibility that iLevelTarget never gets set (e.g. if the bottom of the ocean is too shallow).

[xylar]

Quick thought though, do we maybe want to compare refBottomDepth(iLevel-1) and refBottomDepth(iLevel) and pick the one closer to 100m? Maybe not too important though.

If we want the closest layer, we should be looking at refZMid, rather than refBottomDepth. We don't want the layer whose bottom is closest to the given depth, we want the layer whose middle is closest. (We also need to keep in mind that I think the sign of refZMid is the opposite of refBottomDepth).

[vanroekel]

great point @xylar -- I'm happy with your fix you propose just above.