Motivating failure (tmp_local_artifacts/global_run/FINDINGS.md,
2026-09-11): a global tripolar aquaplanet NaN’d at outer step 7. The
ONLY diagnostic on offer was a bare non-finite-face count — “producer
0/0 upstream, investigate”. A human needed several runs + a
bisection to find the actual cause: nu_h = 2.0e4 with dt = 900 s
at the ~3 km polar cells gives a viscous-diffusion number of 1.65
against an explicit-Laplacian bound of 0.125 — 13x over — and
bound_kh (the per-cell runtime clamp that would have protected
against exactly this) was never enabled.
Nothing at configure time said any of that. This module is the fix:
a small set of checks run once, AFTER the real per-cell metric
arrays exist (ocean_metrics_t, filled by configure_ocean_metrics
+ land-masked by configure_ocean_land_mask) — never off the
nominal &grid_nml dx/dy, which are DEGREES on spherical/tripolar
grids and are in any case the NOMINAL spacing, not the smallest
actual cell (a spherical/tripolar grid’s smallest cell can be an
order of magnitude below nominal near a pole).
Each check reports the computed number, the limit, the knob(s)
responsible, and a concrete fix — never a bare “X exceeded”.
Severity: a VIOLATED HARD STABILITY BOUND (viscous CFL, the tracer
diffusive number) is an ERROR — returned via the P0 ierr status
(OCEAN_STATUS_ERR_SETUP), never error stop (this module always
has ierr to report through; configure_ocean_metrics et al. do
the same). A MARGINAL/QUALITY issue (the Munk-layer resolution
criterion, ah_max silently clamping nu_h) is a WARNING — logged,
run proceeds. The viscous-CFL check is itself downgraded from ERROR
to an informational WARNING when the run already carries automatic
runtime protection (bound_kh, or the stress_tensor operator’s
own always-on per-cell CFL limiter) — the raw nu_h*dt/dx^2 number
is no longer what the kernel actually uses in that case, so a hard
configure-time failure would be a false positive (see
docs/CLOSURE_MATRIX.md / rdb_ocean_horizontal_viscosity.F90
module header for bound_kh / stress_tensor semantics).
Nodes of different colours represent the following:
Solid arrows point from a submodule to the (sub)module which it is
descended from. Dashed arrows point from a module or program unit to
modules which it uses.
Where possible, edges connecting nodes are
given different colours to make them easier to distinguish in
large graphs.
Nodes of different colours represent the following:
Solid arrows point from a submodule to the (sub)module which it is
descended from. Dashed arrows point from a module or program unit to
modules which it uses.
Where possible, edges connecting nodes are
given different colours to make them easier to distinguish in
large graphs.
Variables
Type
Visibility
Attributes
Name
Initial
real(kind=wp),
private,
parameter
::
DIFFUSIVE_NUMBER_LIMIT
=
0.5_wp
Two-axis explicit forward-Euler Laplacian stability bound
kappa_h*dt_therm*(1/dx^2+1/dy^2) <= 0.5 — unchanged from the
existing rdb_config.F90 check this module absorbs (only the
length scale changes: the real per-cell metric minimum, not
nominal dx/dy).
real(kind=wp),
private,
parameter
::
MUNK_MIN_CELLS
=
2.0_wp
Minimum number of grid cells the Munk sidewall boundary layer
delta_M = (nu_h/beta)^(1/3) must span; below this the wall
carries grid-scale (2-delta) noise instead of a resolved
boundary-layer profile (see
validation_examples/ocean/acc_channel/acc_channel.nml, where
this exact criterion is documented and was hand-derived).
real(kind=wp),
private,
parameter
::
SIGMA_STIFFNESS_LIMIT
=
0.2_wp
Terrain-following STIFFNESS (slope) parameter bound
rx0 = |H_a - H_b| / (H_a + H_b) <= 0.2 over every face joining
two wet columns, where H is the COLUMN the sigma coordinate
divides into nz layers — under an ice shelf that is the WATER
column b - z_draft, not the bathymetry.
This is the classical σ-coordinate criterion: Beckmann &
Haidvogel (1993), J. Phys. Oceanogr. 23, 1736-1753, §2c, who
introduce r = |Δh|/(2h̄) (algebraically the same number) and
smooth their seamount to r <= 0.2; the “hydrostatic
consistency” condition of Haney (1991), J. Phys. Oceanogr. 21,
610-619, is the same statement. It bounds the σ pressure-gradient
truncation, whose amplitude goes as the CUBE of the interface
offset Δe between neighbouring columns
(a_peak = N²·Δe³/(6·dx·H̄), derived in
validation_examples/ocean/ice_shelf_cavity/README.md), so a
factor 2 in rx0 is a factor 8 in spurious acceleration.
WARNING, never an error: a violated rx0 is not an instability
on its own — with N² = 0 the truncation is identically zero at
any rx0 — and plenty of useful runs are forced hard enough,
damped hard enough, or short enough not to care. What it says is
that the run’s spurious PGF force is NOT small, so a quiescent or
long integration over that geometry will measure the truncation
rather than the physics.
It fires on healthy shipped cases, by design. Measured over
validation_examples/ocean/: 9 of 72 namelists trip it, on four
distinct geometries. The double_gyre"spoon" continental
slope reads 0.348 and neverworld2’s shelf reads 0.893
(seamount_obc_baroclinic sits just over at 0.235), and all of
them run for hundreds of days — because they carry
nu_h = 10000 m² s⁻¹, i.e. a constant lateral-viscosity floor
big enough to arrest a steady spurious force at a/r instead of
integrating it. That is the correct reading of the warning on a
forced configuration, and it is worth saying once at configure.
The three ice_shelf_cavity/ files are QUIET (rx0 ≈ 0.015:
flat bed, and the only tilted boundary is a 13.8 m lid step).
Motivating failure:
validation_examples/ocean/isomip_plus/ocean0_idealised_draft.nml
carries rx0 = 0.73 at the ISOMIP+ trough sidewall (a 23 m water
column beside a 146 m one across one 2 km face) and goes
non-finite at day 3.2 with nothing in the log at configure time.
real(kind=wp),
private,
parameter
::
VISCOUS_CFL_LIMIT
=
0.125_wp
Single-axis viscous-diffusion stability bound nu_h*dt/dx_min^2.
Taken directly from the constant this codebase ALREADY uses at
runtime for exactly this quantity: the bound_kh per-face clamp
(rdb_ocean_horizontal_viscosity.F90) limits the harmonic
viscosity to bound_coef*0.125/(dt*(idx^2+idy^2)), documented
there as “~1/4 of the forward-Euler stability limit” of 0.5 (the
same two-axis sum-form bound rdb_ocean_hdiff_tracer.F90’s
kappa_h check already uses, see DIFFUSIVE_NUMBER_LIMIT
below). This audit collapses that two-axis form to the single
worst axis (nu_h*dt/dx_min^2 rather than
nu_h*dt*(1/dx_min^2+1/dy_min^2)) so the reported number matches
the plain “viscous CFL” a user computes by hand, at bound_coef=1
parity with the runtime clamp’s own margin — i.e. this check trips
at exactly the nu_h that would need bound_kh’s protection.
Two-axis explicit forward-Euler diffusive number
kappa_h*dt_therm*(1/dx_min^2+1/dy_min^2), conservatively
evaluated at the SAME worst-case dx_min on both axes (matches
the “use the minimum cell” instruction; exact on an isotropic
worst cell, strictly more conservative than using the true
per-axis pair). dx_min<=0 returns 0.
Munk boundary-layer width delta_M = (nu_h/beta)^(1/3).
beta<=0 (f-plane — no meridional PV gradient, no Munk
boundary layer) or nu_h<=0 returns huge(1.0_wp) (no
constraint — never trips the >= 2-cell criterion).
Largest dt (s) that keeps ocean_viscous_cfl_number at or
below limit, at fixed nu_h/dx_min — the “reduce dt below
…” half of the audit’s suggested fix.
Largest nu_h (m^2/s) that keeps ocean_viscous_cfl_number at
or below limit, at fixed dt/dx_min — the “reduce nu_h
below …” half of the audit’s suggested fix.
nu_h*dt/dx_min^2 — the single-axis viscous-diffusion stability
number checked against VISCOUS_CFL_LIMIT. dx_min MUST be the
smallest actual cell edge in the domain (e.g.
metrics_dx_min/ocean_stability_min_cell), never a nominal
&grid_nml dx/dy (degrees on non-Cartesian grids). dx_min<=0
(degenerate/unset grid) returns 0 (no constraint expressible).
Arguments
Type
Intent
Optional
Attributes
Name
real(kind=wp),
intent(in)
::
nu_h
real(kind=wp),
intent(in)
::
dt
real(kind=wp),
intent(in)
::
dx_min
Return Value
real(kind=wp)
Subroutines
public pure subroutine ocean_sigma_stiffness_worst(nx, ny, i0, i1, j0, j1, column, wet, rx0_max, i_at, j_at, is_x, h_thin, h_thick, n_over, n_face)
Worst (largest) ocean_sigma_stiffness over every face joining
two WET columns inside [i0,i1] x [j0,j1], with its location,
its two column thicknesses, and how many faces are over
SIGMA_STIFFNESS_LIMIT.
Run all configure-time stability checks. Must run AFTER
configure_ocean_metrics + configure_ocean_land_mask (needs
the real filled ocean_metrics_t), before ocean_state_enter_data.
ierr present -> OCEAN_STATUS_ERR_SETUP on any hard-bound
violation (never error stop); warnings always just log,
whatever ierr does. Rank-0-only logging (mirrors every other
configure_ocean_* info/warning line); the ERROR path itself
always fires (every rank must agree the config is broken).
Reference column thickness (m) at T cells, ghosts included —
bt_work%bt_H_ref, which is b - z_draft afloat and 0
where grounded. Present ⇒ the terrain-following stiffness
check (Check 5) runs; absent ⇒ it is skipped, which is what a
caller with no barotropic datum yet should do.
Worst-case (smallest) delta_M/dx ratio over the physical
domain, honouring a latitude-varying beta under
coriolis_scheme='planetary' on a non-Cartesian grid (beta =
2*omega*cos(lat)/R, maximal — hence delta_M MINIMAL, the
worst case — at the most equatorward row) and a constant beta
(&ocean_topo_nml coriolis_beta) everywhere else. ratio_min is
huge(1.0_wp) (no constraint) when beta<=0 everywhere (an
f-plane run has no Munk boundary layer to resolve).
Smallest actual cell edge over the WET PHYSICAL domain (excludes
ghosts and land), taken over BOTH dxT and dyT, with its (i,j)
location and which axis (is_x) it came from — for actionable
messages (“near j=110”). Host-side, configure time; the grid
sizes here are at most a few 10^5 cells (a global tripolar
config), trivial to scan once.