Bugs in pkg/layer: Does pkg/layers work in pressure coordinates? and more.

[mjlosch]

• [ ] Comparing pkg/layers output for meridional overturning circulation computations between pressure and z-coordinates, I have the impression that pkg/layers code does not work properly in pressure coordinates. In the code, there are references to w-points above and below, directions from "warm" to "cold", which has the opposite meaning in pressure coordinates than in z-coordinates and to the surface level as k=1, e.g.

https://github.com/MITgcm/MITgcm/blob/f1a296844badbe342fd00bcb35c8f59a9a016e49/pkg/layers/layers_thermodynamics.F#L67-L87 or https://github.com/MITgcm/MITgcm/blob/f1a296844badbe342fd00bcb35c8f59a9a016e49/pkg/layers/layers_fluxcalc.F#L476-L485 or https://github.com/MITgcm/MITgcm/blob/f1a296844badbe342fd00bcb35c8f59a9a016e49/pkg/layers/layers_fluxcalc.F#L582-L586

• [ ] Further, I have the impression that the following lines of code are not doing what they are supposed to do (replacing the .EQV. operator, which is most likely my own fault):

https://github.com/MITgcm/MITgcm/blob/f1a296844badbe342fd00bcb35c8f59a9a016e49/pkg/layers/layers_fluxcalc.F#L765-L767 and should probably be

           IF ( ( xx(n).GE.xx(1).AND.x(i,j).GE.xx(km) ).OR. 
      &         ( xx(n).LT.xx(1).AND.x(i,j).LT.xx(km) ) ) THEN 

[jm-c]

@mjlosch I have 2 comments:

1. There might be some problems with layers_thermodynamics.F (e.g., issue #60) but this piece of code is not used to compute the meridional overturning in density coordinate (routine is empty when LAYERS_THERMODYNAMICS is "#undef"). So, in some sense, it's a different issue.
2. The layers_fluxcalc.F has been used to get the transport in potential density coordinate or in temperature coordinate, with different ordering direction (temp increasing from bottom to top but density decreasing from bottom to top), but some comments might not be accurate. Also, I vaguely remember (but could be wrong) that it has also been used for the atmosphere in P-coordinate. This does not mean the problem you saw is not related to p-coord but its not that the implementation is completely missing.

[jm-c]

The default value for layers_krho (formerly layers_kref), which set which ref. pressure to use to define the potential density, is set to 1, independently of the choice of vertical coordinate. So, if I want to use potential density reference to the surface $\sigma_0$, I can keep the default (=1) for Z-coord. but will need to specify the Nr value for P-coord.

@mjlosch which rerefence pressure are you using for the potential density ?

[mjlosch]

for z-coordinates I use

 layers_krho(1) = 37,

which is for the ecco-50-layer grid at rF(37) approximately -1914.2 meters 9.81.1035=1934.5 dbar`

for p-coordinates this corresponding level is

# with 50 levels rF(15) = 1856.325596093750 dbar
 layers_krho(1) = 14,

in both cases I use:

 layers_bounds(1:89,1) = 0.00,
    30.0000, 30.5000, 31.0000, 31.5000, 32.0000, 32.5000, 33.0000, 33.5000,
    33.6105, 33.7210, 33.8316, 33.9421, 34.0526, 34.1631, 34.2736, 34.3842,
    34.4947, 34.6052, 34.7157, 34.8262, 34.9367, 35.0473, 35.1578, 35.2683,
    35.3788, 35.4893, 35.5999, 35.7104, 35.8209, 35.9314, 36.0419, 36.1443,
    36.2255, 36.2915, 36.3512, 36.4083, 36.4614, 36.5097, 36.5529, 36.5892,
    36.6207, 36.6476, 36.6697, 36.6873, 36.7022, 36.7155, 36.7278, 36.7386,
    36.7487, 36.7578, 36.7657, 36.7730, 36.7794, 36.7853, 36.7908, 36.7959,
    36.8009, 36.8057, 36.8104, 36.8149, 36.8195, 36.8240, 36.8283, 36.8325,
    36.8368, 36.8409, 36.8451, 36.8492, 36.8530, 36.8568, 36.8605, 36.8644,
    36.8682, 36.8720, 36.8756, 36.8789, 36.8821, 36.8854, 36.8862, 37.2322,
    37.5781, 37.9241, 38.2701, 38.6161, 38.9621, 39.3080, 39.6540, 40.0000,