rdb_ocean_gm Module

Gent & McWilliams (1990) / Griffies (1998) skew-flux thickness diffusion. The eddy-induced (“bolus”) overturning is realized as a layer THICKNESS flux (never an explicit velocity), applied as its OWN sequential operator on the CURRENT thickness, after the dynamics of every outer step (rdb_continuity :: continuity_gm_apply, driven by rdb_ocean_dyn :: run_gm_step) — mirroring how MOM6 applies its thickness_diffuse step after the split-RK2 dynamics step, updating h in place.

Why sequential, not folded

The per-face availability cap h_avail = A·(h − H_VANISHED)/(4·dt) bounds four faces’ outflow by h − H_VANISHED, so h_new >= H_VANISHED — but ONLY for the h it was computed from. Until 2026-10 this slot computed uhD from the stage-entry h and FOLDED it into the resolved continuity sweeps (gm_fold_x/y), so the cap bounded the stage-entry thickness and the resolved outflow came on top: on the 1-degree Southern Ocean (z*, open steps) an 8.6 cm partial bed cell with two open faces onto fillers was drained by GM at the cap on all four faces and the resolved flux took it to −8.2e-4 m. Computed here from the thickness the dynamics LEFT, the cap bounds what is actually there (prototype: python_prototypes/gm_sequential/gm_sequential.py).

At each u/v face the eddy streamfunction is Sfn_unlim = -(KhTh * dy_Cu) * Slope using the STORED isopycnal slope from the slopes slot. The per-layer bolus transport is the vertical difference of the limited streamfunction, formed by the column recurrence (surface->bed, uhtot=0 at surface): uhD(k) = max( min(Sfn_in_H - uhtot, h_avail_i), -h_avail_{i+1} ) uhtot = uhtot + uhD(k) so Sum_k uhD = 0 (closes at the bed) — mass/tracer conservative.

Bottom-blocking (MOM6’s rule for the GM streamfunction) acts on the UNLIMITED streamfunction first: no transport from a donor layer that lies entirely below the receiving column’s bed, and a share scaled by the fraction above it for a donor layer that straddles it (gm_block_below_bed). On an open z-level step this is what keeps GM from pouring deep water through the step into the fillers below the shallow column’s bed.

Three limiters, in order: (1) safe-streamfunction blend toward a column-spread return flow where slope > slope_max; (2) mass- availability rsum bound (the conservation guard keeping each layer

= H_VANISHED); (3) per-layer donor cap. The donor side is keyed on the SIGN of uhtot (column i when uhtot<=0, else i+1).

KhTh is a 2D face field, constant-filled from &ocean_gm_nml khth (VarMix/MEKE seam populates it later via +=), CFL-clamped per face. gm_src(nx,ny) carries the GM PE release (>= 0 for a stable tilted column) for future MEKE coupling.

Bottom-up convention (k=1 bed, k=nz surface; interface K=1 bed, K=nz+1 surface); the slopes slot zeroes the slope at both caps so the streamfunction vanishes there and the recurrence closes.

Partial-step z-level faces (&vcoord_nml zfixed_closed_faces)

Under z_fixed a face column is not the whole water column: a layer that is an inert FILLER on either side (inside the bed or the ice draft) is a WALL for that layer at that face (metrics%open_u/open_v == 0). Continuity applies that mask to the resolved flux, and GM runs on its own (it does not pass through that mask), so GM must build its overturning on the OPEN part of each face column itself. With the knob on, each face first marks its open layers, ok(k) = open(k) .and. live(h_W(k)) .and. live(h_E(k)) (rdb_vl_is_live, the one vanished-layer predicate), and the recurrence then:

• gives a not-ok layer ZERO transport and zero availability (h_avail = 0, so it adds nothing to rsum): the streamfunction is carried UNCHANGED across it, so it is 0 at the BOTTOM of the open column (it starts at 0 at the bed, and the bed-side closed / filler layers cannot change it);
• closes the column into the TOPMOST open layer k_top_open (uhD(k_top_open) = -uhtot) instead of k = nz, so the streamfunction is also 0 at the TOP of the open column — the free surface, or the ice base with fillers above it. The full-column closure into nz would push the residual through a filler under the ice, or through a face whose top layer is closed.

Hence uhD == 0 on every closed face-layer and every filler, and Sum_k uhD = 0 at every face (the open-column integral IS the column integral). With nothing closed ok is all-true, k_top_open = nz and the arithmetic is the full-column form operation for operation. The slopes slot zeroes slope / N^2 at every interface not strictly inside a face’s open column on this path, so gm_src (the MEKE source) sees only the open column as well. Knob OFF ⇒ use_open = .false., ok all-true, the masks are never indexed ⇒ byte-identical.

Non-finite guard: the streamfunction limiters are min/max clamps, which nvfortran’s relaxed FP lowers NaN-blind (CLAUDE.md) — a NaN slope would come out as a plausible ±rsum transport. A non-finite slope is therefore read as ZERO slope (no GM at that interface) before any clamp sees it; finite inputs are untouched.

Default off (&ocean_gm_nml enable=.false.) => slot allocated but gm_compute_transports no-ops => bit-identical. Refs: Gent & McWilliams (1990); Gent et al. (1995); Griffies (1998).