Spurious values of temperature in sigma layers

[willaguiar]

I was looking at the temperature (thetao) values of mom6-panan, online binned into sigma coords (online monthly means). It seems like the binning of the temperature is a little wacky ( i.e., it values like -300oC or + 5000oC in some places, in all three resolutions, Check below).
I thought of mentioning that here, so people are aware in case they use thetao(sigma).

Tagging diagnostics list issue #54

[adele-morrison]

Have you checked the thickness at these locations? They might be >0 m but effectively 0 m. The minimum layer thickness is 1.0E-10 m.

[willaguiar]

ok, that makes sense. Taking an average dxt ~ 3000m, and dyt ~6000m for panan 1/20th that means the minimum volume value should be about ~0.0018 m^3? I tried filtering cell volumes <0.15 m^3, but the spurious temperatures still persists.
I can only get rid of them if I filter volumes <5e+9. Even filtering 1e+9 is not enough

[adele-morrison]

Do you have a map to see where they are? And is there a thickness variable?

[willaguiar]

I had a look on the locations of these values for the year 2000 (direct model monthly mean output) , they seem to be mostly along continental slope around Ross Sea:

I located the negative spurious values by filtering out any T> -10C (the minimum temperature I found in z* layers was ~-1.9C). When I look at the thickness of the sigma cells in these spurious values they seem to be quite big [d] :

[angus-g]

Just a quick update: I re-ran the first month of 2000 and saw spurious temperature values. I then ran the first hour, outputting each timestep. I can still see interesting temperatures in one of the timesteps here:

>>> d.thetao_rho.isel(time=3, xh=2185, yh=246).values[65:76]
array([-1.6350667e+00, -1.6472415e+00, -1.6616650e+00, -1.6761010e+00,
       -1.6831371e+00,  4.6796357e+03,  4.6796357e+03,  4.6796357e+03,
       -1.6841160e+00, -1.7331796e+00,  5.5577391e-01], dtype=float32)

I saved the corresponding thicknesses in this timestep:

>>> d.h_rho.isel(time=3, xh=2185, yh=246).values[65:76]
array([ 2.4866812e-01,  3.8814417e-01,  4.2784742e-01,  4.8299605e-01,
       -1.0547950e+03,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
        1.0758734e+03,  0.0000000e+00,  1.0592498e+03], dtype=float32)

It looks like some of the density bins are missed, and something about this is messing up the interpolation of temperature?

>>> d.rhopot2.isel(time=3, xh=2185, yh=246).values[:48]
array([1036.7155, 1036.7178, 1036.7192, 1036.7205, 1036.7214, 1036.7224,
       1036.7234, 1036.7244, 1036.7256, 1036.7329, 1036.7396, 1036.743 ,
       1036.768 , 1036.8319, 1036.9534, 1037.0736, 1037.0994, 1037.1187,
       1037.123 , 1037.1276, 1037.1301, 1037.1316, 1037.1318, 1037.1309,
       1037.1284, 1037.1244, 1037.1193, 1037.1132, 1037.107 , 1037.1007,
       1037.0933, 1037.0829, 1037.0729, 1037.0642, 1037.059 , 1037.0564,
       1037.0547, 1037.0504, 1037.0428, 1037.0363, 1037.0334, 1037.0359,
       1037.0425, 1037.0475, 1037.0516, 1037.056 , 1037.0636, 1037.067 ],
      dtype=float32)

There is a bit of an inversion in rho2 in the column, which could be related? The model-space temperature and thickness values look normal, otherwise.

[willaguiar]

Thanks for looking into that Angus!
Interesting that the spurious values are exactly the same. Density bins [65:76] are between 1037.02 and 1037.13 so right around the density inversion, which I guess it would explain the spurious temperatures if the model is interpolating to get the values in the missing density bins.

Although, I thought MOM6 binned values by averaging over the pre-defined density bins instead of interpolating them (similar to what we do offline for OM2)?

[angus-g]

Although, I thought MOM6 binned values by averaging over the pre-defined density bins instead of interpolating them (similar to what we do offline for OM2)?

No, I think it's done by polynomial interpolation. Because this is output on the rho2 diagnostic grid, it uses the coord_rho ALE code to generate the grid. For each column, it'll build a polynomial reconstruction and interpolate. I think for diagnostic coordinates, it defaults to PPM H4 interpolation (the model coordinate would default to P1M H2).

There is also this rather ominous warning in the get_polynomial_coordinate docstring 👻

!! IT IS ASSUMED THAT THE PIECEWISE POLYNOMIAL IS MONOTONICALLY INCREASING.
!! IF THIS IS NOT THE CASE, THE NEW GRID MAY BE ILL-DEFINED.

There are a couple of diagnostics which use a different method for determining density, which may indeed be binning... They are h_rho, uh_rho, etc. But there's no temperature on these "pure potential density coordinates".

[willaguiar]

OK, that makes sense.. do we have a list/way to figure out which outputs are averaged and which are interpolated? (probably any output that is not multiplied by h_rho?). Cause I imagine any rho2 output that is interpolated should have the same problem in regions/times of density reversal.

@aekiss , tagging you in case this matters for OM3 too

[angus-g]

do we have a list/way to figure out which outputs are averaged and which are interpolated

Everything on the ocean_model_rho2 component in the diag_table will be using the ALE method. The only outliers are the _rho-suffixed variables like h_rho, uh_rho, vh_rho (and some GM ones) that are on ocean_model.