Initialization of z-grid: non-monotonous increase of zt

[fabien-roquet]

The procedure to compute zt, dzw and dzt in the function u_centered_grid of numerics.py can lead to non-monotonous zt grid. The issue was already there in the pyOM fortran code.

The constraint dzt_i = zw_i - zw_i-1, dzw_i = zt_(i+1)-zt_i , zw_i = 0.5(zt_i+zt_(i+1)) does not test/enforce whether zt is increasing monotonously. A better method would be to enforce that zt are at mid-depth of zw points:
dzt_i = zw_i - zw_i-1, dzw_i = zt_(i+1)-zt_i , zt_i = 0.5(zw_i+zw_(i+1))

[dionhaefner]

In general we like to stay compliant with PyOM whenever feasible. Can you say more, for which values of dzt did you encounter this?

[Titouan-Moulin]

We tried to use these levels

 dzt = [525.40186,  474.85938,  425.53467 , 378.23193,  333.62646,  292.21765,
 254.344    ,220.16626 , 189.69604  ,162.82532 , 139.35278  ,119.02173,
 101.539795  ,86.59766  , 73.899414  ,63.152832  ,54.094727,  46.48535,
  40.109863 , 34.77881  , 30.331055,  26.62793   ,23.543945 , 20.984375,
  18.857422,  17.092773  ,15.62793  , 14.417969 , 13.413086  ,12.583008,
  11.894531  ,11.325195  ,10.853516  ,10.463867,  10.142578   ,4.966797]

computed using Madec & Imbard (1996) function.

[dionhaefner]

I see. I can reproduce the problem, but I'm a bit stumped regarding what to do here.

• If we introduce a new way to compute vertical levels, we break backwards-compatibility with all existing setups and PyOM.
• If we do nothing, we may get silently wrong results like in your case.

Maybe the simplest (but also unsatisfying) fix would be to detect this case and throw an exception. Do we have any intuition when this case arises? Does this put a bound on how fast levels are allowed to vary?

[fabien-roquet]

Could be a question for Carsten Eden?

[dionhaefner]

Perhaps. Getting things fixed upstream is good; on the other hand it doesn't change the fact that changing this logic will silently invalidate all model setups that are currently in use.

We can sure explore fixes under the hood, but will probably have them be opt-in (i.e., flip a setting to use a different method to compute z levels, or allow them to be set manually). Any preference?

[Titouan-Moulin]

I think a condition to have a monotonous zt could be on:

dyt[1:]-np.diff(yt)/2
OR
k_indices = np.arange(len(dyt)-1)
dcond = -(-1) ** k_indices * (dyt[0] + 2 * np.cumsum((-1) ** (k_indices+1) * dyt[:-1]))

(we have dcond=dyt[1:]-np.diff(yt)/2)
dcond should be constant in sign +or - (depending if zt are positive or negative)
But this condition can still lead to "not smooth" zt, so a better condition could be dcond monotonous and constant in sign

[dionhaefner]

If you want to submit a PR that implements a better scheme to generate z levels I'm happy to review it (and provide tooling to make this conditional / hidden behind a setting).