Procedures

accel_visc_rem post-apply reweight: vel = snap + rem·(vel − snap) — the whole explicit-tendency sum accumulated since the snapshot is attenuated by the per-layer viscous remnant (linearity ⇒ identical to weighting each tendency individually, MOM6 u = u_init + dt·visc_rem· (CAu + PFu + diffu)). Faces with rem == 1 (the unconditioned source=1.0 init, and every face before the first vdiff fills the producer) are SKIPPED, not rewritten — snap + 1·(vel−snap) is not an FP identity, and the skip keeps rem≡1 bitwise inert. Friction-dominated near-massless layers (rem → 0) keep their entry velocity. Public only for the unit-test suite. Under mem:separate this do concurrent runs on the device-resident arrays in production (all mapped on the ocean state); unit tests must map their own arrays explicitly. The masked write stays inside the do concurrent body (legal — no cross-iteration dep).

accel_visc_rem stage-entry snapshot: snap = vel, device-side. Public only for the unit-test suite. Under mem:separate this do concurrent runs on the device-resident arrays in production (both are mapped on the ocean state); unit tests must map their own arrays explicitly (!$acc enter data copyin / update self).

Accumulate one RK2 stage’s zonal mass flux into the windowed accumulator with the ½ RK2 weight baked in: uhtr += 0.5·mass_flux_x·dt. mass_flux_x is already area-weighted (m³/s = u·h_face·dy_cu); the product is m³.

Meridional analogue of accumulate_flux_x.

Adcroft reciprocal: 1/x, but 0 -> 0 (zero-width faces give zero inverse, no NaN/Inf). Single source for every metric inverse (D4).

Register the four flat knobs of one &ocean_dataovr_nml tag. e is target, intent(in) for the same reason cfg is in every register_ocean_* above: the schema stores a pointer to the target and assigns through it at parse time, and the actual argument is always a component of the target cfg, so the association outlives this call.

Append a message to a growable error list.

Add the ice-shelf top-drag tendency into the already-summed slow forcing. Separate from sum_slow_tendencies_into_F_slow (rather than a sixth term in it) for one reason: the top-drag slot is OPTIONAL all the way down the driver chain, and the sum above must stay a single unconditional kernel with no present branch inside its do concurrent.

.true. iff the anisotropy-direction mode has an implemented direction tensor. Only mode 0 (grid-relative (n1,n2) = aniso_dir, a constant tensor) is built; the MOM6 flow-aligned modes are not ported. Drives the configure-time fail-loud guard in validate_config so a requested-but-unimplemented mode aborts the run instead of silently falling back to the grid-i default.

ieee_is_nan OR-reduced across the five stats scalars.

Append one value token to the values buffer.

Replace the bt mode in the per-layer face velocities with the barotropic-substep end-step value, adding Δu·wt_k to every layer, Δu = u_bt_end − u_bt_at_n − dt·F_bt_u (same for v). Split-explicit convention (Hallberg 2009): momentum uses the END-of-step barotropic velocity; layer continuity earlier used the time-mean transports. Both legs of the corrector use the same end-step anchor (mismatched anchors overshoot the gravity-wave phase speed). Also rescales h_layer uniformly so the column total matches H_ref + η_end. hTr is deliberately NOT rescaled (would break exact tracer mass conservation; T = hTr/h drifts by O((η_end−η*_slow)/H) per step).

Post-parse Cartesian domain sizing (MOM6 GRID_CONFIG=”cartesian” + AXIS_UNITS / LENLON / LENLAT). When &ocean_grid_nml len_lon/len_lat are set (> 0) on an ocean Cartesian grid, DERIVE the uniform &grid_nml dx/dy (metres) from the domain extent, the physical cell count (nx/ny = MOM6 NIGLOBAL/NJGLOBAL) and axis_units, exactly as MOM6’s set_grid_metrics_cartesian: degrees : dx = rad_earth · len_lon · π/180 / nx (arc length, no cos(lat)) km : dx = 1000 · len_lon / nx meters : dx = len_lon / nx and likewise dy from len_lat / ny. This runs BEFORE grid init, so the metrics, barotropic CFL and Coriolis all see the derived metres. Default (len_lon <= 0) leaves dx/dy untouched ⇒ bit-identical.

Configure the ocean diagnostic set from the unified &ocean_diag_nml diags selection string. Single production entry point: parses spec once, registers the canonical defaults (consulting the spec for per-diagnostic :off skips and :cadence / :op / :coord overrides), then registers any spec entries that name a derived-catalog diagnostic (with their overrides).

Apply the pre-factored tridiagonal (from build_factorize_tracer_matrix) to one tracer: form T = hTr/h, run the Thomas RHS sweep + back-substitution against the stored pivots, and reconstitute hTr = T_new * h_layer. Operates on concentration so it conserves cell-centred T; h_layer and the factored a/b/c are read-only (shared across every tracer).

Stamp src * wet_mask(i,j) onto the lowest massive layer of a tracer’s hTr array (first k with h_layer > h_min, scanning k = k_bot(i,j)..nz from the first live layer up), mirror into the matching budget contributor. Flat-impl over plain allocatables — the outer subroutine reaches ms%tracers(idx)%hTr on the host before calling this.

Public only for the unit-test suite (no production module imports it); ignore when developing production code in other modules. Helper: write per-layer density from layer_rho_init(:) when the caller has set at least one non-sentinel value. Sentinel is < 0; we expect exactly nz_ml non-sentinel entries listed in k=1..nz_ml order (bed → surface).

Truncate face velocities to |u| ≤ maxvel. MOM6’s MAXVEL analogue. Called once at the end of each outer step (after the RK2 average so the clipped state is what gets carried into the next stage). No-op when maxvel <= 0.

Set south and north wall-face per-layer v. Mirror of apply_zonal_baroclinic for the y-direction.

Post-anomaly nudging for the north edge. “Incoming” at north = v_wall < 0 ⟹ tau_in.

Post-anomaly nudging for the south edge. “Incoming” at south = v_wall > 0 ⟹ tau_in.

Post-anomaly nudging for the east edge. “Incoming” at east = u_wall < 0 (westward, into domain) ⟹ tau_in.

Post-anomaly nudging for the west edge. “Incoming” at west = u_wall > 0 (eastward, into domain) ⟹ tau_in.

Consume one key = value [value...] starting at token ti; advance ti past the consumed tokens. Validates structure, unknown keys, duplicate keys, and applies via parse_tokens.

Orlanski radiation for the east open edge (Orlanski 1976). Interior face index I = i_e-1 (1st), I-1 = i_e-2 (2nd). Outward normal = +x. dhdx = u(i_e-1) - u(i_e-2) (eastward gradient).

Orlanski radiation for the north open edge (Orlanski 1976). Interior face J = j_n-1 (1st), J-1 = j_n-2 (2nd). Outward normal = +y. dhdx = v(j_n-1) - v(j_n-2) (northward gradient).

Orlanski radiation for the south open edge (Orlanski 1976). Interior face J = j_s+1 (1st), J-1 = j_s+2 (2nd). Outward normal = -y. dhdx = v(j_s+1) - v(j_s+2) (southward gradient).

Orlanski radiation for the west open edge (Orlanski 1976). Interior face index I = i_w+1 (1st), I-1 = i_w+2 (2nd). Outward normal = -x. dhdx = u(i_w+1) - u(i_w+2) (westward gradient). rx updated in-place; u_prev NOT updated here (done by the caller after all edges are set).

Ice-shelf-cover twin of apply_surface_restore_2d_impl: the open-water factor 1 - cover_frac composes multiplicatively with wet_mask, so a covered column receives EXACTLY zero restoring — and, because the same factor multiplies the budget mirror, exactly zero restoring shows up in the heat/salt surface budget there too. Separate _impl, not an in-loop present() test (house idiom).

Stamp the per-step restoring increment dt·p·(target - surf)· wet_mask onto the top layer (k = nz) of a tracer’s hTr array and mirror it into the matching budget contributor. Explicit-shape dummies so NVHPC stdpar compiles device kernels against static bounds.

Wet/dry variant of apply_surface_src_2d_impl: the DYNAMIC cell wet mask composes multiplicatively with the static one, so a dynamically dry column (total depth below &ocean_wetdry_nml dry_depth) receives NO surface flux — heating a mm-scale residual sliver would blow its temperature up (docs/ocean_wetdry_plan.md §4.4). Separate _impl (not an in-loop optional test): the knob-off path keeps the original kernel untouched, byte-identical.

Wet/dry NOBUDGET twin — see apply_surface_src_2d_nobudget_impl and apply_surface_src_2d_dyn_impl.

Stamp inv_scale · Q_field(i,j) · wet_mask(i,j) onto the first LIVE layer (k_top(i,j)) of a tracer’s hTr array, mirror into the matching budget contributor. Explicit-shape dummies so NVHPC stdpar can compile device kernels against static bounds.

Byte-for-byte copy of apply_surface_src_2d_impl with the budget dummy and its accumulation line removed — the pseudo-salt mirror of salinity’s surface flux, which by contract (budget_id = TRACER_BUDGET_NONE) must not touch salt_budget_surface. Separate _impl, not an in-loop present(budget) test (house idiom, see apply_surface_src_2d_dyn_impl’s docstring) — this keeps the production S/T impl untouched and the increment hTr receives bit-identical to salinity’s.

The two cell-centred surface kernels that do NOT route through the assembler’s Q_heat / Q_salt, with their ice-shelf-cover dispatch. Shortwave penetration reads a pristine q_sw component (or moves a lump the masked deposit never added) and restoring forms its flux in-kernel from the live SST/SSS, so each needs the cover factor of its own; everything else the atmosphere contributes is already masked inside ocean_surface_flux_assemble.

Ice-shelf-cover twin of apply_sw_penetration_impl: the open-water factor 1 - cover_frac composes multiplicatively with wet_mask into the column irradiance I0, so a fully covered column neither removes the surface lump nor deposits a profile — it is left EXACTLY untouched. Separate _impl, not an in-loop present() test (house idiom, see apply_surface_src_2d_dyn_impl) — the cover-off path keeps the original kernel byte-identical.

Per-column two-band shortwave redistribution. Explicit-shape dummies so NVHPC stdpar compiles device kernels against static bounds. Difference form (I(d_top) - I(d_bot)) — no division by h, so a vanishing layer (h → 0 ⇒ d_top == d_bot) absorbs zero automatically with no guard. The transmission T(d) is the shared sw_transmission (same-module ⇒ inlined by NVHPC), so the deposition and the boundary-layer coupling cannot diverge.

Post-parse time-unit fixup shared by the file + string config paths. t_end, status_interval, ocean_diag%dt_out are given in time_unit from &time_nml; multiply through to SI seconds. dt_fixed/dt_max stay in seconds; default time_unit="s" gives a 1× factor (bit-identical to before).

Post-RK2 velocity housekeeping: the advective-CFL truncation (E7) followed by the absolute maxvel cap. Called once at the end of each outer step (after the RK2 average, before the ALE remap) so the carried-forward / remapped velocity field is bounded. Replaces the bare apply_maxvel_clamp call at both ocean_dyn_step / ocean_dyn_step_split sites.

h_layer(k) += dt * (w(k) - w(k+1)) per cell. Bed and surface interfaces feed through whatever w the caller set (zero for the bed by default after compute_w_from_continuity with enforce_bed_bc = true). No h-floor — plain accumulation; the caller is responsible for guarding against negative thickness if the prescribed w + dt is large.

Set west and east wall-face per-layer u to Flather mean + zero-gradient baroclinic anomaly (radiating edges) or uniform clamped_u (CLAMPED edge). Explicit-shape dummies; one do concurrent per edge (j outer, k inner).

Byte footprint of an allocatable array (0 when unallocated), so a gated-off slot naturally contributes nothing to a counted total. The COUNTED companion to the runtime’s MEASURED free-memory delta: each *_state_t sums its own arrays through this, the driver reports the total before enter_data (works on CPU builds too, where there is no device query), and mem_log_device_actuals reconciles the count against the measured mapping — a state array added without a matching bytes() term makes the measured map exceed the count and self-announces the drift (the anti-rot guard for the otherwise rot-prone per-slot count).

Copy the target field into h_layer on grounded columns only — non-grounded columns’ h_new was never written, and their old h_layer is already the exact target.

Counted allocatable footprint of the ocean C-grid barotropic slot.

Tear down host allocations. Call exit_data first — destroying a still-device-mapped state leaks device memory silently.

Attach the C-grid barotropic allocatables to the device. The parent struct is mapped by the orchestrating routine. Read+write fields use copyin; pure-workspace fields use create.

Reverse of enter_data — copy out prognostic fields, drop scratch. Caller detaches the parent struct after this returns.

Allocate every C-grid barotropic array zero-filled. East-face arrays are (nx+1, ny), north-face (nx, ny+1). Scalars (manning_n, coriolis_f) are populated by state_init_from_config afterwards.

Public only for the unit-test suite (no production module imports it); ignore when developing production code in other modules. Linearized barotropic-substep kernel. Forward-Euler steps the barotropic state (eta, ubt, vbt) at dt_inner for n_steps against the constant slow forcing force_u, force_v (each at the same C-grid location as ubt, vbt). Accumulates the per-step running sums into ubt_sum, vbt_sum, eta_sum; at the end divides by n_steps and stores the time-mean back into bt_eta, bt_ubt, bt_vbt for the caller.

Nonlinear barotropic substep. Same time-mean accumulator pattern as barotropic_substep_linear but the per-substep dynamics include the three barotropic nonlinearities that matter for MOM6-grade physics:

Interior (normal-width) entry point for the nonlinear barotropic fast loop. Unpacks the bt_work fast-loop arrays and forwards them to barotropic_substep_nonlinear (bt_halo=0), so the ~20-array plumbing lives here once instead of at every call site. The optional bc/t/eta_forcing propagate by absence (F2018 15.5.2.13), so this single entry reproduces the former present()-branch call variants bit-for-bit. See bt_wide_substep for the wide-halo march-in twin.

Counted allocatable footprint of the barotropic fast-loop work state (BTCL_u/v derived-type coeffs excluded) slot (0 when unallocated).

Attaches component arrays only — no bare copyin(this) (the polymorphic stack descriptor caused AMD libomptarget cross-slot overlap). The workstate descriptor reaches the device via the root copyin(state%ocean) (bt_work is inline all the way up); the per-slot copyin was a competing second mapping that cost ~36% of the fast loop.

Allocate the 2D barotropic-substep arrays. Pass nz_ml to also allocate the split-driver slow-tendency accumulators; omit when only the 2D barotropic substep is needed (unit tests).

Translate a raw nc_check-style status (0 = ok) from one of the nc_* reader calls in load_bathymetry_into_array into the caller’s ierr contract: .true. on success; on failure, .false. with ierr = OCEAN_STATUS_ERR_IO when ierr is present (closing ncid first, when given, so a mid-read failure does not leak the file handle), or error stops with the SAME generic text nc_check itself would have used had ierr never been threaded through — keeps the legacy (no ierr) behaviour byte-identical while unblocking the ierr-present return path (F1/F2 of the P0.1 review).

Fill ghost-cell bathymetry by constant extrapolation from the nearest interior cell. Deliberate duplicate of rdb_bathymetry::fill_bathymetry_ghosts_array — see the module docstring for why this module cannot use that one.

Approximate median via a full sort — b is a setup-time array (called once per create(), never per-step), so O(n log n) is fine; this exists purely to make the sign-error message name a representative depth rather than minval/maxval (which a single outlier cell would distort).

In-place sign-normalise b (any shape — interior or full array, the caller decides what it passes) to Roundabout’s positive-down depth convention, then validate on the NORMALISED array: zero wet cells (b > LAND_DEPTH_THRESHOLD) is rejected as OCEAN_STATUS_ERR_BATHYMETRY_SIGN, naming the median depth and the convention that was requested, so the caller can see at a glance that the OTHER convention was probably meant. NOT a “no negatives” check (see module docstring) — a majority-negative but non-empty wet fraction is accepted (legal under wet/dry).

Cell concentrations of one tracer through THE vanished-layer column rule (rdb_vl_column_conc): hTr/h on a live layer, the donor’s concentration on a filler, never a quotient by a near-zero or negative thickness.

Per-face body of vdiff_set_viscous_bbl (see there). x_face selects the u-faces (nx+1, ny) (cells i-1, i) or the v-faces (nx, ny+1) (cells j-1, j).

Will &ocean_vdiff_nml bbl_glue actually be ON after setup?

.true. only for a bottom-drag variant with a real kernel (BDRAG_LINEAR / BDRAG_QUADRATIC). BDRAG_INVALID returns .false.. The single gate validate_config consumes (PR-6).

Map a logical cap location to geographic (lat, lon) in degrees. lam_deg : pseudo-longitude (geographic lon at the join ring), any real (wrapped internally); the i-direction. s : cap-row fraction, 0 = join ring (lat = phi_join), 1 = fold line; the j-direction. phi_join: join latitude (deg); cap covers lat > phi_join. lon_pole: longitude (deg) of the first cap pole; partner at lon_pole + 180.

Latitude (deg) of the two cap poles — they sit ON the join latitude, so this is simply phi_join.

.true. when both background schemes are selected at once.

HORIZONTAL (cross-face) Boole quadrature of the layer pressure increment dpa = g * int rho' dz — the face integral the FV pressure-gradient contour needs (Adcroft, Hallberg & Harrison 2008 §3; Yung, Hallberg, Adcroft & Morrison 2026 §2.4).

The cross-face Boole quadrature of boole_dpa_face for a CONSTANT-BY-LAYER (PCM) T/S column, with MOM6’s near-bottom mass-weighting of the interpolated T/S (MOM6 int_density_dz_generic_pcm, intx_dpa).

boole_dpa_face with the Wright (1997) vertical rule boole_dpa_intz_layer_wright at each of the five sub-columns – identical sub-columns, weights and summation order, no eos_t handle. 15 inline Wright density evaluations per face per layer (the end points are the columns’ own dpa).

5-point Boole-quadrature density-anomaly integral over one layer (Adcroft, Hallberg & Harrison 2008; White, Adcroft & Hallberg 2009). Returns the layer-integrated pressure-anomaly increment dpa = g * int rho' dz and its first moment (from the layer TOP edge inward) intz_dpa = 0.5 * g * dz^2 * bracket.

boole_dpa_intz_layer specialised to Wright (1997): the same five sub-points and weights, with the density written inline (wright_rho, the expression eos_density_point evaluates) – no eos_t handle, no per-point variant dispatch. The caller selects this twin ONCE, outside its loops.

Boole weights of the five sub-point density anomalies r5 (top to bottom): dpa = g*dz*<rho'> and the first moment from the top.

The five sub-point (T, S, p) triples of the in-layer Boole rule, top (n = 1) to bottom (n = 5), for the per-EOS twins – the same points boole_dpa_intz_layer writes out in place.

Linear-exact one-sided edge pair for a BOUNDARY layer (k=1 or k=nz), where a centred slope has no second neighbour.

Linear-exact half-jump across a BOUNDARY cell (k=1 or k=nz), where a centred stencil has no second neighbour.

Smallest barotropic substep count n_inner such that the substep dt_outer/n_inner satisfies the 2-D external-gravity-wave CFL: dt_bt <= cfl_safety * l_cfl / c_ext, with l_cfl the 2-D CFL length (metrics_bt_cfl_length) and c_ext the external wave speed sqrt(g*H_max). MOM6 set_dtbt analogue, now with the cross-direction term included (the legacy estimate used a 1-D length and under-counted n_inner by ~sqrt(2) on square cells, leaving the effective 2-D CFL at ~0.92 for cfl_bt_safety=0.65 — on the edge of the forward-backward scheme’s stability).

Smallest n_inner >= 1 with dt_outer/n_inner <= dt_bt_safe, for an ALREADY-LIMITED safe barotropic substep dt_bt_safe (s) — the per-wet-cell minimum bt_cfl_dt_wet returns, reduced across ranks. bt_auto_n_inner is this with dt_bt_safe formed from a single (c_ext, l_cfl) pair.

Per-WET-CELL external-gravity-wave CFL limit (MOM6 set_dtbt):

Single source of truth for the BT march-in exclusion set (mirrors the explicit-bt_halo > 0 fail-loud checks in validate_config). Reports whether ANY exclusion is active and names the first one (for the AUTO resolution log). When none is active, reason is “”.

compute_bt_rem under &vcoord_nml zfixed_closed_faces: the same H/(H + r·hbbl·dt_inner) with H = Σ_k h_face·open (the OPEN-column centred face depth, the weight face_depth_mean_* uses). A face whose every layer is closed has H = 0 and keeps bt_rem = 1, exactly like a dry face on the original path (it carries no barotropic transport: dy_cu_bt = 0 there).

compute_bt_rem_wave_drag under &vcoord_nml zfixed_closed_faces: MULTIPLIES H/(H + r_H·dt_inner) into bt_rem with the OPEN-column face depth H = Σ_k h_face·open (the bt_rem_open_impl depth). H <= 0 (every layer closed) ⇒ unmodified, MOM6’s guard.

Counted allocatable footprint of the wide-halo BT shadow state (0 when unallocated, i.e. whenever &ocean_bt_nml bt_halo = 0).

Offset-copy normal-width input arrays into the wide shadow arrays. Dispatches to the non-polymorphic _impl body to avoid the class-box GPU descriptor issue (same pattern as enter/exit_data).

Non-polymorphic copy_in body. Offset = bt_halo: w_X(iw, jw) = X(clamp(iw-off), clamp(jw-off)) over the FULL wide extent — the inner band is a direct offset copy; the outer bt_halo ring is a clamped-index (constant-extrapolation) fill. The ring fill matters: without it the ring carries stale end-of-fast-loop values from the previous stage (H_ref = 0, eta from t-1), which at a PHYSICAL (non-seam) edge is never refreshed by any exchange and free-runs an inconsistent zero-depth integration that blows up in O(25) outer steps. At an MPI seam the ring is immediately overwritten with true neighbour data by entry_exchange, so the clamped fill only governs physical edges — the same sane ghost-band construction the v1 normal-width path gets from its own ghosts. Scratch / accumulator arrays (w_eta_new, w_ke, w_eta_sum, …) do not need copy-in — the substep initialises them.

Offset-copy wide output arrays back to the normal-width arrays. Dispatches to the non-polymorphic _impl body.

Non-polymorphic copy_out body. X(i,j) = w_X(i+off, j+off) over the full normal index range. Outputs: time-mean eta/ubt/vbt, uhbt/vhbt, *_end snapshots.

Deallocate all wide state.

Attach all wide arrays (and wide metrics leaf arrays) to the GPU present table. The containing ocean_dyn_t is already mapped by the caller; this routine attaches the components.

Non-polymorphic enter_data body (avoids class-box GPU descriptor issue).

One wide grouped exchange (eta+ubt+vbt) + wide singles for the other 7 input arrays (H_ref, ubt_prev, rem_u, force_u, vbt_prev, rem_v, force_v). Fills the entire wide ghost band before the fast loop. Counter effect: +1 bt_group, +1 centre_2d, +3 face_x_2d, +3 face_y_2d. Dispatches to the non-polymorphic _impl body.

Non-polymorphic entry_exchange body.

Detach all wide arrays from the GPU present table.

Non-polymorphic exit_data body.

Allocate the wide shadow state. Builds grid_w (same nx_phys/ny_phys as grid, nghost = grid%nghost + bt_halo), fills wide metrics via the same formula generator, fills wide f_corner. grid_config must be GRID_CONFIG_CARTESIAN or GRID_CONFIG_SPHERICAL; supergrid/tripolar are excluded at configure time.

Wide-halo (march-in) entry point for the nonlinear barotropic fast loop. Unpacks the wide shadow arrays (w_*, wide grid/metrics, wide f_corner) and forwards them to barotropic_substep_nonlinear with bt_halo = bt_wide%bt_halo, so the ~20-array plumbing lives here once rather than at the call site. bc propagates by absence. Tides (eta_forcing) are a configure-time exclusion on the wide path, so none is forwarded. Interior twin: barotropic_substep_nonlinear_interior.

Per-cell T-area (m²) for physical integrals: cached metrics areaT once set_area has run, else uniform grid%dx·grid%dy.

Total KE: Σ over masked interior of 0.5·h·(u²+v²)·areaT·weight, u/v averaged from C-grid faces to centres. Public only for tests.

Total mass: Σ over masked interior of h_layer·areaT·weight over k. Units: m³. Public only for the unit-test suite.

Total tracer content: Σ over masked interior of hTr·areaT·weight over k. hTr is thickness-weighted, so this is ∫(tracer·volume), the conserved quantity. Public only for the unit-test suite.

Build the backward-Euler tridiagonal per cell column and run the Thomas forward factorization — the tracer-INDEPENDENT half of the vertical-diffusion solve (depends only on kv/h/dt). Run once per stage; every registered tracer then reuses the factored coefficients via apply_factored_tracer.

ONE coalesced pass: mask(i,j,1) = 1.0 iff any layer of column (i,j) is strictly below the floor; n_grounded counts them (explicit OpenACC reduction — a sum() on a present-mapped array would run host-side under NVHPC non-managed mode and read the stale shadow).

Shared prefix of rdb_ocean_create_from_string AND rdb_ocean_create_pending (P2.5): allocate a handle, parse + validate the config, check dt_fixed. Does NOT touch g_handle_live — callers set it themselves once they know which of the two flows they are in (immediate complete_ocean_create vs staying pending for geometry injection).

Register the validated groups + the still-external groups onto schema, capturing defaults from cfg.

target. Non-grounded columns are SKIPPED — their h_new is never read downstream (every consumer falls back to h_old via the mask).

Build layer-mean T/S concentrations (c = hTr/h) from extensive tracer content + pre-remap thicknesses, for the VCOORD_RHO density inversion. Per-layer rdb_vl_conc (a filler reads hTr/h, its donor’s concentration by I1′; a zero-thickness layer reads 0). Flat-impl, explicit-shape; one cadence-bounded launch per remap.

Layer-centre GEOPOTENTIAL depths (positive-down from the z = 0 datum) from a column of layer thicknesses. Bottom-up: k=1 bed, k=nz_ml surface. z_ctr(k) = z_top + sum_{k'=k+1..nz_ml} h(k') + 0.5*h(k).

Convert a bind(c) character(kind=c_char) buffer + explicit length into a Fortran allocatable string. No null-termination assumption — c_len is authoritative (matches the recovered precedent’s convention).

Name of canonical-catalog entry i (1-based). Public for the P7 discoverability C ABI and the unit-test suite.

Number of names in the canonical-diagnostic catalog (the static set register_default_diags MAY register — some entries are gated, see canonical_diag_gate_hint). Bounds for canonical_diag_catalog_name’s index argument. Public for the P7 discoverability C ABI (rdb_ocean_canonical_catalog_size) and the unit-test suite.

The namelist gate whose closure suppresses canonical diagnostic name, for the “you asked for this and did not get it” warning. "" for the four unconditional diagnostics (SSH/u/v/KE) — a hint of "" means “this one should have registered; that is a bug, not a config”. MUST be kept in lock-step with register_default_diags’s gates; the count of non-empty hints must equal the number of conditionally-registered canonical diagnostics (test_ocean_diag/diag_gate_hint_covers_every_gate is the lock — any PR that adds a new gated canonical diagnostic must add its hint here or that test fails). This duplicates the gate NAMES only, never the gate LOGIC — the single if in register_default_diags remains the one place the gate is evaluated.

Comma-separated list of catalog diagnostic names (for error text).

NO ICE OVER LAND (design rule R2): force z_draft = 0 on every column the bathymetry already calls land (b < LAND_DEPTH_THRESHOLD), and report how many were touched so the caller can log it.

Convert ONE shelf-box bound from metres to grid coordinate units, leaving the CAVITY_BOUND_INF “no limit” sentinel alone. Without the guard the sentinel would be scaled by the degrees-per-metre factor on a spherical grid and come out as a finite (if absurd) bound — harmless numerically, but it would stop meaning what it says, and the next reader would have to re-derive that.

Interfacial buoyancy flux B_b (m^2/s^3), NEGATIVE = stabilising — Yung et al. (2025) eq. (6) p. 5831,

The volume_compensation = "uniform_open_ocean" sink: remove dw metres of the top layer from every wet cell the ice does NOT cover, the removed parcel carrying that cell’s own T and S so no concentration there changes.

Apply the compensation sink’s top-layer ratio to one PASSIVE tracer’s load. Unconditional (scale = 1 off the sink), so there is no mask branch inside the kernel.

The uniform per-unit-area thickness the uniform_open_ocean sink removes: vol_melt/area_open, and exactly zero when there is no open ocean to remove it from (a fully ice-covered domain — the sink then does nothing rather than dividing by zero, and the volume stays in, which the mass budget reports honestly as a growing total).

Count the INTERIOR columns the cavity grounds: wet bed (b >= LAND_DEPTH_THRESHOLD) but water thickness b - z_draft < h_min. Interior-only, because the ghost band carries extrapolated bathymetry and would bias the fraction the grounded_max_frac sanity bound is taken against.

Count cells whose top-of-column pressure does NOT carry the isostatic ice load — the melt path’s guard that configure_ocean_cavity (the sole ms%p_top producer) ran, and ran before this check.

Count ice-covered WET columns sitting at exactly f = 0 — the configure-time domain check for exchange_law = "hj99".

The BAROTROPIC DATUM: bt_H_ref = b - z_draft on a column that has water under the ice, and exactly 0 on one that is GROUNDED (b - z_draft < h_min, i.e. land by the very rule seed_wet_mask_impl applies).

Max violation of the counted-once invariant (I) in METRES of reference depth, over the WET columns: max |bt_H_ref - (b - z_draft)| where b - z_draft >= h_min.

a - b with the cavity missing-value convention.

Normalise a freshly loaded draft field onto Roundabout’s convention (DEPTH, positive down, >= 0).

Configure-time guard: every draft entry is finite and >= 0.

Is the configured ice-shelf draft UNIFORM over the whole array?

Exchange-velocity dispatch with the Coriolis parameter taken from the bundle (par%f_cor). Thin wrapper over cavity_exchange_velocities_f, which holds the dispatch.

Exchange-velocity dispatch — ONE argument list for every law, so a later do concurrent kernel dispatches with a single select case and no reshaping. l_plus is the viscous Obukhov scale of the CURRENT iterate: the single scalar that carries the stratification feedback for every implicit law. Pass CAVITY_L_PLUS_NEUTRAL for the neutral / unsuppressed evaluation.

Thickness-weighted mean of (T, S, u, v) over far_depth METRES below the ice base, with a PARTIAL last layer.

v1 ice-cover fraction: BINARY, 1 wherever there is any draft. An area-blended calving front (a partially covered cell) is a melt-physics decision, not a geometry one, so it is deliberately left for the slice that needs it — the field exists now only so that slice does not have to re-open the metrics lifecycle.

p_ice_ref = (rho_ref*GRAVITY) * z_draft (Pa).

Fill the two OWNED surface-flux components from the solved interface. FULL OVERWRITE of the whole plane, never += — the same rule the sea-ice coupler’s fillers follow, and for the same reason (a += ratchets across outer steps with no bound).

Holland & Jenkins (1999) eqs. (14)-(18) p. 1792:

Yung et al. (2025) “StratFeedback”, eqs. (7)-(8) p. 5832:

The three heat fluxes of (E2), W/m^2:

(c_i_eff, kh) for the requested ice-conduction mode — the two numbers through which every mode enters the quadratic.

Viscous Obukhov scale L+ = L/delta_nu with delta_nu = nu/u*, i.e. L+ = -u*^4/(nu*kappa*B_b) — Yung et al. (2025) eq. (5) p. 5831; the same definition in Vreugdenhil & Taylor (2019) eq. (27) and Rosevear et al. (2022) eqs. (6)+(8) p. 2592. POSITIVE for melting. Returns CAVITY_L_PLUS_NEUTRAL for a vanishing buoyancy flux.

Is this trial L+ the neutral / unsuppressed limit?

Does this law’s (gamma_t, gamma_s) depend on the interfacial buoyancy flux — and therefore on the melt rate it produces? Yung et al. (2025) p. 5833 states the consequence: “Since the transfer coefficients depend on L+, which in turn depends on melt rate via surface buoyancy forcing, iteration is required for convergence of the three-equation parameterisation solution.”

Solid-ice thickness rate (m/s) from the canonical mass flux, m_ice = m_mass/rho_i. REPORTING ONLY — this is Jenkins, Nicholls & Corr (2010)’s a_b convention.

Freshwater-equivalent thickness rate (m/s) from the canonical mass flux, m_weq = m_mass/rho_fw. REPORTING ONLY — this is ISOMIP+’s m_w (Asay-Davis et al. (2016) eq. (24) p. 2485), the number their figures are in.

Copy a 2-D cavity field into the diag buffer, writing IEEE NaN wherever the column was not solved. Every cavity fill that is a plain read of a slot array routes through here, so the missing-value convention cannot drift between them.

The real-mass top-layer source: add the meltwater VOLUME to h_layer at the first LIVE layer k_top(i,j), replace the virtual salt flux the assembler stamped with the real advective salt w*s_ice, and add the enthalpy dh*T_b — mirroring both tracer increments into the existing surface budget contributors.

The pseudo-salt mirror of cavity_mass_apply_impl’s SALT increment, with no budget accumulation — budget_id = NONE, so it must not touch salt_budget_surface.

Is this step’s clamped-withdrawal tally a FAIL-LOUD condition?

The two INTERIOR integrals the mass budget and the compensation sink need, reduced on device in one pass:

Data-parallel driver: solve n independent columns.

Masked 2-D driver: form the friction velocity and solve the interface on every ICE-COVERED column of an (nx, ny) plane, leaving the rest untouched at exactly zero.

Scalar entry point returning only what a coupling seam consumes: the interface state, the canonical melt mass flux and the ocean -> interface heat flux. A thin wrapper over cavity_solve_melt that keeps the solution BUNDLE inside the callee, so a do concurrent over columns needs no derived-type local(...) clause at all — each iteration writes its own array elements.

cavity_melt_point plus the two EXCHANGE VELOCITIES the solve converged on.

cavity_melt_point_gamma with the Coriolis parameter passed as a per-column SCALAR — what cavity_melt_columns_2d calls, so no device code has to build a modified ocean_cavity_exchange_t (see cavity_exchange_velocities_f).

Is this step’s status tally a FAIL-LOUD condition?

Dimensional Obukhov length L = -u*^3/(kappa*B_b) (m), POSITIVE for a stabilising (melting) buoyancy flux. McPhee, Maykut & Morison (1987) p. 7029; the same scale appears as Yung et al. (2025) eq. (5) p. 5831. (Holland & Jenkins (1999) uses L_O in their eq. (18) but never defines it.)

Outer-iteration residual G(x) = ln(L+_new(x)) - x, with the destabilising branch folded in: a non-positive, non-finite or sentinel L+_new means the buoyancy flux at this iterate is destabilising (or zero), which every law treats as L+ = +infinity, so G = +infinity. Mapping it that way keeps the bisection bracket valid instead of taking log() of a negative number.

Resolve the &ocean_cavity_melt_nml gamma_s “unset” sentinel to the ISOMIP+ default gamma_t/CAVITY_GAMMA_RATIO_ISOMIP (Asay-Davis et al. (2016) Table 4 p. 2483; the 35 is Jenkins, Nicholls & Corr (2010) p. 2309).

The documented SAFE STATE returned on every non-OK status: zero melt, interface salinity equal to the far field, interface temperature on the liquidus there. m_mass is EXACTLY zero, so a failed column contributes nothing to a heat/salt budget rather than contributing a plausible wrong number.

The two sides of (E3), in (g/kg)*kg/m^2/s:

cavity_mask_impl with a constant multiplier folded in.

Define every component of the solution bundle. Called FIRST by cavity_solve_melt, before any early return, because the type carries no default initialisers (see its docstring: they would bar it from a do concurrent local(...) clause on gfortran).

cavity_solve_melt_f with the Coriolis parameter taken from the bundle (par%f_cor) — the scalar entry point every host caller and the kernel suite use.

Solve the three-equation system with any implemented exchange law, and return the full interface state plus the fluxes a coupling seam will consume.

Evaluate the whole interface at a trial x = ln(L+): exchange velocities at that stratification, the closed-form three-equation solve, the buoyancy flux the answer implies and the L+ it re-diagnoses. The fixed point of x -> ln(L+_new) is the solution of the implicit system.

Reduce the per-column status plane into four counts, ON DEVICE.

Ice temperature actually used: exactly zero when the mode ignores it, so an unset or stale T_ice cannot leak into an insulating run through L_eff.

Closed-form solve of (E1)-(E3) on the linear liquidus carried by eos. Returns the interface state and the canonical melt mass flux (kg/m^2/s, > 0 melting). See the derivation block above for the root selection, the cancellation-safe quadratic and the pre-solve melt/freeze branch.

The TRIMMED initial surface (MOM6 TRIM_IC_FOR_P_SURF, trim_for_ice): the free-surface anomaly eta_trim that puts a loaded column top exactly where the displaced water’s own weight equals the load, so the initial state is at rest under the MOM6 barotropic split (&ocean_bt_nml bc_pgf_forcing).

Two-equation variant: the interface salinity is the FAR-FIELD salinity, S_b = S_w exactly, i.e. the gamma_s -> infinity limit of the three-equation form (NOT the gamma_s -> 0 limit, which sends T_b -> T_w and the melt rate to zero). Holland & Jenkins (1999) section 2b(2) p. 1791; Jenkins, Nicholls & Corr (2010) eq. (6) p. 2302, where the single transfer coefficient is Gamma_TS ~ 0.006 and the freezing point is evaluated at the far-field salinity.

Interface friction velocity (m/s),

water = b - z_draft — the reference water-column thickness the datum, the layer split and the wet-mask seed all work on. Kept as one named routine so the three call sites cannot drift.

sum_k sum_ij f(i,j,k)*area(i,j) over the physical interior.

Per-cell advective CFL from the C-grid face-velocity pairs and metric inverses: cfl = (|u_c|·idx + |v_c|·idy)·dt, with u_c/v_c the face averages. !$acc routine seq so it inlines into the reduction loops below — same module ⇒ NVHPC keeps it inlined (and it’s a status-cadence cold path regardless). Sole home of the CFL formula, shared by the gated + un-gated compute_max_cfl loops.

Existence check for one tag’s file. NOT pure — inquire is an I/O statement, which is exactly why this is a subroutine setting has_error rather than another predicate: the pure ones above stay pure and independently testable, and this keeps the filesystem dependency in one visible place. Silent for a tag that names no file.

Ordering guard: update_2d/_3d must be called with the same t update_all most recently refreshed this field’s bracket with (skipped for STATIC — its bracket never changes, so “freshness” is meaningless). Catches a consumer calling update_2d/_3d before update_all has run this step, which would otherwise silently blend a stale bracket.

&ocean_isopycnal_nml check_h_positive guard. Abort on the FIRST negative layer thickness, naming the pipeline stage that produced it, the offending (i,j,k), and that column’s full thickness profile.

Fail-loud guard shared by every per-step/fill accessor: id must name an active field of the right rank, and the caller’s dest shape must match what was declared at registration.

&vcoord_nml remap_check_preconditions guard (audit findings V5, V6).

Fail loud on a cyclic group with no period. The reader makes the same check per field; catching it once here names the namelist group the user actually edited.

Flag unsupported namelist syntax embedded in a bare token.

Fail-loud TRIPWIRE for invariant I1′ — h_layer <= H_VANISHED ⇒ hTr = h_layer·c_live (the donor live layer’s concentration; hTr = 0 in a column with no live layer) for every registered tracer. Gated on &vcoord_nml check_vanished_content (default .false.), which is the knob the stability suite turns on for the cases that actually have vanishing layers.

Gate: enabled AND inside the step window.

Interior argmax of |arr|: first-encountered strict maximum over i ∈ [i0, i1], j ∈ [j0, j1], all k (the ghost ring is excluded by the caller’s bounds). Host-side (debug-window cadence only).

Sample the BT substep’s 2D in/out fields at the fold seam: the fold delta is bt_ubt_end - ubt_at_n - dt*F_bt_u, so these four (+ eta) name which INPUT went non-finite when the fold’s loud count fires — the “what fed it” the nan-catch cannot see.

Hot-face anatomy at a phase seam: interior argmax |u| and |v| with the local thickness pair (donor/receiver cells), the per-face viscous remnant, and the column context (thickness of the layer below/above at the max face). The forensic question this answers: WHICH face takes the explicit dt·F kick, is it an outcrop edge (massive|vanished thickness pair), and is visc_rem actually small there (i.e. would MOM6’s attenuation have caught it)? Row format (grep “HOTFACE”): HOTFACE s u|v (i,j,k) val hL hR rem h_dn h_up Debug-window only (chksum_active gate); D→H of the two face fields per sample — affordable at window cadence, not at production cadence.

Map a grid-location tag (LOC_{H,U,V,Q}) to the field’s 2D extents on the Arakawa-C grid — u/v faces carry the extra wall-normal row/column, the corner both. This is the single place the C-grid extent convention lives; call sites pass loc.

Physical-cell index range for a grid-location tag: the ghost ring excluded. Face locations carry one MORE physical entry than cell centres in their staggered direction (an x-face array spans nx_phys + 1 faces), which is why LOC_U/LOC_Q extend i1 by one and LOC_V/LOC_Q extend j1.

Emit one CHKSUM row (write(*,…) like the KE_ATTR probe — probe output bypasses the logger by design: greppable, no prefix, survives logger-level filtering).

Sample the prognostics (h, u, v, every tracer) at a phase seam. Call AFTER the phase named by label; drains device queues first so the async apply chain has landed. Bounds come from the location-aware API (LOC_H / LOC_U / LOC_V), not hand-written.

2D twin of chksum_stats_3d (BT work fields).

Device-side sum/min/max/nonfinite over an explicit-shape 3D field. No present clause: present_or_copyin reads the device copy in production and copies-in host data in unmapped unit tests (same convention as the BT fold’s loud count).

Clamp an open-area fraction into [0,1], NaN-safely.

Flush + close the NetCDF file. Idempotent.

Refresh the BAROTROPIC face widths from the LIVE layer thicknesses when &vcoord_nml zfixed_closed_faces is on.

Tokenize a group body up to and including the terminating ‘/’. Token kinds are encoded by their text: “=” , “,” , a quoted or bare value, an identifier, or special markers for unsupported syntax which are flagged here. The ‘/’ ends collection.

Abort all MPI processes with given error code

Broadcast a single real(wp) scalar from compute rank 0 Uses compute communicator so I/O server ranks don’t participate.

Return the compute-only communicator Falls back to comm_world() if comm_env_init has not been called (e.g. tests that use raw MPI_Init).

Return this process’s rank in the compute communicator

Return number of compute ranks (excludes I/O ranks)

Finalise MPI. No-op if the comm-env was never initialised, so it is safe to call unconditionally at program exit – the shared per-test main does this for the few tests that init MPI, and it costs nothing for the rest (and on serial builds via the stub).

Return the cached global (world) communicator.

Phase 1: Initialise MPI, cache world rank/size Call this before reading config. Idempotent: a second call is a no-op, so multiple entry points / testdrive cases can call it.

Return world rank of this node’s I/O server (-1 if none)

Return .true. if this rank is a dedicated I/O server

Return world ranks of compute processes on this node

Return number of compute ranks on this node

TEST/TOOLING-ONLY. Pops the override pushed by comm_env_push_compute_comm, restoring comm_env_compute_comm() to the real compute communicator. Fails loud on pop-without-push (a bug in the caller, not a state this seam should absorb silently). Does NOT finalize the popped communicator – the pusher owns that (e.g. freeing a self-comm built for the push).

TEST/TOOLING-ONLY. Temporarily overrides what comm_env_compute_comm() returns to comm, until the matching comm_env_pop_compute_comm().

Return this process’s world MPI rank

Phase 2: Assign roles (I/O server vs compute), bind GPU Call this after reading config, before decomp/solver init.

Return total number of MPI processes (including I/O ranks)

Shared tail of rdb_ocean_create_from_string AND rdb_ocean_create_finalize (P2.5): engine_setup -> resolved-n_inner check -> engine_enter_data -> finalize handle bookkeeping. Consumes whatever geometry h%engine%staged_* fields carry (empty/unset for the single-call path — byte-identical to before P2.5). On any failure, destroys the WHOLE handle (F9) and clears the single-live-handle guard; c_handle is c_null_ptr and h is null on return in that case, matching rdb_ocean_destroy’s out-null convention.

Per-face multiplicative damping factor for the BT-substep velocity update (linear-drag branch): bt_rem_face = Htot_face / (Htot_face + r·hbbl·dt_inner) applied as ubt_new = bt_rem_u·(ubt_old + dt_inner·forces) each inner step. When the bt_substep_drag knob is off this must NOT be called and the workspace stays at 1 (no-op, bit-identical).

PR-2 (bt-rem-from-av-rem): build bt_rem_u/v from the SAME viscous remnant the layered momentum solve uses, MOM6’s barotropic solver. Dispatched the same way as compute_bt_rem — a RESETTER, mutually exclusive at configure with bt_substep_drag (D2, double-counted bed drag) and with bt_halo > 0 (validate_config) — so this and compute_bt_rem/ reset_bt_rem never both run for the same stage; src/core/ ocean/README.md’s “exactly one resets, everything else MULTIPLIES” contract gets this as its third resetter.

MULTIPLIES the Egbert & Ray (2001) / Jayne & St Laurent (2001) linear (Rayleigh) barotropic wave drag into bt_rem_u/v: bt_rem_u *= Htot_face / (Htot_face + lwd_drag_u·dt_inner) lwd_drag_u/v is a static, face-resident piston velocity [m/s] built once at configure by configure_ocean_wave_drag. Uses the IDENTICAL Htot_face expression as compute_bt_rem (reuse, not a second H_tot). Composes with substep_drag exactly as MOM6 composes lin_drag_u with the viscous remnant. Htot_face <= 0 ⇒ leave bt_rem unmodified (MOM6’s guard).

Device kernel for the per-layer side-drag Rayleigh rate. See ocean_channel_drag_compute_tendencies for the derivation. wet_q is (nx+1, ny+1); dyCu is (nx+1, ny) (u-face length normal to the zonal flow); dxCv is (nx, ny+1).

HBBL-distributed bottom drag. Mirrors MOM6’s LINEAR_DRAG and quadratic-with-HBBL formulations: the drag stress is spread across the bottom hbbl metres rather than dumped into the bed-most layer.

SSH anomaly = 0.5·(bt_eta_end + bt_eta) − eta_PF: the part of η the BT substep produced beyond what the slow PGF saw. Zero at steady state.

Faithful port of MOM6’s PressureForce_FV_Bouss per-layer PGF for the Boussinesq + per-layer Rlay path.

FV_MOM6 pressure gradient, constant-by-layer (PCM) T/S, density at the IN-SITU pressure — MOM6 PressureForce_FV_Bouss with RECONSTRUCT_FOR_PRESSURE = False (int_density_dz_generic_pcm).

FV_MOM6 pressure-gradient with in-layer T/S reconstruction.

Reduced-gravity / gprime PGF for NK = 2.

Face-centred depth-weighted column averages of pbce (gtot_E/W/N/S). Wall cells fall back to pbce(:,:,nz). By construction Σ_k h_face(k)·(pbce(k) − gtot_face) = 0 per column, making the bc-PGF Δu correction depth-mean zero – UNDER THE SAME h_face/weight derive_bt_from_layers and apply_bt_correction’s folds use.

Per-face upstream column-sum thickness h_face_up_x/y(I,j) = Σ_k h_layer(I_upstream,j,k) used by the BT chain when use_upstream_h_face = .true.. First-order upwind pick by face-velocity sign (sampled at the top of the outer step). Wall faces use the single available cell. No-op when the knob is off.

Σ wet_T·areaT, Σ ci·areaT and Σ (mice/ICE_RHO_ICE)·areaT over PHYSICAL cells (ghosts excluded) — ci/mice from the two-mode per-cell gather (ice_cell_concentration_impl convention, inlined; test_ocean_ice_diags pins the fills’ copy of the same math). One pass, three reduction(+:) accumulators.

EFP twin of compute_ice_totals. 2D-only (no k-slab blocking needed – one “slab” per accumulator), so a single efp_summands_guard call suffices. THREE SEPARATE single-pass reductions (one per accumulator, 7 reduction scalars each) rather than one kernel combining all 21 – ice diagnostics are a status-cadence cold path (SS11.11: “the EFP path costs ~nz times more kernel launches… unmeasurable at status cadence”), so the extra category-sum pass for ci/hi is free.

max (|u_c|·dt·idxT + |v_c|·dt·idyT) over PHYSICAL cells (ghosts excluded). idxT/idyT are metric inverses (= 1/dx,1/dy on uniform).

Per-layer pressure-anomaly gravity coefficient (m/s²): the response of layer k’s pressure to a unit change in η. Montgomery form, bottom-up convention (k=1 bed, k=nz surface): pbce(:,:,nz) = g·ρ_ref/ρ_0 do k = nz-1, 1, -1 g_prime_K = g·(rho_layer(k+1) − rho_layer(k))/ρ_0 pbce(:,:,k) = pbce(:,:,k+1) + g_prime_K·(e_top_of_k − e_bed)/H Uniform-density column ⇒ pbce−gtot ≡ 0 ⇒ bc-PGF correction a no-op. Reads pgf%e_face; requires pgf%variant == OPGF_VARIANT_FV_MOM6 (other variants don’t fill e_face). validate_config and configure_ocean_pgf refuse correction_bc_pgf with any other form, so the error stop below is a backstop for direct callers.

Σ h_layer(i,j,k)·areaT(i,j) over PHYSICAL cells (ghosts excluded). Explicit OpenACC reduction — the sum() intrinsic on a present-mapped array silently runs host-side under NVHPC non-managed mode and returns the stale host shadow.

EFP twin of compute_total_h: order-invariant fixed-point Sigma h_layer*areaT over PHYSICAL cells. K-SLAB BLOCKED: a host loop over k, one device reduction(+:e1..e6) per slab, then a host-side efp_carry combining the slab into the running total – keeps each device reduction block within EFP_MAX_SUMMANDS (SS3.3/SS6.3 of the plan; MOM6’s i/j block-partition arithmetic is NOT ported – the k-slab is simpler and sufficient at EFP_PREC_WIDTH = 36). Ghost exclusion + extent clamping copied VERBATIM from compute_total_h – a divergence here would silently change what is summed between the FP and EFP paths.

Σ 0.5·h·(u_c²+v_c²)·areaT over PHYSICAL cells (ghosts excluded), using cell-centred face averages. Direct OpenACC reduction.

EFP twin of compute_total_ke. See compute_total_h_efp for the k-slab blocking design; face-averaging + extent clamping copied verbatim from compute_total_ke.

Σ hTr(i,j,k)·areaT(i,j) over PHYSICAL cells (ghosts excluded). Same explicit-reduction pattern as compute_total_h.

EFP twin of compute_total_tracer. See compute_total_h_efp for the k-slab blocking design; ghost exclusion + extent clamping copied verbatim from compute_total_tracer.

Fill ms%w_interface by integrating the horizontal- continuity residual upward from the bed. Eulerian z:

Set logger verbosity from string

Bake the C3 nodal/astronomical correction into one OBC edge, then fail loud (rank-0 error + error stop) if any of the edge’s constituent frequencies matches no tide-catalog entry, else log the resolved constituent → (name, f_c, V_c+u_c) map on rank 0. Host-side setup wrapper around the pure obc_tide_nodal_fill — the pure helper signals the failure via fill_ierr; the loud policy lives here.

Populate ocean_state%bc edge tags + per-edge data values from cfg%ocean%bc (read from &ocean_bc_nml). Run this after configure_ocean_bt_split and before ocean_state_enter_data so the BC state is set before the first dyn step.

Barotropic-substep correction knobs (MOM6 frhatu h-weighting, bc-PGF retro-correction, bt_rem_u drag damping, visc_rem joint weight), the BT_cont_type / upstream-PPM h_face workspace allocations, and the rank-0 PGF/BT configuration log lines.

Auto-derive the barotropic substep count n_inner from the external gravity-wave CFL (MOM6 set_dtbt) when requested, then latch the mode-split reference column depth bt_H_ref from the seeded bathymetry. cfg is intent(inout) because auto_n_inner writes cfg%ocean%bt%n_inner.

Build the static ice-shelf cavity LOAD field metrics%p_ice_ref = (rho_ref*GRAVITY)*z_draft (Pa), ASSEMBLE it into the top-of-column pressure multilayer_state_t%p_top, and assert the counted-once datum invariant.

Configure the ice-shelf basal-melt slot (&ocean_cavity_melt_nml, P2b): copy the knobs onto the slot’s three flat parameter bundles, resolve the gamma_s sentinel, build the per-column Coriolis array the hj99 law needs, and CHECK that the interface pressure the liquidus will read has actually been loaded.

Build the static partial-step z-level FACE-CLOSURE mask (&vcoord_nml zfixed_closed_faces; Adcroft, Hill & Marshall 1997; Losch 2008 §2.1 for the ice-shelf cavity).

Copy the &ocean_conv_nml knobs onto ocean_state%vmix – every field, no dead-config gaps. Convective adjustment adds no new slot / allocatable (it lives on the already-unconditionally- allocated vmix), so unlike configure_ocean_tidal_mixing there is no separate enable latch to set on a distinct sub-object; vmix%conv_enable IS the latch. No mutual exclusion with KPP / EPBL: convection masks against whichever BL depth is live this stage (rdb_ocean_dyn.F90 vmix_apply_in_stage).

Copy the &ocean_ddiff_nml knobs onto ocean_state%vmix – every field, no dead-config gaps. Like convection, double diffusion adds no new slot (it rides the unconditionally-allocated vmix and the already-mirrored vmix%eos); vmix%ddiff_enable IS the latch. Folded into vmix_split_kd_heat_salt, so it needs the interior closure pipeline (use_closure) and thermodynamics (it reads T/S and the EOS alpha/beta). enable=.false. => bit-identical.

Bottom-drag variant + coefficients, the continuity PPM positivity guard, and the BT-budget diagnostic probe — plus their rank-0 log lines. All knobs default off/zero (bit-identical to pre-knob nmls).

Copy the &ocean_epbl_nml knobs onto the EPBL slot — every field, end to end (don’t repeat the KPP ri_crit/c_vt2 dead-config gap). Also fills f_centre from the same beta-plane parameters the Coriolis slot uses, copies the EOS hookup, validates the configuration, and resolves the EPBL-vs-KPP mutual exclusion.

Wire surface wind stress, horizontal-viscosity coefficients, and the Coriolis beta-plane + PV-scheme variant from cfg into the ocean slots.

Copy the &ocean_foxkemper_nml knobs onto the MLE slot (B5). Validates: B5 reads epbl%mld, so EPBL must be enabled; and the resolution_taper hook is a hard error until B2 lands.

Along-coordinate tracer Laplacian coefficient (&ocean_hdiff_nml kappa_h). Scalar copy onto ocean_state%hdiff_tracer, mapped with its parent slot at ocean_state_enter_data — no explicit !$acc update needed as long as this runs before that (it does; see rdb_driver.F90). Default kappa_h = 0.0 leaves the kernel’s short-circuit intact ⇒ bit-identical.

Fill ms%k_bot / k_bot_u / k_bot_v — the shared index of the first LIVE layer counting UP from the bed, and the field every bed-side consumer reads instead of spelling 1. The bed-side mirror of configure_ocean_k_top.

Fill ms%k_top / k_top_u / k_top_v — the shared index of the first LIVE layer counting down from the top, and the field every top-side consumer reads instead of spelling nz.

Copy the &ocean_kappa_shear_nml knobs onto the kappa-shear slot — every field, end to end (don’t repeat the KPP dead-config gap). Fills f_centre from the same beta-plane parameters the Coriolis slot uses, copies the EOS hookup, and validates. No mutual exclusion: kappa-shear is an interior closure that coexists with KPP / EPBL and PP81/background.

Derive the static C-grid land masks from the seeded T-cell wet_mask and zero the 6 face metrics at land faces (metrics_apply_land_mask). Run AFTER configure_ocean_metrics (the metrics + inverses must exist) AND configure_ocean_bc (the periodic / fold flags drive the halo-aware mask derivation), but BEFORE ocean_state_enter_data (the host edit is what the GPU copyin captures).

Flow-aware lateral-viscosity closure (Leith / Smagorinsky + biharmonic Smagorinsky_AH), the vertical coordinate (VCOORD_* code + z_fixed reference depth), and the PP81/KPP vertical-mixing switches — with their rank-0 log lines.

Fill the MEKE slot’s cell-centre |Coriolis| from the same metrics_fill_coriolis path VarMix / EPBL use (handles beta-plane AND spherical), so beta = |grad f| for the Rhines length is live when alpha_rhines > 0. The MEKE scalar knobs are copied earlier by ocean_state_copy_config; this only fills f_centre. Host loop before enter_data. No-op when MEKE is disabled.

Fill the ocean_metrics_t slot per cfg%ocean%grid%grid_config, then single-source the inverses + hvisc ratio bundle (metrics_finalize). Must run BEFORE ocean_state_enter_data (the host fill is what the GPU copyin captures). For “spherical”, the &grid_nml dx/dy are reinterpreted as dlon/dlat in degrees.

Configure the atmospheric surface-pressure loading slot (PR-17): copy enable, take ρ₀ from ocean_state%eos%rho0 (the single ρ₀ of record — NOT a namelist knob), and allocate the seam fields. Host-side setup — runs BEFORE ocean_state_enter_data. enable=.false. => no-op, bit-identical. The uniform-p_surf inert warning is emitted in validate_config; here we only log the enable on rank 0.

Pressure-force variant, the reference densities (rho0 / rho_ref, both from the single configured ρ₀ — &ocean_ic_nml rho_0 via eos%rho0), the reduced-gravity (gprime / gfs_scale) knobs, the bathymetry copy into the PGF slot, and the matching barotropic fast-loop gravity g_bt for the gprime / FV_MOM6-reduced-GFS paths.

Configure porous barriers (&ocean_porous_nml, Adcroft 2013). Grows the ocean_metrics_t porous arrays to full face size and fills the STATIC along-face d_min/d_max/d_avg statistics; the layer-averaged open fractions themselves are recomputed on the device every RK2 stage (they depend on the interface heights).

Fan the ONE configured Boussinesq reference density out to every remaining slot that carries its own rho0 copy.

Populate ocean_state%sponge’s per-cell idamp_h/idamp_u/ idamp_v maps from cfg%ocean%sponge (&ocean_sponge_nml) + the already-configured ocean_state%bc edge tags. Run AFTER configure_ocean_bc (reads the edge tags + has_* flags) and BEFORE ocean_state_enter_data.

Copy the &ocean_tidal_mixing_nml knobs onto the tidal-mixing slot — every field, end to end (no dead-config gaps). Seeds the prescribed bottom energy field e_in (v1 uniform path), copies the shared EOS hookup for the N^2 buoyancy derivatives, and validates. No mutual exclusion: tidal mixing is an interior closure that coexists with KPP / EPBL / PP81 / kappa-shear.

Configure the equilibrium body-force tide slot (C1): parse the constituent list + reference dates, fill the astronomy catalog (phase0, nodal f/u), and build the (nx,ny,3) spatial-structure arrays from metrics%geolatT/geolonT. Host-side setup — runs AFTER configure_ocean_metrics (lat/lon must be filled) and BEFORE ocean_state_enter_data. enable=.false. => no-op, bit-identical.

Ice-shelf TOP drag (&ocean_tdrag_nml, Phase 4a): copy the variant + coefficients onto the slot, project the static cell-centred metrics%cover_frac onto the velocity FACES, and enforce the ONE-C_d rule against the melt slot.

&ocean_tracers_nml scalar knobs (PR-7): the ideal-age Dirichlet surface value and its vintage-mode exponential growth rate. enable_ideal_age itself is latched earlier by ocean_state_copy_config (before init(grid), since it gates tracer-slot allocation) — these two scalars gate nothing, so they land here in the post-init configure pass. Both default to 0 ⇒ bit-identical.

Build the STATIC VarMix grid terms (f2_dx2_*, beta_dx2_*, l2_*) on the varmix slot once at configure time, from the filled curvilinear metrics + the cell-centre Coriolis magnitude. The slot scalars are copied earlier by ocean_state_copy_config; this only fills the static arrays (the per-step Res_fn/SN/assembly fire in the dyn step). No-op when VarMix is disabled.

Thermodynamics on/off, velocity-truncation clamp (MAXVEL), DIRECT_STRESS surface-stress distribution, KV_ML_INVZ2 surface-band viscosity, HARMONIC_VISC face-thickness mean, and the DT_THERM thermo/tracer cadence — with their rank-0 log lines.

Configure the barotropic linear (Rayleigh) wave-drag piston-velocity maps (Egbert & Ray 2001; Jayne & St Laurent 2001) — the bulk energy sink for the barotropic tide, MOM6 BT_LINEAR_WAVE_DRAG. Builds a host-only h-point r_h(nx,ny) map (uniform scalar or a resolved-bathymetry roughness proxy), scales it, averages h->face into bt_work%lwd_drag_u/v, and leaves the arrays unallocated when the knob is off (bit-identical). MUST run AFTER bathymetry is set (ocean_state%barotropic%b) and land masking (configure_ocean_land_mask) and BEFORE ocean_state_enter_data — the !$acc enter data copyin in barotropic_workstate_enter_data carries these host-filled values to the device (CLAUDE.md gotcha (2): arrays mapped create do not carry pre-map host values, so this ordering is load-bearing).

Copy the &ocean_wavespeed_nml knobs onto the wave-speed slot (B1). Diagnostic, no mutual exclusion: cg1/Rd read rho_layer directly. Fills f_centre + the static beta_centre = |grad f| field via fill_coriolis_centre + build_static — the SAME metrics_fill_coriolis path VarMix/MEKE use (planetary on spherical/tripolar, beta-plane bit-identical elsewhere) — rather than the legacy hard-coded beta-plane set_f_centre. Copies rho0 from the EOS slot for the Boussinesq gprime. Must run AFTER configure_ocean_metrics (metrics filled + finalized).

Dynamic wetting/drying (docs/ocean_wetdry_plan.md): copy the &ocean_wetdry_nml knobs onto the BT workstate, allocate the wd_* workspaces (lazy — absent when the knob is off, so the default path carries no new arrays and stays byte-identical), and seed the hysteresis wet mask from the seeded bathymetry (barotropic%b — the same array configure_ocean_bt_split later latches into bt_H_ref; bt_eta is still 0 here). Must run BEFORE the restart read (so wd_wet_dyn is allocated + registered when the registry walk runs and a warm restart overwrites the seed with the saved front state) and BEFORE ocean_state_enter_data (host seeding; the enter_data walk attaches whatever is allocated).

Resolve the VCOORD_Z_FIXED nominal layering onto the vcoord slot — and the VCOORD_ZSTAR one, which is the same nominal profile (MOM6 z* dilates it per column; ocean_vcoord_zstar_target): z_fixed_h_ref = &ocean_topo_nml max_depth (the uniform max_depth/nz spacing — the default, byte-identical), or, under &vcoord_nml z_fixed_profile = "list" | "tanh", the stretched per-layer tables z_fixed_zi / z_fixed_dz built by rdb_vcoord :: z_fixed_nominal_dz (z_fixed_h_ref then becomes the profile’s total depth).

Populate the isopycnal rho_target(0:nz_ml) interface densities for VCOORD_RHO / VCOORD_HYCOM. rho_target(0) is the surface (lightest) interface, rho_target(nz_ml) the bed (densest).

Emit one MOM6-style console block from pre-computed totals. The caller has already area-weighted + reduced each scalar over its own state; this routine only latches the t=0 reference (first call), formats the lines, and runs the panic guards. When a budget with active=.true. is passed, the Error column becomes the boundary-flux + surface-source-corrected residual (see conservation_budget_t) and the outflux / source are shown.

Test-only (no production caller): unsplit apply, paired with continuity_compute_fluxes as the split path’s reference oracle. Per-layer forward-Euler thickness update.

Forward-Euler step: h <- h - dt * flux_h. The full split-explicit RK2 scheme (Phase 4) wraps two of these calls around an RK2 averaging pass; for Phase 2 this single-stage step is enough to exercise the kernel under the lake-at-rest, Gaussian-hump, and mass-conservation tests.

Apply the meridional (y-flux) thickness update on top of the zonally-updated state: h(i,j,k) ← h(i,j,k) - dt · (Φy(i,j+1,k) - Φy(i,j,k)) · iareaT Adds the y-divergence to flux_h_layer so the field ends the split step holding the total horizontal divergence that the vertical-advection kernel consumes (w_interface(k+1) = w(k) - flux_h_layer(k)). Φy carries dx_cv; iareaT = inv_dy on uniform metrics.

Apply the zonal (x-flux) thickness update: h(i,j,k) ← h(i,j,k) - dt · (Φx(i+1,j,k) - Φx(i,j,k)) · iareaT Overwrites flux_h_layer with the x-divergence so the companion meridional apply can accumulate the total. Φx is the width-weighted transport (m³/s); iareaT closes the divergence to a per-area rate (= inv_dx on uniform metrics).

Counted allocatable footprint of the continuity-PPM slot (0 when unallocated). One arr_bytes term per array — add a term here when a new allocatable joins the type.

Test-only (no production caller): the unsplit reference path, kept as the oracle the split production path is checked against. Multilayer counterpart to continuity_compute_fluxes_barotropic: identical PPM reconstruction + upwind face pick + flux divergence, lifted per-layer. Each k-slice is independent (the PPM stencil reads only the same k), so the do-concurrent kernels parallelize over (k, j, i) simultaneously for GPU occupancy.

PPM face reconstruction + per-face mass flux + cell-centred flux divergence for the barotropic C-grid state.

Bare copyin(this) removed (stack-descriptor map → AMD cross-slot overlap; see ocean_surfstress_enter_data). The face buffers attach below; ct-descriptor presence (so DCs touching ct%h_face_left_x%data don’t per-launch memcpy) comes from the root copyin(state) in ocean_state_enter_data. A V100 A/B with copyin(this) gone is bit-identical and faster overall, so the root copy fully covers it.

Gent-McWilliams thickness diffusion as its OWN sequential operator: move h_layer AND every tracer by the bolus transport gm%uhD/ gm%vhD, which gm_compute_transports has JUST filled from this same, untouched h_layer with this same dt.

Allocate the 4 face-reconstruction scratch buffers sized at (nx_face, ny_face, nz). Default nz=1 covers the barotropic kernel; passing nz_ml sizes them for the multilayer kernel without forcing a separate init routine. Ocean init passes state%multilayer%nz_ml when the multilayer state is in play.

Meridional (y-only) PPM reconstruction + per-face mass flux. Mirror of continuity_zonal_flux, with the same optional vhbt transport-constraint renormalisation. Writes ms%mass_flux_y_layer. Walls at j=1 and j=ny+1 zeroed. In the Lie split this runs after the zonal apply, so it reconstructs against the already-updated ms%h_layer.

Test-only (no production caller): continuity-only split wrapper; production runs the tracer-interleaved continuity_tracer_step_split. Directionally-split (Lie) PPM continuity step over dt:

Phase-2 (6b) windowed horizontal tracer-advection drain.

Production entry point for the directionally-split continuity + tracer step. Interleaves the two so the CWC discrete theorem holds in the split form:

Zonal (x-only) PPM reconstruction + per-face mass flux on the multilayer C-grid. Companion to continuity_meridional_flux for the directionally-split (Lie) continuity step. Writes ms%mass_flux_x_layer and leaves mass_flux_y_layer / flux_h_layer untouched. Wall faces at i=1 and i=nx+1 are zeroed (closed-wall BC).

Return the default remap procedure pointer for an output vgrid. Used by the registration path to attach the right conservative remap when a diagnostic selects a non-layer output coordinate. Null for LAYER / unknown (no remap needed).

Shared device copy buf = src with shape clipping — every direct-copy derived field routes through here.

src -> dst flat copy on the device. Bare-array shim avoids the deep struct deref inside do-concurrent (tracer hTr lives in an array-of-derived-types registry).

Device pass: per-face PV/∇KE work split, three regions. Faces attributed to their west/south cell’s region mask.

Per-layer forward-Euler velocity update. no_wait (optional, default .false.): when .true. the apply DC loops are issued on OpenACC queue 1 and the routine returns WITHOUT syncing, so a batched caller (run_stage_split velocity-apply chain) can pipeline the whole additive apply sequence and !$acc wait(1) ONCE. Default ⇒ self-contained blocking apply (historical, safe for non-batched callers — e.g. the unsplit run_stage). Not pure because of the async/wait directives; still functionally pure.

Forward-Euler velocity update from the tendencies the compute step wrote into pv_flux_x / pv_flux_y. u_face_x(i, j) <- u_face_x(i, j) + dt * pv_flux_x(i, j) v_face_y(i, j) <- v_face_y(i, j) + dt * pv_flux_y(i, j) Split-explicit RK2 (Phase 4) wraps a pair of these around an RK2 averaging pass.

Counted allocatable footprint of the Coriolis-advection slot (0 when unallocated). One arr_bytes term per array — add a term here when a new allocatable joins the type.

Per-layer Coriolis + advection dispatcher. Reads this%pv_variant and routes to the matching kernel body:

Sadourny (1975) energy-conserving Coriolis + horizontal- momentum-advection form on the barotropic C-grid state:

Public only for the unit-test suite (no production module imports it); ignore when developing production code in other modules. Per-layer PV-conserving Coriolis + horizontal-advection tendency in the Arakawa-Hsu (1990) form (“HK correction”). The wider 3-corner PV stencil at each face suppresses the spurious Hollingsworth-Källén instability that biases the simpler Sadourny 2-corner form at eddy-resolving resolutions.

Per-layer Sadourny Coriolis + advection tendency. Same algorithm as the barotropic counterpart, lifted with a k-axis on every loop. Each k-slice is independent (ζ stencil only reads same-k velocities; KE at centre only reads same-k face values), so the do-concurrent kernels parallelise over (k, j, i) for full GPU occupancy.

Faithful MOM6 SADOURNY75_ENERGY (Sadourny 1975 energy-conserving) per-layer Coriolis + horizontal-advection tendency. This is the TRANSPORT form: the absolute-vorticity flux is the potential vorticity q = (f + ζ)/h_at_corner times the layer MASS TRANSPORT (vh/uh), so the discrete Coriolis term produces zero net domain kinetic energy (energy-conserving). The default enstrophy form (_sadourny, (f+ζ)·v) only matches this under uniform thickness.

Allocate the 4 scratch buffers sized at (nx_face / ny_face / corner, nz). Default nz=1 covers the barotropic kernel; passing nz_ml sizes them for the multilayer kernel. Same backward-compatible pattern as continuity_init.

Populate f_corner with a beta-plane profile f(y) = f_0 + beta * (y - y_ref)

BOUND_CORIOLIS abs_vort recovery: return (f+ζ) at corner (ic,jc) by multiplying the PV q_val back by the corner thickness Pass 2 divided abs_vort by — recovering (f+ζ) to round-off. Recomputes the SAME wet-area-weighted 4-cell hm_num/hm_den, then reconstructs the effective corner thickness with the SAME formula (and floors) the active corner_h variant used in Pass 2: cell_mean (default): h_corner = max(hm_num/max(hm_den,H_DIV_EPS), CORIOLIS_H_MIN_PV); av = q·h_corner. mom6_area (use_mom6_ch): Pass 2 formed q = abs_vort·hm_den/ (hm_num + PV_VOL_NEGLECT), so the consistent inverse is av = q·(hm_num + PV_VOL_NEGLECT)/hm_den (hm_den guarded by H_DIV_EPS for the fully-land corner, where q≡0 ⇒ av≡0 anyway). Without matching the variant the recovered abs_vort would be slightly inconsistent with how q was made whenever BOTH bound_coriolis and corner_h=”mom6_area” are on — visible only in the truly-vanishing- thickness limit. Chosen over a persistent abs_vort buffer so the default-off knob costs ZERO memory.

.true. iff all four cells contributing to a corner sample are wet. The masks are real 0/1 (ocean_metrics_t%wet_T), so the test is a mid-point comparison rather than an equality.

Horner cos Taylor on |r| ≤ π/4. 9 even-power terms.

Self-contained, non-erroring dimension validator — the single-rank (i_offset_global = j_offset_global = 0) testable twin of register_common’s inline checks; mirrors zinit_dims_ok. Opens filename, inspects var’s rank and lengths, and returns .false. on any mismatch, missing variable, or I/O error — NEVER aborts, so test code can probe the false branch without a subprocess/death-test harness (the house convention — see zinit_dims_ok).

Bracket search + blend weight for a query time t (already in file-time units/offset — the caller applies t_offset/t_scale before calling). nt == 1 is degenerate: always returns n0 = n1 = 1, w = 0, never out of range.

Host-ONLY NetCDF slab read for file record rec -> fld%f0 (default) or fld%f1 (into_f1 = .true.), applying scale/add_offset once at read. Uses the module-level host workspace (registration/bracket-advance-time only — never a per-step allocation). Deliberately carries NO !$acc directive: called both before enter_data (STATIC field, at registration) and after it (LINEAR/CYCLIC bracket advance) — the caller pushes to the device itself, only when that is actually correct (data_input_refresh_brackets).

Host-side registry walk (outer shim) for field id.

Translate a namelist/registration-time string into a DATA_TIME_* code. Fail-loud on anything not covered by data_input_time_mode_is_implemented — an unrecognised tag must never silently fall back to a default mode.

.true. iff tag (case-insensitive) names a shipped time mode. Drives the fail-loud dispatch in data_input_time_mode_from_string — house idiom, see lateral_closure_is_implemented.

CF units attribute (“seconds since …”, “hours since …”, …) -> a multiplier converting the raw file time axis to seconds. Only the leading unit word matters (no calendar, no reference date — see module docstring). ok = .false. for an unrecognised leading word; caller decides the fallback.

.true. when at least one tag names a file, i.e. enabling the group would actually do something.

.true. unless a tag names a file without naming the variable inside it. A blank file (tag not file-driven) is valid, and a var with no file is harmless — the tag simply never registers — so only the one combination is rejected.

.true. when evap/lprec are file-driven but the surface-flux component set they write into is switched off. Note the polarity: this reports the PROBLEM, not validity — named for how it reads at the call site.

.true. unless time_mode='cyclic' was requested without a positive cycle_period. Checked only for the cyclic mode; cycle_period is ignored by linear/static.

Days since the astronomical origin 1900-01-01 00:00 UT.

Choose px, py to minimise halo communication cost

Public only for the unit-test suite (no production module imports it); ignore when developing production code in other modules. Convert global physical indices to local physical indices

Initialise decomposition for a given rank

Initialise decomposition from config, with optional auto-factoring

Public only for the unit-test suite (no production module imports it); ignore when developing production code in other modules. Convert local physical indices to global physical indices

Log the process grid + this rank’s subdomain shape. Call from rank 0 after decomp_init_from_config so a run’s decomposition is visible at startup (the halo perimeter is what auto_factor minimises, printed here as the per-rank interior:ghost ratio cue).

Compute rank from process grid coordinates (row-major)

Populate bt_eta, bt_ubt, bt_vbt from the current multilayer state. bt_H_ref must already be set. bt_eta = Σ_k h_layer − H_ref; bt_ubt = Σ_k(u·h_face)/Σ_k h_face. Face thickness averages the two abutting columns (wall faces use the single cell). With use_upstream_h_face, the interior face-h is the first-order upwind pick — consistent with compute_h_face_upstream.

Public only for the unit-test suite.

The DERIVED_REQ_* code a catalog entry is gated on, -1 for an unknown name. The DECISION register_derived fails loud on, exposed as a lookup so the suite can assert it without provoking error stop — the repo’s standing pattern for testing a fail-loud rule.

Public only for the unit-test suite.

Fail loud on a derived diagnostic whose prerequisite knob is off. Same three-step shape register_derived uses for an unknown name — error ring, logger, error stop — so the C ABI and the console both see it.

Sign of det(M(lam)) on rows 2..kc (Sturm/Hallberg recursion).

Active OpenACC device number, or -1 on CPU builds. Kept private — only the budget logger needs it.

.true. iff a density-coordinate diagnostic selection (global vgrid='density', or a per-diagnostic :density/:rho attribute anywhere in diags) has a usable rho_levels axis: at least one entry, strictly increasing (invert_density_targets assumes a monotone light->dense target list), and within the declared MAX_OCEAN_DIAG_Z_LEVELS array bound. Unlike sigma/z*, density bins have no auto-fill, so an unset axis must abort rather than silently size to zero. When density is not requested at all this is unconditionally .true. (no constraint).

The [diag] console line’s min / max / mean for one diagnostic buffer, over the finite cells only.

Translate a raw nc_check-style status (0 = ok) from one of the nc_* calls in open_stream into the caller’s ierr contract: .true. on success; on failure, .false. with ierr = OCEAN_STATUS_ERR_IO when ierr is present (closing ncid first, when given, so a mid-create failure does not leak the file handle), or error stops with the same generic text nc_check itself would use had ierr never been threaded through — byte-identical legacy behaviour when ierr is omitted. Mirrors rdb_bathymetry’s bathy_io_ok (P0.1 F1).

Index-based bounding box (test-only). Weight = 1 inside the closed interval [i0,i1] × [j0,j1], 0 outside. Indices 1-based and inclusive; ghost cells count if inside the box.

Counted allocatable footprint of one region mask (0 when unallocated). weight is device-mapped by ocean_diag_enter_data_impl, so it must carry a term — see diag_var_bytes, which folds this in per masked diagnostic.

Whole-grid mask — all weights = 1. Equivalent to “no mask”; exists for explicit registration and global-integral tests.

Horizontal section (test-only): one cell row at j = j_row, spanning i ∈ [i0,i1]. Integrates transport across a constant-y line (v_face_y_layer at j_row = meridional throughflow).

Vertical section (test-only): one cell column at i = i_col, spanning j ∈ [j0,j1]. Mirror of h_section for zonal throughflow.

Query the vanished-masking mode (used by the registration path to tag non-layer diagnostics with has_missing for the NetCDF writer).

Whole-array min / max / sum of a diagnostic buffer in ONE pass, over the FINITE cells only.

Counted allocatable footprint of ONE registered diagnostic (0 for every buffer that is unallocated).

Map the &ocean_diag_nml vgrid string to a DIAG_VGRID_* tag. Unrecognised (schema-validated upstream) falls back to LAYER. Moved from rdb_driver (P2.4) so engine_configure_diag can use it without a driver dependency; re-exported by rdb_driver for any existing caller of that name.

Map the &ocean_diag_nml output_precision string onto a NetCDF element type for the diagnostic DATA variables.

Build + solve the tridiagonal system per face column. u_face is either u_face_x_layer (x_face = .true., shape (nx+1, ny)) or v_face_y_layer (x_face = .false., shape (nx, ny+1)). h_layer is cell-centred (nx, ny, nz) and we average across the face direction to get the face thickness. kv_centre is the same cell-centred diffusivity field used by the tracer kernel — we average it across the face to get a face-located value at each interface.

Low-level geometry check used inline by register_common: .true. iff the variable has the expected rank and its two validated-length dims are each at least as long as required (decomposition-safe: a subdomain slab only needs offset+extent to fit, not to equal the file’s global length). Never aborts.

Resolve a group name to known / external / unknown and parse (or skip) its body up to the terminating ‘/’.

True iff z_draft is the FULL ghosted (nx_total, ny_total) array the seeders index with (ng+i, ng+j) — i.e. it matches h_layer’s horizontal extent. Guards the (1, 1) placeholder ocean_metrics_t allocates when the cavity is off from ever being read as a field.

Conservative upfront availability limiter on the accumulated window transports uhtr/vhtr, applied BEFORE drain_reconstruct_hprev.

Scale each interior x-face transport by its INFLOW-receiving cell’s factor (drain_avail_limit step 2). Face i between cell (i-1) and cell (i): uhtr(i)>0 ⇒ receiver i, uhtr(i)<0 ⇒ receiver i-1.

Meridional analogue of drain_avail_scale_x. Face j between cell (j-1) and cell (j): vhtr(j)>0 ⇒ receiver j, vhtr(j)<0 ⇒ receiver j-1.

dst = src (explicit-shape device copy).

Concentration field Tr = hTr / max(hprev, DRAIN_MIN_H) for the WENO drain (the CW path builds the same field as the first loop of drain_parabola_*; factored out so the WENO path can reuse it without the parabola coefficients).

MOM6 hup/hlos/min_h two-test limiter on the zonal face transport (volume units). Face i between cell (i-1) and cell (i). Positive flow (uhr_x(i) > 0), donor = cell (i-1): hup = areaT(i-1)·hprev(i-1) − areaT(i-1)·min_h hlos = max(0, −uhr_x(i-1)) (already-committed outflow via the donor’s OTHER (west) face) cap when (hup−hlos)−uhr < 0 AND 0.5·hup−uhr < 0. Negative flow mirror, donor = cell (i).

Meridional analogue of drain_limit_x. Face j between cell (i,j-1) and cell (i,j); positive donor = cell (i,j-1).

Rebuild the per-cell zonal CW PPM parabola from the CURRENT Tr = hTr/hprev (V2 — per pass). Interior cells (3..nx-2) use the limited PPM edges; the 2-cell boundary band falls back to PCM (aL=aR=Tr ⇒ swept reduces to the donor value, 1st order), matching tracer_advect_zonal_one_impl’s near-wall band. TODO(MOM6-fidelity): MOM6 advect_tracer keeps full PPM up to the wall (dropping to PCM only at genuine local extrema / zero mask2dCu faces), so we are 1st-order in the 2 cells nearest a true WALL where MOM6 is PPM-with-mask (periodic seams are fine — the wrap restores full PPM). Tracked divergence; revisit if near-wall tracer diffusion matters. a6 = 6·Tr − 3·(aL+aR). Mirror-T at land neighbours (C2); bit-identical for all-wet (wet_T≡1).

Meridional analogue of drain_parabola_x. aL = south-edge, aR = north-edge value of each cell. Mirror-T at land neighbours (C2); bit-identical for all-wet.

hprev = max(0, areaT·h_end + div(uhtr,vhtr)) · iareaT, then the vanishing-layer hatch hprev += max(0, 1e-13·hprev − h_end) (Adcroft & Hallberg 2006; reuse of VANISHING_LAYER_TOL thinking).

Re-weight a tracer’s thickness-weighted content onto a new layer thickness, holding the CONCENTRATION fixed:

drain_rescale_hTr + the closed-budget fill. The concentration hold / un-hold pair is NOT content-conserving cell by cell – hTr := Tr·h moves Σ areaT·hTr by Σ areaT·Tr·δh, which is only zero when Tr is uniform – so both halves have to be recorded or the closed budget is only valid on window boundaries and wobbles at every mid-window report. Recording both makes them cancel exactly, since the un-hold is the arithmetic inverse of the accumulated hold.

fld = fld - sub (uhr -= uhh).

MOM6 swept-average CW parabola flux for the zonal faces. Face i, donor = cell (i-1) if uhh>0 else cell (i). Per-pass Courant CFL = |uhh| / (areaT·hprev) on the donor, clamped [0,1]. uhh >= 0: F = uhh·( aR − 0.5·CFL·((aR−aL) − a6·(1 − ⅔·CFL)) ) uhh < 0: F = uhh·( aL + 0.5·CFL·((aR−aL) + a6·(1 − ⅔·CFL)) )

WENO analogue of drain_swept_flux_x: swept-average zonal face flux from the concentration field tr using the WENO rung ladder. recon is TRACER_RECON_WENO5/7/9 (1/2/3); the internal rung_max is recon + 1 (WENO5→2, WENO7→3, WENO9→4).

Meridional analogue of drain_swept_flux_x.

Meridional analogue of drain_swept_flux_x_weno.

hprev(i) -= (uhh(i+1) - uhh(i))·iareaT (volume div → thickness).

Meridional analogue of drain_update_h_x.

hTr(i) -= (F(i+1) - F(i))·iareaT (zonal flux divergence; dt is already baked into uhh⊂uhtr, so no dt here).

drain_update_tracer_x + the closed-budget fill: the SAME increment written to hTr is accumulated (times w) into budget_adv, so the console out term sees the horizontal tracer transport the windowed drain performs. Without this the drain moves tracer that ms%*_budget_horiz_adv never records, and the Heat/Salt Error column has to fall back to raw drift (which then reports a live surface flux as a “leak”). w is DRAIN_BUDGET_POST_AVERAGE_WEIGHT for every drain call site.

Meridional analogue of drain_update_tracer_x.

Meridional analogue of drain_update_tracer_x_budget.

Periodic wrap (+ north fold) of a cell-centred drain field. no_wait (optional): when .true. AND not folding, the periodic wrap is issued async on queue 1 without syncing, so a caller can batch several independent wraps (e.g. the pal/par/pa6 parabola triple) and !$acc wait(1) once. Ignored when fold_n (the fold reads the wrapped field, so the periodic wrap must complete first).

Periodic wrap (+ north fold) of an x-face drain field (nx+1,ny,nz).

Periodic wrap (+ north fold) of a y-face drain field (nx,ny+1,nz).

fld = 0.

Compute-rank entry point. Runs the Arakawa C-grid + continuity-PPM split-RK2 ocean dynamical core in src/core/ocean/.

Run the full compute-rank lifecycle for the ocean regime.

rdb --validate-only: everything driver_run does before the first step that can REFUSE a configuration — the full engine_setup sequence, so every fail-loud configure-time check in the configure_ocean_* stages (rdb_ocean_setup.F90), the stability audit and the IC seed — and nothing after it: no device mapping, no stepping, no output files (the diag stream is not opened; see engine_setup’s validate_only).

Bin 1 is NOT bounded by efp_carry (only bins 2..6 are – see the module docstring), so the double-precision transport needs its own guard: after summing nranks local values, the partial sum in bin 1 must stay <= 2**53 / nranks for the MPI_SUM-on-doubles combine to remain exact (the analogue of MOM6’s prec_error = huge(1_int64) / num_PEs, but for the int64-as- double transport rather than int64 transport). .false. ⇒ the caller must error stop, never silently proceed.

Renormalise bins 6..2 into (-2**P, 2**P), propagating the excess into the next-more-significant bin, without changing the represented value. Bin 1 is left untouched (unbounded by construction; see the module docstring on EFP_MAX_RANKS). Mirrors MOM6 carry_overflow, which loops EFP_DIGITS..2 for the same reason.

Greedy sign-magnitude fixed-point decomposition of r into six int64 bins of weight pr(n), n = 1..6. Unrolled (no loop, no array indexing over pr/I_pr) so this is a template a !$acc routine seq in-module duplicate can mirror exactly (see rdb_ocean_console_stats::efp_decompose_impl, which pins against this procedure in test_efp_impl_matches_canonical).

In-module !$acc routine seq duplicate of rdb_efp::efp_decompose – six SCALAR bin outputs (not an int64(6) array) so the call sites below can accumulate directly into reduction(+:e1..e6,epoison) clauses (OpenACC has no portable array reduction). Duplicated rather than called from rdb_efp because NVHPC’s device codegen does not inline a pure !$acc routine seq helper across a module boundary (CLAUDE.md Gotchas); test_efp_impl_matches_canonical (RDB_ENABLE_TESTING- gated) pins this copy bin-for-bit against the canonical procedure over the same magnitude table – compute_ice_totals’s docstring documents the identical pattern for ice_cell_concentration_impl.

real64 -> efp_t. Wraps efp_decompose; unlike the pre-fix version, the NaN / overflow flags are NOT discarded – they seed a%poison (nonzero iff is_nan .or. is_ovf), so a poisoned summand still taints every later efp_plus/efp_to_real even though the flags themselves aren’t returned here (callers needing the raw flags call efp_decompose directly).

Inverse of efp_to_transport: unpack a flat real64(EFP_TRANSPORT_WIDTH*n) buffer (post-collective, still exact integers as doubles) back into efp_t values, converting each bin AND the summed poison counter back to int64 and carrying the bins. ok = .false. iff any unpacked double is not an exact integer (would indicate the transport-exactness bound was violated) – the caller (halo_allreduce_efp_list) turns that into a fail-loud error stop, never a silent truncation.

Exact bin-wise integer subtraction, regularised. See efp_plus for why regularisation (not mere carry) is required for bit-for-bit order invariance, and for the poison propagation (additive here too – subtraction of a poisoned operand is still poisoned, never “cancels” back to clean).

Exact bin-wise integer addition, REGULARISED (not merely carried). Order-invariant BIT-FOR-BIT: efp_plus(a,b)%v == efp_plus(b,a)%v, and a running fold acc = efp_plus(acc, x_i) over any permutation of the x_i converges to the SAME raw bins (not merely the same reconstructed real) – test_efp_order_invariant asserts this on %v(:) directly, per the plan’s §9.2.

real64 = efp_to_real(efp_minus(a, b)) – the difference of two ~1e21-scale EFP totals resolved to the EFP quantum (2**-3P), NOT to ulp(1e21) as a double subtraction would give. This is the fix documented in the plan’s SS2.2: an implementer who converts both operands to real64 FIRST and subtracts loses the whole benefit of this module.

efp_carry plus: force every bin to share the overall sign, so a single well-conditioned FP accumulation (efp_to_real) can form Sum pr(n)*e(n) without alternating-sign cancellation error. Mirrors MOM6 regularize_ints.

Fail-loud (never silent) guard: a k-slab’s physical cell count must not exceed EFP_MAX_SUMMANDS, else a bin could overflow int64 before the next efp_carry. 1.34e8 is an 11500^2 single-rank layer – unreachable today, but the check costs one comparison at status cadence (CLAUDE.md: an unchecked bound “is a silent-corruption path exactly like the NZ_STACK_MAX one”).

efp_t -> real64. Regularises a LOCAL COPY of a (never mutates the argument) – pure with intent(in), per FORTRAN_STYLE.md’s “default new procedures to pure” (MOM6’s EFP_to_real instead mutates its intent(inout) argument).

Pack a list of efp_t values into a flat real64(EFP_TRANSPORT_WIDTH*n) buffer for a single collective (halo_allreduce_efp_list). Each bin, PLUS the poison counter, is transported as an EXACTLY-representable double (see EFP_MAX_RANKS): buf((i-1)*EFP_TRANSPORT_WIDTH + n) = real(list(i)%v(n)) for n = 1..EFP_DIGITS, and buf((i-1)*EFP_TRANSPORT_WIDTH + EFP_DIGITS+1) = real(list(i)%poison). Summing poison through the SAME MPI_SUM collective as the bins is what makes a poisoned rank’s contribution reach every other rank identically – see efp_t’s docstring.

One <Label> : <total> Error <residual> conservation line. With budget_on the residual closes the budget — (total − ref) + out − src — and the cumulative outflux / surface source are appended; without it the residual is raw drift total − ref and the line is byte-identical to the pre-budget format. label is a 4-char tag (Mass / Salt / Heat) so the colons align.

Flat-impl of the I1′ sweep for ONE tracer. Explicit-shape dummies so NVHPC does not walk a descriptor per launch.

Register the default diag-manager variable set and open the per-rank NetCDF stream. No-op when diagnostics are disabled. open_stream_file = .false. (the --validate-only path) runs the level/selection checks but creates neither the output directory nor the stream. Moved verbatim from driver_run_ocean’s private configure_ocean_diag helper (P2.4) so the API/bench setup paths can reach it too; F5 residual — now takes optional ierr instead of error stopping (see rdb_ocean_diag’s set_output_*_levels for the leaf-level conversions this threads through).

Map the ocean state onto the device, then (when bt_halo > 0) allocate + attach the wide-halo BT march-in workspace, then warm up the device-path halo exchanges (UCX cuda_ipc handle open, outside any timed region — see driver_run_ocean’s “Content-safety argument” comment, carried over unchanged).

Unwind device residency. Idempotent-adjacent: caller checks device_mapped (mirrors rdb_handle’s ocean_handle_t).

Host-side setup: decomposition -> grid -> god state -> IC seed -> restart (optional) -> the 21 configure_ocean_*-family stages -> ghost wraps -> halo init -> land mask -> wave drag/porous -> sea-ice IC. Does NOT map the state onto the device — call engine_enter_data next. Preserves driver_run_ocean’s exact stage order; see the per-stage comments below (carried over from the driver almost verbatim).

Advance one outer step at fixed dt, starting from simulation time t (consumed by tide/astro forcing inside ocean_dyn_step_split; the caller advances its own t_current by dt afterwards — see the module docstring for what stays caller-side). Wraps ocean_dyn_step_split (rdb_ocean_dyn.F90:1986) plus the per-step halves of the setup stages that must run BEFORE the dyn-core advance: file-forcing update/apply, boundary-data refresh, porous-area refresh.

Second half of one outer step: derive Q_heat/Q_salt from the surface-flux component set (no-op unless &ocean_forcing_nml enable_components) and, when diagnostics are configured, run one ocean_diag_t%step. Call immediately after engine_step and, when sea ice is enabled, engine_step_ice too — see engine_step’s docstring for why this is a separate call.

P2.4b: sea-ice per-step physics — ocean-side frazil accumulation, EVP dynamics (every outer step, independent of the thermo cadence: the ice’s own fast/slow split) plus the resulting ice->ocean stress coupling, and, at thermo cadence (engine%state%dyn%is_thermo_step()), category transport followed by the atmospheric-forcing / basal-flux / frazil-uptake / column-thermo / snowfall / brine / heat / shortwave / ITD chain. Transcribed VERBATIM from driver_run_ocean’s former inline sea-ice block — the internal call order is load-bearing (see the “MANDATED ORDER” comment below, itself carried over unchanged) and is NOT reordered here. Returns immediately, a no-op, when &ocean_ice_nml enable = .false. — bit-identical to before this phase.

Close the diag NetCDF stream (if one was opened), release the process-global ocean-halo and fold-exchange module state, then release host-side allocations. Call after engine_exit_data.

Populate the catalog at runtime (procedure pointers can’t be a parameter constructor). Idempotent.

Create path via libc mkdir(2) if it does not already exist (single path component — not a recursive mkdir -p; a missing grandparent still fails, loudly, at the inquire re-check below). Exists to close the headline P7 crash: &ocean_diag_nml/ &output_nml default output_dir = "./output" (rdb_config.F90:2316), so a first-time create() against a fresh checkout with no ./output/ directory would otherwise reach nc_create_file -> nf90_create -> ENOENT -> error stop, killing the whole host process.

Write the named liquidus coefficient SET onto the EOS handle. Called once, at configure time, from the earliest configure_ocean_* stage — before the flat-POD handle is copied onto the vmix / EPBL / kappa-shear / tidal-mixing slots and before ocean_state_enter_data, so every copy and every device kernel taking eos_t by value sees the configured set (same contract as rho0 / p_ref; no !$acc update device is owed).

Thermal-expansion and haline-contraction coefficients of the ACTIVE equation of state at a point, in the SAME DIMENSIONAL convention the eos_t members alpha_T / beta_S carry:

Variant dispatch on bare 3D arrays — the state-agnostic body of ocean_eos_compute (the shim hosted in rdb_ocean_eos_compute), which forwards the multilayer state’s registry arrays here. Kept free of state-type dependencies so the dispatch stays callable from a bare array context. All impls evaluate at the SINGLE, HORIZONTALLY UNIFORM reference pressure eos%p_ref (&ocean_eos_nml p_ref, default 0 ⇒ surface/potential density) with the H_VANISHED vanishing-layer fallback to eos%rho0.

Density sensitivities ∂ρ/∂T and ∂ρ/∂S of the ACTIVE equation of state at a point — the signed twin of eos_buoyancy_coeffs, which is where the per-variant closed forms live:

Scalar density evaluation at a point — the same formulas the 3D eos_*_impl kernels apply, exposed for finite-difference verification of eos_specvol_derivs and for host-side diagnostics. Takes the shared eos_t handle by value. The else is unreachable-by-contract (see eos_validate).

eos_density_point AND eos_specvol_derivs at the same point from ONE evaluation of the active EOS. The isopycnal-slope and Redi builders need both (drho/dX = -rho^2 * dSV/dX, locally referenced); calling the two routines separately evaluated the EOS twice – under Roquet, two full roquet_spv_point calls, the first of which threw its derivatives away. Each branch uses the same expressions as the two routines it fuses, so rho, dsv_dt and dsv_ds are the numbers they return.

Seawater freezing point T_f (degC) at salinity S and pressure p — the ocean-side prerequisite for the sea-ice port (PLAN_SEA_ICE.md, PR 1). Same point-function style as eos_density_point (flat-POD eos_t by value, device- callable), plus elemental so callers can evaluate whole salinity arrays in one reference.

Linear two-tracer EOS, flat-impl form:

Roquet et al. (2015) SpV EOS evaluated at a single reference pressure p_ref (a SCALAR by design — see the horizontal-uniformity contract in eos_compute_arrays). Same outer-shim signature as eos_linear_impl/eos_wright_impl — bare 3D arrays, model (PT, SP) tracers, vanishing-layer fallback to rho_0.

Analytic specific-volume sensitivities dSV/dT and dSV/dS (SV = 1/rho) at a point. Needed by the EPBL energy bookkeeping (pressure-weighted PE-per-unit-tracer-change weights) and kappa-shear buoyancy — dSV/dX = -(1/rho^2) d(rho)/dX.

Host-side fail-loud gate over the device-supported variant set. Device point routines (!$acc routine seq) cannot error stop, so membership in the device-callable set is guaranteed HERE at configure time; the device else branch is then unreachable-by-contract.

Wright (1997) rational EOS evaluated at a single reference pressure p_ref, which is a SCALAR by design — see the horizontal-uniformity contract in eos_compute_arrays. Same outer-shim signature as eos_linear_impl — bare 3D arrays, vanishing-layer fallback to rho_0.

FV-Wright PGF column sweep. Top-down per-column traversal that simultaneously produces the hydrostatic pressure stack p_edge_out and the in-situ density rho_insitu_out at each layer centre.

Per-column EPBL solve. One do concurrent (j, i) with the serial work in k inside (j -> i -> k ordering); ALL sweep state is carried in scalars (design doc D6) — the only column arrays are the six iteration-invariant workspaces filled by the prep sweep.

Run EPBL over the domain: fill this%kd_int (interface diffusivity) and this%mld. Call at thermo cadence with the thermo dt. Outer shim: dereferences the tracer-registry hTr arrays + the 2D Q_heat/Q_salt forcing fields on the host (array-of-DT and allocatable indirection blocks NVHPC device codegen), then forwards to the column kernel which reads q_T_kin(i,j) / q_S_kin(i,j) per column inside the DC.

mstar = (mechanical TKE available for entrainment) / u^3. Schemes: constant; OM4 Ekman/Obukhov balance; RH18 fits. All followed by the optional convective reduction (mstar_conv_adj in [0,1]; the u=0 corner multiplies by (1 - adj), matching the reference behaviour).

Turbulent Langmuir number La = sqrt(u*/u_s_SL): the Stokes drift averaged over the surface layer of thickness zsl under the Phillips spectrum, in the singularity-safe form (Breivik et al. 2016 with Webb & Fox-Kemper 2015 directional spreading). No MIN_LANGMUIR / LA_DEPTH_MIN floors — those belong to MOM6’s profile-averaging wave paths, not LF17.

LF17 statistical wave state from the water-side friction velocity alone: COARE 3.5 fixed-point inversion u* -> U10 (Edson et al. 2013), then the Pierson-Moskowitz-based surface Stokes drift and Phillips peak wavenumber (Li & Fox-Kemper 2017; Breivik et al. 2016). BLD-independent — call once per column, outside the MLD iteration.

Apply the Langmuir enhancement to mstar (Reichl & Li 2019; Li et al. 2016 stability modification). The modified Langmuir number folds the boundary-layer stability regime in via Ekman / Obukhov / MLD length-scale ratios, split by the sign of the surface buoyancy flux; all-zero lac coefficients give La_mod = La exactly.

Fold the EPBL diffusivity into the vmix interface fields. Called EVERY stage (the interior closure rewrites kv/kt each stage; kd_int itself refreshes at thermo cadence). Interior interfaces only — k=1 (bed) and k=nz+1 (surface) stay at the closed-boundary zero in both source and target.

Mixing-length shape factor in [translay_scale, 1]: 1 at the surface, decaying to the transition-layer floor at the MLD. shaped = .false. (no MLD iteration) returns 1.

Equilibrium argument V_c (radians) for every catalog constituent at day number dnum. s,h,p,N are taken in radians (deg-mod-360 -> rad); the result is left un-modded (cos/sin are periodic). The ±pi/2 diurnal signs are load-bearing.

Reset the ring to empty.

Number of live messages in the ring (0..ERROR_RING_SLOTS).

Message at ring-relative index i (0 = most recent, 1 = next-most-recent, …). Empty string if i is out of [0, error_ring_count()).

Push one message; error_ring_get(0) is always the most recent. Overwrites the oldest slot once full.

fxoc *= mask_u/evp_sub_steps, fyoc *= mask_v/evp_sub_steps (:1415-1441).

mask_t: wet_T inside the physical domain AND in every ghost band that is not pinned; a ghost band on a physical, non-periodic edge (land_*) is 0 (SIS2 mask2dT semantics — pins a wall edge’s ghosts to land, matching SIS2’s own domain-edge convention even when a driver run’s wet_T ghost happens to read 1). An MPI seam ghost follows wet_T, which the ocean setup exchanged; the local periodic wrap (wrap_x/_y, only on an axis the halo does not own) then overwrites a single-rank periodic band, so on one rank the result is the same wrap-or-0 mask as before. mask_u(i,j) = mask_t(i-1,j)*mask_t(i,j), mask_v ditto in y, mask_q = product of the 4 surrounding mask_t (SIS2 mask2dBu).

u_tmp = ui (full array — the v-momentum MUST read pre-update u, D4). Explicit do concurrent element copy rather than a bare whole-array assignment (repo convention for device-resident arrays — see rdb_ml_dynamics’s h_layer0 save pattern).

Interior copy of the gathered mis/mice/ci, masked to mask_t (defence-in-depth beyond the caller’s own wet_T gate); ghost rows/cols zeroed (the periodic wrap that follows fills them).

mi_u(i,j) = 0.5*(mis(i-1,j)+mis(i,j)); mi_v(i,j) = 0.5*(mis(i,j-1)+mis(i,j)) (SIS2 :967-974, rdb index translation §1). Array-edge faces (i=1/i=nx+1, j=1/ j=ny+1) have no neighbour on one side; mis at those ghost rows/cols was already periodic-wrapped or zeroed, so a naive mis(i-1,j)/mis(i,j) read is always in-bounds here EXCEPT at the two hard array edges themselves — those faces are handled explicitly.

pres_mice = p0_rho*exp(-c0*max(1-ci,0)) (:878); dxharm = 2*dxT*dyT/(dxT+dyT); del_sh_min_pr = 2*del_sh_min_scale*dt^2 / (Tdamp*dxharm^2) guarded on dxharm > 0 (:880-890).

PR 36: PROJECT_ICE_CONCENTRATION (SIS2 SIS_dyn_cgrid.F90:1064- 1077). ci_proj = ci*exp(-dt_cum*sh_dd) then pres_mice = p0_rho*exp(-c0*max(1-ci_proj, 0)). ci_proj is a local() scalar, NOT an array: SIS2 materialises it only for the sigI/sigII/find_ice_strength diagnostics Roundabout does not have (documented divergence). del_sh_min_pr is NOT recomputed here (it has no ci dependence, evp_pres_mice_impl above). ci_proj is deliberately unclamped above 1 – max(1-ci_proj, 0) already saturates the effect at p0_rho, and for dt_cum*|sh_dd| > 709 (an unreachable regime in any sane run) exp overflows to +Inf, max(1-Inf, 0) = 0, exp(0) = 1 – IEEE launders the overflow to exactly the correct saturated value, so no guard is needed (SIS2 has none either).

q(ic,jc) = f_corner*tot_area / (Σ areaT*mis over the 4 cells + tot_area*m_neglect) (:977-982); mi_ratio_A_q via ice_evp_mi_ratio_point (requirement 5). 4 T-cells around corner (ic,jc): (ic-1,jc-1) (ic,jc-1) (ic-1,jc) (ic,jc) (SW/SE/NW/NE, §1). Array-edge corners (no T-cell on one side) get q=0/mi_ratio=0 (land-corner convention — consistent with mask_t=0 beyond the array edge).

sh_Dt / sh_Dd at cells (:1053-1061).

sh_Ds at corners (:1045-1050) — requirement (4): the SINGLE scalar no-slip factor (2-mask_q) on the WHOLE combined strain. Computed over the interior+1 ring (ic,jc in [1,nx+1]x[1,ny+1] — the full corner array; out-of-band neighbours contribute 0 via zero ghost velocities at the hard array edges, never per-term mirroring).

str_s relax (:1137-1143). Corners in [1,nx+1]x[1,ny+1]; the 4 surrounding T-cells at an array-edge corner are handled by zeta’s own ghost values (zero-mass ghost cells => zeta=0 there, contributing nothing) — no special-case branch needed.

str_d/str_t semi-implicit relax (:1124-1134), non-weak_low_shear branch only.

PR 36: the FINAL CFL clip (SIS2 :1443-1500) – TRUNC_BACKOFF (0.95) back-off instead of the exact bound, PLUS a count of the ice-bearing faces it touched (mi > m_neglect, SIS2 :1466,1469 – massless faces clip silently, matching SIS2: counting them would flood the driver’s warning with meaningless ice-free clips at every margin). Not a do concurrent: reductions use !$acc parallel loop reduction(...) (ice_compress_impl is the local precedent for a reduction that also mutates the arrays it walks). Same bound algebra as evp_truncate_velocity_impl, duplicated rather than shared: the in-loop variant runs evp_sub_steps (432 by default) times per outer step and must NOT carry a reduction (each would be a device->host sync); this variant runs once and must. CLAUDE.md’s “duplicate explicitly” rule – merging the two costs 432 syncs per outer step.

PR 36: the shared CFL-clip algebra – the transport-CFL bound on the ice velocity (SIS2 SIS_dyn_cgrid.F90:839-870, the in-loop half at :1338-1361, the final half at :1443-1500; this routine is the counting-free, caller-chosen-backoff form both reuse; evp_truncate_final_impl below wraps it with the 0.95 back-off and the mi > m_neglect count).

u-momentum (:1172-1231, requirement 1: fxic_now carries the FULL str_t force term). Loop over u-faces ng+1..ng+nxp+1 x ng+1..ng+nyp — each iteration writes only its own face.

v-momentum (:1257-1334, mirror of u). D4: reads u_tmp (the PRE-update u), never the just-updated ui. Minus on the str_t divergence term (:1263-1267).

Counted allocatable footprint of the EVP scratch (0 when unallocated, i.e. whenever &ocean_ice_nml dynamics is off).

Reverse of init (host deallocation). Idempotent — safe on an already-clean workspace. Device unmapping is exit_data’s job (call it first).

Attach the scratch to the device with create — every array is pure scratch, fully (re)written before it is read on every ice_evp_dynamics call (masks rebuilt, cell fields filled, the subcycle kernels overwrite the rest), so create (not copyin) is correct. One-level components of the type(...) dummy — same shape as ocean_sea_ice_enter_data_impl, no associate-leaf needed.

Reverse of enter_data_impl — pure scratch, delete throughout.

Allocate the EVP scratch for an (nx,ny) T-cell domain, zero-init. Called ONCE per workspace lifetime — from ocean_sea_ice_init (gated on dynamics) and from the EVP test harness. nx/ny declared before the arrays that use them is unnecessary here (all allocatables), but the extents are cached for destroy.

Periodic ghost-wrap for a corner-staggered field (e.g. str_s), shape (nx_total+1, ny_total+1). No corner-wrap helper exists in rdb_ocean_periodic (only centre/face_x/face_y) — this is the EVP-local twin, same two-pass (x-then-y) structure.

SIS2 :899-907 — zero ice velocities where BOTH neighbouring cells are massless (or the face is masked/land).

Zero the subcycle-averaged ice->ocean stress accumulators via an explicit do concurrent device kernel (F1): fxoc/fyoc are copyin-mapped device-resident arrays, so a host = 0.0_wp would zero only the HOST copy and leave the device copy carrying the prior call’s average (the accumulate below would then converge to S/(N-1) instead of S/N). Runs on the device-present arrays; inert no-op on host builds.

del_sh / zeta (:1082-1095). shear_at_T averages the 4 surrounding corner sh_Ds values.

East neighbour rank with periodic wrap-around.

West neighbour rank with periodic wrap-around.

MOM6 fac_fn: the WENO-Z nonlinear factor (1+tau/b)^2, EXACT-clamped (not an additive epsilon) to a dominant value when the smoothness indicator b is degenerate (|b| <= eps*tau), so a locally constant stencil takes over without a 0/0. Shared by weno3/5/7.

face_depth_mean_u with MOM6 wt_u weighting (&ocean_bt_nml forcing_visc_rem): the weight is h_face·visc_rem(k) instead of h_face, so layers the implicit vertical-friction solve will immediately damp (grounded sliver stacks under the BBL glue, visc_rem → 0) contribute nothing to the barotropic forcing. Without this the spurious grounded-layer PGF’s depth-mean drives the fast loop ballistically even after the layer velocities themselves are glued (PGF_BUG.md §9). Denominator falls back to zero-output on an all-remnant-zero column (the substep should not force an immobilized column).

Symmetric v-face counterpart of face_depth_mean_rem_u — same MOM6 wt_u floor, ieee_is_finite-guarded the same way.

Depth-average a u-face 3D field, weighted by the per-layer face thickness frhat_h_face_step returns (&ocean_bt_nml frhat_scheme; FRHAT_ARITHMETIC = the plain mean of the two abutting cell columns’ h_layer values, today’s default; the argument is otherwise unchanged). Writes to a 2D field at the same u-face shape. Wall faces (i=1, nx+1) fall back to the single available cell (frhat_h_face_step degenerates there for either scheme — see its docstring).

Symmetric v-face counterpart of face_depth_mean_u.

Public only for the unit-test suite (no production module imports it); ignore when developing production code in other modules. Face-thickness for the vdiff implicit operator.

Push msg to the error ring, log it (unchanged behaviour), then either set ierr = code and return (caller must return immediately after this call) or error stop msg when ierr is absent. See the module docstring for the full contract.

Layer ideal age = tracers(idx_age)%hTr / h_layer (s). Read-out only; registered when &ocean_tracers_nml enable_ideal_age (idx_age > 0).

Invoke the var’s fill (and remap if non-LAYER) so output_buffer holds the current sample. The accumulation / log emission step consumes output_buffer after this.

Fill ghost cell bathymetry by constant extrapolation from the nearest interior cell. This ensures boundary flux computations see a consistent bottom elevation across the ghost-interior interface. Works on a target 2D array directly so both the coastal and ocean barotropic slots can reuse it.

The historical constant-piston glue (bbl_per_face = .false.).

Device-side scalar fill. Used by reset_accumulator and (when idx_temperature / idx_salinity is unset) the manager could also seed output_buffer via this; today only the accumulator reset goes through here.

Per-column CAVITY_MELT_* status code, as a real. 0 is OK; anything else is a column that took the kernel’s documented zero-melt safe state, and the console’s warning line says how many there were. Masked like the rest, so the plane cannot be read as “everything is fine” over open ocean.

Fill a cell-centre |Coriolis| array via metrics_fill_coriolis (D7). beta_plane (default) is BIT-IDENTICAL to the legacy EPBL / kappa-shear set_f_centre. The corner output is discarded here (filled into local scratch).

Fill a C-grid corner Coriolis array via the single generator-driven routine metrics_fill_coriolis, honouring &ocean_grid_nml coriolis_scheme (D7). beta_plane (default) is BIT-IDENTICAL to the legacy coriolis_adv_set_beta_plane; planetary uses the metrics geography. The centre output is discarded here (filled into local scratch).

Haline exchange velocity gamma_S (m/s) — see fill_exch_vel_t.

Thermal exchange velocity gamma_T (m/s) the solve CONVERGED on — not Gamma_T·u* re-derived here. Under hj99/yung25 the two differ by the stratification suppression.

Periodic-x wrap + north-fold (scalar copy) of the static corner Coriolis array for a tripolar grid. f is reflection-invariant (same latitude at the conjugate corner), so negate=.false.

2D counterpart of fill_getter_3d.

Shared tail of every 3D P2 getter: base-address c_loc, actual array extents (never assumed from grid metadata — read straight off the pointer), and the generation stamp. arr must already be pointer-associated with a live, non-empty state array.

Reservoir-based ghost fill, east edge.

Sign-switch ghost fill for the east edge.

Reservoir-based ghost fill, north edge.

Sign-switch ghost fill for the north edge.

Reservoir-based ghost fill, south edge.

Sign-switch ghost fill for the south edge.

Reservoir-based ghost fill, west edge. hTr_ghost = tres_w(j, k, it) * h_ghost — unconditional.

Sign-switch ghost fill for the west edge.

Zero-gradient fill of h_layer ghosts at open-ish edges.

S* = S_far − S_b (g/kg) — the salinity contrast the haline exchange acts on. Positive under melting (meltwater freshens the interface).

Total sea-ice concentration (0..1) at T-centres — the two-mode per-cell gather of ice_cell_concentration_impl (rdb_ice_state), inlined (see fill_ice_conc_thick_impl). Registered by register_default_diags only when &ocean_ice_nml enable. Public only for the unit-test suite.

Shared conc/thick device kernel. Inlines the two-mode gather of ice_cell_concentration_impl (rdb_ice_state — convention of record; test_ocean_ice_diags pins the copies equal): ncat==1 legacy lumped (per-CELL m_ice, ci = 0/1), ncat>1 SIS2 ITD (ci = min(1, Σ part_size), mice = Σ part_size·m_ice). emit_thick=.false. ⇒ buf = ci; .true. ⇒ buf = mice/ICE_RHO_ICE (grid-mean thickness, m). Scalar flag branch is constant-folded on the device — one kernel, no scratch companion.

Sea-ice drift speed |u_ice| at T-centres (m/s) — C-face pairs averaged to centre, 2-D twin of the ocean fill_ke stencil. Ice off / not init ⇒ zeros (never registered by default; opt-in via &ocean_diag_nml diags).

Grid-mean sea-ice thickness (m) at T-centres: mice/ICE_RHO_ICE (MOM6 effective-thickness convention) — the two-mode per-cell gather of ice_cell_concentration_impl, inlined (see fill_ice_conc_thick_impl). Registered by register_default_diags only when &ocean_ice_nml enable. Public only for the unit-test suite.

Eastward sea-ice velocity at T-centres (m/s) — u-face pair averaged to centre, 2-D twin of fill_u_centre. Ice off / not init ⇒ zeros.

Northward sea-ice velocity at T-centres (m/s) — v-face pair averaged to centre, 2-D twin of fill_v_centre. Ice off / not init ⇒ zeros.

EPBL interface diffusivity (m^2/s). Buffer is layer-shaped (nx, ny, nz); we emit the value at the BOTTOM interface of each layer (kd_int(:, :, k) convention), losing only the identically-zero surface interface.

Kappa-shear interface diffusivity (m^2/s). Buffer is layer-shaped (nx, ny, nz); we emit the value at the BOTTOM interface of each layer (kd_int(:, :, k) convention), losing only the identically-zero surface interface.

KE per cell on the C-grid: KE = 0.25·(u_W² + u_E² + v_S² + v_N²), the discrete C-grid KE-density (consistent with the Coriolis KE_ARAKAWA stencil). Public only for the unit-test suite.

Depth-integrated kinetic energy at each cell: KE_total(i, j) = Σ_k 0.5 · h_layer(k) · (u_c² + v_c²)

Broadcast the scalar viscosity kappa into the interface-located workspace. Boundary interfaces (bed at k=1, surface at k=nz+1) are forced to zero to match the closed BCs the column solve assumes; the original scalar kernel hard-coded those BCs via α_1 = 0 and β_nz = 0.

Basal melt mass flux (kg m-2 s-1), positive = melting.

Basal melt rate in the ISOMIP+ reporting unit (m yr-1 of ice-equivalent freshwater) — see melt_m_per_yr_factor.

Mixed-layer depth via the de Boyer Montégut threshold.

Per-column scan from the first LIVE layer (k_top, which is nz_ml unless a rigid top has vanished the layers above it) toward the bed; first layer whose ρ exceeds (surface ρ + MLD_DENSITY_THRESHOLD) marks the MLD as the cumulative h-sum above it. No crossing → MLD = full column depth. Threshold met at the surface itself → MLD = 0.

EPBL active-mixing-layer depth (m) into the k=1 plane.

Layer pseudo-salt = tracers(idx_pseudo_salt)%hTr / h_layer (psu). Read-out only; registered when &ocean_tracers_nml enable_pseudo_salt (idx_pseudo_salt > 0). Mirrors fill_age.

Pseudo-salt deviation D = pseudo_salt - S (psu): a direct, measured proxy for how far the passive-tracer transport path has drifted from the active-tracer (salinity) path. Gated on BOTH indices being registered.

In-situ density per layer — direct read of the EOS slot (driver must have run ocean_eos_compute this step; not re-invoked here), with the NaN missing-data sentinel on every vanished layer.

buf = rho_layer on live layers, IEEE NaN on vanished ones.

Layer salinity, same pattern as temperature.

Interface salinity S_b (g/kg).

SSH = total column thickness minus bathymetry depth, into the k=1 plane of buf. Multilayer or barotropic path per use_multilayer. Public only for the unit-test suite.

Interface temperature T_b (degC) — on the liquidus at S_b.

Layer temperature = tracers(idx_temperature)%hTr / h_layer. Tracer-registry indirection dereferenced HOST-side before the flat-impl kernel (array-of-DT deep deref blocks NVHPC device codegen). Public only for the unit-test suite.

T_f(S_far, p_top) (degC) — the in-situ freezing point of the FAR FIELD at the ice base. Distinct from tbdry, which is the freezing point of the INTERFACE salinity S_b; their difference is the whole three-equation correction, so shipping both makes it visible.

T* = T_far − T_f(S_far, p_top) (degC) — the far-field temperature above the IN-SITU freezing point at the interface pressure. Positive drives melting. This is the single number the melt rate is roughly linear in, so it is the first thing to look at when a melt rate looks wrong.

Upwind-aware tracer ghost fill for south and north open edges.

Zero-gradient CONCENTRATION fill of a tracer’s hTr ghosts at open-ish edges over the FULL cross-extent (ghost×ghost corners included): hTr_ghost = (hTr_int / h_int) * h_ghost. Runs BEFORE the upwind-aware per-edge fill, so corners keep this zero-gradient value (no corner T/S blow-up under inflow). Same gating/corner coverage as fill_h_layer_ghosts.

Upwind-aware tracer ghost fill for west and east open edges. West inflow : u(i_w,j,k) > 0 → ghost = clamped_tr * h_ghost West outflow : u(i_w,j,k) <= 0 → ghost = interior (zero-gradient) East inflow : u(i_e,j,k) < 0 East outflow : u(i_e,j,k) >= 0

Shared flat-impl for any tracer concentration field.

Depth-integrated zonal transport (m²/s): T_x(i, j) = Σ_k h_c · u_c (cell-centre form)

Depth-integrated meridional transport (m²/s) — mirror of x.

C-grid u-face → cell centre by simple 2-point average.

Melt friction velocity u* (m/s), max(sqrt(C_d(u²+v²+u_tide²)), u*_min). This is the MELT law’s u* and nothing else: it drives no momentum drag, and KPP/EPBL do not read it.

Zero-gradient fill of per-layer face velocities into the OPEN-edge ghost region, corners included. u_layer (nx_total+1, ny_total, nz); v_layer (nx_total, ny_total+1, nz). x-pass over the full cross extent; y-pass clipped to the physical span unless the adjacent x-edge is open. Gated on open-ish tags ⇒ WALL/PERIODIC ⇒ no DC ⇒ bit-identical.

The model’s OWN relative vorticity ζ = ∂v/∂x − ∂u/∂y, at the same C-grid corners coriolis_adv computes it at (circulation/area form, C1 slip factor), averaged corner -> cell centre. See fill_vorticity_z_impl for the formula and why it replaced a centred T-point stencil that differenced ocean velocity straight against the zero stored at land faces (a spurious no-slip-like vorticity sheet at every coast, even under the free-slip default).

Cell-centred relative vorticity that MIRRORS the dynamics’ own corner ζ bit-for-bit — the shared rdb_rvc_zeta_corner body (src/shared_module_utilities/rdb_rel_vort_corner.inc), the same circulation-form formula coriolis_adv_compute_tendencies Pass 1 evaluates inline for its q_corner scratch (src/core/ocean/kernels/coriolis_adv/rdb_coriolis_adv.F90), incl. the C1 slip factor: no_slip=.false. (default, free-slip) masks a land corner (wet_q=0) to zero rel-vort; .true. (no-slip) gives it the 2-wet_q image-vorticity value instead.

Water-column thickness under the shelf, D = bt_H_ref + bt_eta (m) — the ONE expression every consumer of the column depth uses, which is exactly what the cavity datum (bt_H_ref = b − z_draft) buys. Masked by cover_frac: the quantity is perfectly well defined in open water, but as a CAVITY diagnostic a domain mean should be the mean cavity thickness, not that diluted by the open basin.

Ice-shelf draft (m, POSITIVE DOWN from the geoid) — geometry, so it needs only &ocean_cavity_dyn_nml and is masked by cover_frac rather than by the melt slot’s active. Beyond the calving front there is no draft, which is missing data, not a draft of zero.

Device-side zero of buf — used when a fill’s input is missing (e.g. tracer not registered) so the fold reads a defined value.

Copy the accumulator into output_buffer. Called at cadence-fire BEFORE log / NetCDF emission consumes output_buffer.

Device-side out = accum * scale_factor. Used for MEAN finalise where scale_factor = 1/dt_accum.

Marginal zonal face area d(uhbt)/du. At u = 0 returns the average of the two cubic-branch slopes (discontinuity harmless — only consumed via max(…, h_neglect)).

Meridional mirror of find_duhbt_du.

Zonal mass transport through a u-face given face velocity u. C¹ continuous in u (cubic near zero, linear saturation beyond).

Resolve a variable id: try the namelist override name first (when non-blank), then the documented fallbacks in order. Error-stops with a descriptive message when none is found.

Meridional mirror of find_uhbt.

Format an integer.

Format a logical as a JSON true/false literal (as opposed to fmt_logical’s Fortran .true./.false.).

Format a logical as namelist .true./.false..

Format a real in a stable, namelist-valid exponential form.

Cell-centred 2D north fold (copy).

Cell-centred 3D north fold (copy).

SW-corner 2D north fold (+ fold-line projection); negate for a true-vector component, copy for a scalar / vorticity.

Fill the north halo of a 2D cell-centred field by the T-fold.

3D T-fold halo-fill — identical per level.

2D Bu-corner fold: north-halo fill + on-line projection.

Fill the north halo of a 2D x-face (Cu) field. Sign-flipped by default (negate absent / .true., the true-vector contract — wind stress, velocity). negate=.false. copies instead: for a SCALAR carried on a u-face (e.g. the viscous remnant visc_rem_u — a fraction, not a flux component), the 180-degree fold rotation still swaps which side of the seam the value sits on, but the value itself does not change sign (see fold_north_corner_2d’s negate for the matching corner-stagger contract).

3D x-face (Cu) north-halo fill, per-level identical. See the 2D twin for the negate (vector vs. scalar) contract.

2D y-face (Cv) fold: north-halo fill + on-line antisymmetric projection at the fold row. Used for the barotropic bt_vbt field in the BT fast loop. negate (default .true., see the u-face twin) selects true-vector (sign flip + self-conjugate zero) vs. scalar (copy + self-conjugate left unchanged, matching fold_north_corner_2d’s scalar contract) — a scalar on a v-face (e.g. visc_rem_v) is the SAME physical attribute of the SAME face seen from both sides of the seam, so the duplicated DOF at the fold line must agree, not cancel.

3D y-face (Cv) fold: north-halo fill + on-line antisymmetric projection at the fold row. See the 2D twin for the negate (vector vs. scalar) contract.

Build tile rx’s send/receive lists for every peer and both column families. Pure: every north-row rank computes every tile’s receive lists itself (decomp_init arithmetic), so no handshake.

Free the plan’s lists and reset it to the empty default. Safe on a plan that was never built.

Every destination column i = 1..ncol of tile rx’s window for a column family: the owning tile own(i) of its mirror, the owner’s local storage column src(i) to read, and the fold-row class cls(i). ncol = w+2ng (T) or w+2ng+1 (U); the arrays must hold at least that many entries.

Storage rows of message row r (1..fold_stagger_nrows) on a north-row tile of nyl physical rows: the sender reads src_row, the receiver writes dst_row. * T, u: r = d = 1..ng: src ng+nyl+1-d, dst ng+nyl+d. * v, corner: r = d+1, d = 0..ng: src ng+nyl+1-d, dst ng+nyl+1+d (d = 0 is the fold-line row, src = dst = ng+nyl+1).

Combine the current output_buffer sample into accumulator per time_op.

Masked fold — multiplies sample by weight(i, j) before folding. MAX/MIN treat weight == 0 as “don’t update” so masked-out cells keep their seed value. Loop bounds clip to whichever extent is smaller (accumulator vs mask) so an undersized mask doesn’t OOB.

Whole-buffer fold without a region mask. One do concurrent per op so the compiler can specialise — case-inside-loop blocks NVHPC device codegen.

Column family of a stagger (FOLD_FAM_T for T/v, FOLD_FAM_U for u/corner).

Rows a stagger moves per column: ng (T, u) or ng+1 (v, corner: the fold-line row is always sent, the receiver uses it only where recv_cls == FOLD_ROW_WEST).

First global cell a and width w of tile rx — the decomp_init split (remainder to the WEST tiles), restated here so the plan stays pure and dependency-free (cross-checked against decomp_init by the unit test).

Tile holding global cell c (1..ni), closed form of the decomp_init split.

x-face (u) 2D north fold (negate).

x-face (u) 3D north fold (negate).

y-face (v) 2D north fold (negate + fold-line projection).

y-face (v) 3D north fold (negate + fold-line projection).

&ocean_forcing_nml enable_components=.true. with &ocean_thermo_nml enable_thermodynamics=.false. is a configuration that validates but does nothing (the component set only feeds the assembler->apply_tracers thermo path) — named predicate per the lateral_closure_is_implemented idiom so validate_config reads as one line.

MOM6 near-bottom hWght mass-weighting fractions of one face (MASS_WEIGHT_IN_PRESSURE_GRADIENT): the distance by which the layer top of one column sits below the other column’s bed, scaled by the squared relative thickness contrast. Zero ⇒ plain linear interpolation (hwt_ll = hwt_rr = 1, hwt_lr = hwt_rl = 0).

Render canonical_diag_gate_hint(name) into a human-readable sentence for the “requested but not registered” warning. Phrased as a QUESTION, never a diagnosis: the hint is a hand-maintained mirror of register_default_diags’s gates (§11.2 of the PR-64 plan) and can drift from the true cause. An empty hint means the diagnostic is one of the four unconditional ones (SSH/u/v/KE) and SHOULD have registered — that is a bug, not a closed gate.

Safe accessor: return clamped_tracer(it) or 0 if unallocated/out of range.

MOM6 bottom-blocking (“Avoid moving dense water upslope from below the level of the bottom on the receiving side”). sfn is the unlimited streamfunction at an interface: the transport of everything BELOW it, > 0 from L to R. Its donor layer is the one just below the interface on the donor side ([e_bot, e_top]).

Fill the 2D face KhTh fields from the per-face base, CFL-clamped per face (native u-face idxCu/idyCu, v-face idxCv/idyCv) and zeroed on wall faces. Base = VarMix khth_ext_* when use_ext, else the scalar khth; the same diffusive-CFL min is applied either way.

Clamp a slope to +/- smax (bounded slope for the PE release). A non-finite slope reads as 0 (no release) rather than being laundered into ±smax by the NaN-blind clamp (CLAUDE.md).

u-face GM streamfunction + bolus-transport column recurrence. Interior u-face (i=2..nx) pairs columns iw=i-1 (west) and i (east). Bottom-up sweep: interior interfaces Kr=2 (bed-most) -> nz (surface-most), uhtot=0 at the bed; the surface BC (Sfn=0 at Kr=nz+1) is closed after the loop by uhD(nz)=-uhtot, giving Sum_k uhD=0 exactly. Interface Kr straddles ka=Kr (above) and kb=Kr-1 (below) and fills LAYER kb; the rsum bound keys on ka, the donor h_frac and per-layer cap on kb.

v-face GM column recurrence — mirror of gm_column_x with the v-stagger. Interior v-face (j=2..ny) pairs columns js=j-1 (south) and j (north). See gm_column_x for the open-column (use_open) rule and the MOM6 nk_linear divergence note.

Flat-impl GM kernel. Passes: CFL-clamp the 2D face KhTh, the u-face and v-face column recurrences into uhD/vhD, the gm_src PE release. Each face’s column sweep is serial in k (the uhtot recurrence) but faces parallelise over (i,j).

Fill uhD/vhD (m^3/s) with the GM bolus thickness transport and gm_src with the PE release, from the CURRENT ms%h_layer and the stored slopes. Run every outer step AFTER the dynamics, and immediately followed by continuity_gm_apply with the SAME dt and the same, untouched h_layer: the availability cap A·(h − H_VANISHED)/(4·dt) bounds that h and no other (see the module docstring). No-op when disabled, uninitialised, or the slopes slot is absent/disabled.

Donor mass fraction h_avail(k)/rsum_above(k) (0 when no mass is available above). rsum_above(k) is the cumulative availability from the surface down to and including layer k, so hf in [0,1].

GM potential-energy release at cell centres for the MEKE seam: gm_src = 1/4 * Sum_k rho0 * (KHSlope^2N^2) * h over the four straddling faces, summed over interior interfaces; slope clamped to slope_max, N^2 floored at 0. gm_src >= 0 for a stable tilted column. An interface value is attributed to the layer below it (kb=K-1); bed + surface carry zero slope/N^2.

max(N^2, 0) for the PE release, with a non-finite N^2 read as 0 (a bare max is NaN-blind under relaxed FP, CLAUDE.md).

.true. when the GM block below this step already calls varmix_compute (GM present + enabled + slopes present, at a thermo step) — so the standalone Gap-1 resolution-function refresh must NOT run it again (avoids a double compute and preserves the GM+VarMix slope ordering bit-for-bit).

Zonal PPM tracer advection of every horizontally-advected tracer by the GM bolus flux uflux (the resolved path’s own kernel, tracer_advect_zonal_one_impl; heat/salt budgets weighted by budget_w). No OBC ghost override: continuity_gm_apply has closed every non-periodic edge face.

Meridional twin of gm_tracer_advect_x.

.true. unless variant == OPGF_VARIANT_GPRIME with nz /= 2 (PR-6 fail-loud). The reduced-gravity gprime PGF hard-writes ONLY k=1 (bottom) and k=2 (top) — with nz > 2 layers k=3..nz carry zero pressure gradient AND the “top” branch lands on layer 2 of nz (the abyss under the bottom-up convention); with nz < 2 the kernel early-returns leaving the whole PGF zero. The guard is gprime-specific — every other variant supports general nz, so this returns .true. for them regardless of nz. Wired into validate_config against cfg%nz_layers (config-time).

Great-circle distance (m) between two geographic points (deg), via the haversine formula (numerically stable for short arcs).

Proleptic-Gregorian Julian Day Number (integer, at 00:00 UT). Fliegel & Van Flandern algorithm.

Add a key to the group (sourced-allocate into a box). Error stop on duplicate key name (programming error).

Index of name in the group’s key list, 0 if absent.

All key names of a group as a char array.

All registered + external group names as a char array.

(Re)allocate both buffers for nslab slabs per peer (exit-delete / dealloc / alloc / enter-create, as ocean_halo_buffers_ensure_nz). Never while a group is open.

OPEN-column upstream face thickness, for upstream_h_face under &vcoord_nml zfixed_closed_faces.

Order-invariant EXACT cross-rank combine of nval EFP values in ONE collective (PR-32). Replaces N separate scalar halo_allreduce_sum calls with one packed allreduce.

MPI_Allreduce with MPI_MAX — max-type reductions are exact in FP, so a global max stays layout-reproducible (ocean-MPI plan D5). Used for auto_n_inner’s global gravity-wave CFL (shared n_inner).

MPI_Allreduce with MPI_MIN for global timestep

MPI_Allreduce with MPI_SUM for CG dot products

Cross-rank int64 sum for the decomposition-invariant chksum bitcount (rdb_ocean_chksum). pic_mpi_lib exposes no integer(int64) allreduce overload (MPI_INTEGER8 reaches only send/recv), so the value rides the EXACT real64 allreduce: a per-field POPCNT sum is bounded by (#elements x 64), which stays FAR below the double-mantissa bound 253 for any realistic grid (253/64 ~ 1.4e14 cells), so both the local->double cast and every partial MPI_SUM are exact — the invariant survives.

Free pre-allocated halo buffers

Pre-allocate halo buffers for all 4 fields

Exchange ghost-cell halos for a single 2D field

GPU-direct halo exchange via CUDA-aware MPI

Exchange ghost-cell halos for a 3D field (all nz layers packed per direction)

GPU-direct halo exchange for a 3D field, batched across layers.

Pack and post non-blocking MPI sends/recvs for all 4 fields

Wait for MPI to complete and unpack received ghost cells

Release the persistent halo buffers. Idempotent.

Release the persistent 3D halo buffers. Idempotent.

Lazy-allocate the persistent send/recv buffers used by halo_exchange_2d_device (and halo_exchange_2d via the same pool). Sized to the per-rank subdomain; resizes on grid change.

Lazy-allocate the persistent 3D send/recv buffers used by halo_exchange_3d_device. Sized to the per-rank subdomain and the largest nz seen so far – subsequent calls with a smaller nz reuse the existing buffer (just write fewer elements), grow-only on nz so an ML run that later does an NH w halo (nz_ml+1) doesn’t free/reallocate. Resizes in full if (nx_total, ny_total, nghost) change.

Resolve an opaque c_ptr to a Fortran pointer and verify it.

Allocate a new ocean handle and return its opaque C pointer. Caller (rdb_ocean_create_from_string) is responsible for eventually passing it to handle_destroy — including on every early-return failure path, so a handle that never reached is_initialised does not leak.

Free an ocean handle. Idempotent on null/already-destroyed/garbage handles — silently no-ops rather than aborting, so a Python __del__ can call this blind. Caller is responsible for having already unwound device residency and called engine_teardown (engine%state%destroy() releases the god-state’s Fortran allocations) — h%state is only a pointer VIEW onto h%engine%state (see ocean_handle_t’s docstring), never separately allocated, so there is nothing to deallocate through it here; this only releases the handle’s own heap block (which takes engine — a value component — down with it). garbage or already destroyed

True when the handle’s metrics carry a full ghosted areaT (every grid the engine builds does); anything else falls back to grid%dx*grid%dy.

.true. iff the configured biharmonic dispatch will produce a STRICTLY POSITIVE dissipation coefficient somewhere. Drives the configure-time fail-loud guard requiring a biharmonic backstop when MEKE backscatter is on — the negative harmonic backscatter feeds a grid-scale mode that only a positive biharmonic can dissipate.

Henyey, Wright & Flatte (1986) JGR 91:8487 latitude dependence of the internal-wave-driven mixing rate, in the SIMPLIFIED constant-N0 form of Harrison & Hallberg (2008) JPO 38:1894 — the in-situ column stratification is replaced by a fixed reference N0, so the factor collapses to a pure function of latitude:

Corner thickness of the HK PV, recomputed with EXACTLY the Pass 2 formula of coriolis_adv_compute_tendencies_hk (wet-area-weighted 4-cell mean, array-edge clamps, CORIOLIS_H_MIN_PV floor), so the closed-face branch knows the h_corner each stored q was divided by without a persistent buffer.

One Arakawa-Hsu pair coefficient (q1 + q2 + q3)/12 with each PV re-evaluated at a corner thickness of at least h_ref/2: q → q·h_X/(h_ref/2) where h_X < h_ref/2, q unchanged otherwise (so an inactive floor returns the unfloored sum bit for bit). h_ref is the larger thickness of the pair’s u- and v-face.

One HOTFACE row (see chksum_hotface docstring for columns).

Add the Smith & McWilliams (2003) anisotropic direction-tensor coefficients onto the co-located isotropic viscosities. The tension (T-cell) coefficient gains kh_aniso·(1−n1n2²) and the shear (Bu-corner) coefficient gains kh_aniso·n1n2². For the default grid-i direction n1n2 = 0 ⇒ T gains kh_aniso, the corner gains nothing — stronger damping of along-i tension. The corner outer ring stays untouched (the divergence stencil never reads it; hvisc_avg_A_face already zeroed it).

Average the per-face harmonic viscosity (ah_face_x at u-faces, ah_face_y at v-faces) onto the T-cell centres (ah_t) and the Bu corners (ah_q). The stress form needs A co-located with the tension (T-cell) and shear (corner) strains; the lateral-mix closure produces A at faces, so this is a 4-point face→cell / face→corner reduction. On a uniform A field every average returns A, preserving the Laplacian reduction.

Pass 1 of BOTH velocity biharmonics (scalar nu_4 and the flow-aware nu4_face_*) under &vcoord_nml zfixed_closed_faces: the intermediate Laplacian lap_u/lap_v with a closed face-layer treated as a FREE-SLIP wall.

Per-cell CFL viscosity limiter (MOM6 BOUND_KH). Clamps the T-cell viscosity to Kh_Max_xx and the corner viscosity to Kh_Max_xy, each derived from the actual discrete stress stencil metrics + dt so the explicit forward-Euler viscous update can never overshoot. Replaces the global ah_max cap. On a uniform square grid Kh_Max = bound_coef·0.25/(dt·(1/dx²+ 1/dy²)).

Flow-aware biharmonic friction (MOM6 SMAGORINSKY_AH analogue). Identical to hvisc_compute_biharmonic_impl except Pass 2 multiplies the second Laplacian by the per-face viscosity nu4_face_x/y instead of the scalar nu_4. The face fields are filled upstream by ocean_lateral_mix_compute_smag_ah, which sets them to C_b · L⁴ · |D| clamped to [nu4_bg, nu4_max]. Pass 2 additionally clamps each face coefficient to the per-face explicit-biharmonic CFL ceiling (hvisc_nu4_cfl_bound, scaled by bound_coef) on top of the static nu4_max floor/ceiling applied upstream — so a strain spike on a fine cell can never violate the local CFL bound.

Constant-coefficient biharmonic friction: applies -ν₄ · ∇²(∇²u) to the face velocities via two chained 5-point Laplacians. Adds into the existing du_visc / dv_visc buffers (which Laplacian friction has already filled), so the caller can run with both nu_h and nu_4 non-zero.

Per-face metric Laplacian × spatially-varying viscosity. Explicit-shape dummies so NVHPC stdpar emits a device kernel without per-launch descriptor walks. See metric_lap_u/metric_lap_v for the curvilinear FV form. bound_kh engages the per-face harmonic CFL ceiling (hvisc_kh_cfl_bound); .false. ⇒ bit-identical.

Per-face metric Laplacian × scalar viscosity. Used when no lateral-mix closure is active — falls back to constant nu_h. When nu_h = 0 the kernel still zeros all interior + boundary cells so the apply step sees a defined state. bound_kh engages the per-face harmonic CFL ceiling (hvisc_kh_cfl_bound); .false. ⇒ bit-identical.

MOM6 thickness-weighted stress-divergence operator. Three phases: (1) tension str_xx at T-cells, (2) shear str_xy at Bu corners, (3) the divergence (1/(h_u+h_neglect))·∂str. wet_u/wet_v/wet_q mask the stress so no momentum is diffused across a coastline; h_neglect = H_VANISHED floors the velocity-point thickness. See module header for the form + the all-wet uniform-h Laplacian reduction.

Fill the T-cell and corner harmonic-viscosity fields with the scalar nu_h (no flow-aware closure active). Constant ⇒ the T/corner averaging is exact, so the all-wet uniform-h reduction to the velocity Laplacian holds bit-for-bit.

C-grid KE budget: each face’s u·du_visc rate is split half to each adjacent T-cell and depth-integrated with ρ_k h_k. Race-free — every (i,j) writes only its own ke_diss(i,j).

Per-face harmonic-viscosity ceiling for the velocity-Laplacian paths (MOM6 Kh_Max_xx analogue, uniform-grid reduction): ν_max = bound_coef · 0.125 / (dt · (idx² + idy²)) — one quarter of the forward-Euler stability limit ν·dt·4·(idx²+idy²) ≤ 2, the same margin MOM6’s harmonic bound uses (“avoid overshoots when bound_coef < 1”). Keeps the FROZEN depth-mean viscous forcing on the barotropic mode (F_bt) out of the phase-reversed anti-damping regime for grid-scale gravity modes (see the bound_kh docstring). Returns huge (no clamp) for a fully-masked face. !$acc routine seq — called from the Laplacian do concurrent.

Per-face explicit-biharmonic CFL ceiling on the biharmonic viscosity ν₄ (m⁴/s). Forward-Euler stability for −ν₄·∇⁴u on a local cell of spacing (dx, dy) requires (established project constant) ν₄ · dt · ((π/dx)² + (π/dy)²)² ≤ 2, so the per-face bound is ν₄_max = bound_coef · 2 / (dt · ((π·idx)² + (π·idy)²)²) where idx = 1/dx, idy = 1/dy are the metric inverses at the face (idxCu/idyCu at u-faces, idxCv/idyCv at v-faces). bound_coef (MOM6 HORVISC_BOUND_COEF, default 0.8) is the CFL safety margin shared with the harmonic hvisc_clamp_A. Returns a huge value (no clamp) when the metric inverses are both zero (a fully-masked land face) so the caller’s min leaves the coefficient untouched there. !$acc routine seq — called from the biharmonic do concurrent.

Surface-layer thickness of the hycom z* nominal floor (m), for the configure banner: the stretched profile’s top entry, else the uniform z_fixed_h_ref/nz, else 0 (the unconfigured sigma fallback of ocean_vcoord_rho_target_column).

Outer shim (outer-shim + flat-impl pattern). No-op when the ice slot is not live (is_init gate) or when ice%ncat == 1 (the legacy-mode bit-identity contract — see module docstring). No dt/eos needed: this is a pure area/mass reshuffle, no thermodynamics.

Device kernel: ONE do concurrent(j, i) over PHYSICAL cells (ghosts excluded), wet_mask > 0.5 inner gate (never a masked DC header), serial cat loops inside — each cell touches only its own (i, j, :) slice, race-free. No per-cell gather into local category-sized arrays (avoids register pressure and any compile-time ncat cap, unlike the ICE_NK_MAX-capped PER-LAYER arrays elsewhere in the ice model) — scalar temporaries only, all declared in local(...).

Outer shim: forward the four seam arrays + four scalars to the device kernel. No registry indirection here, but keep the shim+_impl split for the explicit-shape device-kernel discipline.

Full-array fill (incl. ghosts). Explicit-shape + decl-order (dims first). Runs on device-resident seam arrays (mapped by ice%enter_data).

Bottom-freezing branch of ice_resize_SIS2 (SIS2_ice_thm.F90:1164-1188; prototype sis2_resize.py:20-45). When bmelt < 0 (net upward heat deficit at the base), freeze ocean water onto the bottom ice layer: enth_freeze = min(enthalpy(nk), enth_ocean - min_denth_freeze), m_freeze = -bmelt/(enth_ocean - enth_freeze), mass-weighted mix of the bottom layer’s enthalpy + salinity, bmelt reset to 0.

Bottom melt peel (SIS2_ice_thm.F90:1246-1271; prototype sis2_resize.py:83-111). Same peel as ice_top_melt_peel but from k=nk down to k=0. The prototype’s separate ablation return is dropped (it is h2o_ice_to_ocn restricted to this call; PR 3b can re-derive it if needed).

SIS2 cell_ave_state_to_ice_state (:540). Per cell, per category: pre-floor category 1, optional thin-ice rolling (roll_factor > 0), general floor, then re-derive part_size(c) = mca_ice(c)/m_ice(c) and m_snow(c) = m_ice(c)*(mca_snow(c)/mca_ice(c)) (per-ICE-area). part_size(0) = 1 - Sum_c part_size(c) — MAY be negative here; ice_compress_impl (Phase 4) is what restores >= 0.

SIS2 zonal_mass_flux (SIS_continuity.F90:1064): PPM reconstruction of the category-SUMMED mass htot, ONE total face transport uhtot from the swept-volume parabola integral (volcfl_face, bit-for-bit the SIS2 face expression), then the PROPORTIONATE split uh(c) = uhtot*mca(donor,c)*I_htot(donor) (SIS_continuity.F90:1199-1205, Adcroft reciprocal — I_htot=0 when htot(donor)<=0). Stencil + edge STORAGE CONVENTION + swept-face orientation copied verbatim from continuity_zonal_flux (rdb_continuity.F90:871-948): hl_x_work(i) == h_face_left_x(i) is the value AT east face i from the LEFT cell i-1 (that cell’s OWN downwind edge); hr_x_work(i) == h_face_right_x(i) is from the RIGHT cell i (its OWN left edge). H3-style limited edges + CW84 + ppm_limit_pos (SIS2 runs PPM_limit_pos UNCONDITIONALLY on the PD scheme, so it is not optional here). CFL metric: SIS2’s shipped default vol_CFL=.false. uses CFL = |u|*dt*IdxT(donor) (SIS_continuity.F90:1165) — the T-cell inverse spacing, NOT the dy_cu*iareaT swept-area ratio (== SIS2’s vol_CFL=.true. variant). Identical on uniform Cartesian; correct on spherical / anisotropic.

Meridional twin of ice_cat_flux_x_impl. Same face-indexed edge STORAGE + swept orientation as continuity_meridional_flux (rdb_continuity.F90:1159-1202): hl_y_work(i,j) == north-face h_face_left_y (from the SOUTH cell j-1), hr_y_work(i,j) == h_face_right_y (from the NORTH cell j). CFL uses IdyT(donor) (SIS2 vol_CFL=.false. default).

Two-mode per-cell concentration/mass gather (module docstring convention), shared by rdb_ice_evp (ice_evp_step) and rdb_ice_ocean_coupler (ice_ocean_stress_flux) so neither module depends on the other.

Per-(cell,category) orchestrator: gather (bottom-up -> top-down flip, TRAP #2), optics, conduction (ice_temp_sis2), resize (snow add, bottom-freeze, top/bottom melt peel, rebalance), scatter (flip back).

KEEP IN SYNC with ice_transport_compress_cell (the HOST tested seam twin, directly below). Same excess/ratio/compaction algorithm; this twin exists only for a different argument shape — full device-present arrays + scalar (i,j) indices for the fused device path (fixed-size device locals), vs the twin’s small standalone per-cell arrays for the unit test. Any change to the compaction logic here MUST be mirrored there (the test only drives the twin).

Outer per-cell dispatch: ONE do concurrent(j,i) over physical cells (ghosts excluded), wet_mask > 0.5 inner gate, serial in category within the cell (same shape as ice_adjust_categories_impl). The per-cell algebra is the SAME algorithm ice_transport_compress_cell implements (public, directly unit-testable on small standalone arrays — SPEC §7 test 3’s single-cell hand-check) — but this production kernel does NOT call it with a derived slice: part_size(i,j,:) etc. are NON-CONTIGUOUS sections (the category axis is not the fastest dimension), which a device !$acc routine seq call cannot take safely (would force a compiler temporary — the array-of- derived-type/strided-slice indirection trap, CLAUDE.md memory). Instead the algorithm is INLINED here operating on (i,j,c) triples directly, mirroring ice_adjust_categories_impl’s established pattern. ok is a per-cell-then-reduced flag (D7: SIS2 FATALs on the top-category overflow inconsistency; ported as a !$acc parallel loop reduction instead, since do concurrent cannot itself carry a boolean/min reduction in the house style — CLAUDE.md “Reductions use !$acc parallel loop reduction(...)”).

Outer shim (outer-shim + flat-impl pattern): dereference the tracer registry (ms%tracers(idx)%hTr) on the HOST and forward bare arrays to the device kernel — same rule as ice_frazil_accumulate / ice_frazil_uptake. No-op when either S or T is unregistered.

Device kernel over PHYSICAL cells (ghosts excluded — same physical-cells-only contract as ice_frazil_accumulate_impl). fb/sst_seam/ssurf_seam/tfw_seam are zeroed unconditionally first. The SAMPLE seam (sst_seam/ssurf_seam/tfw_seam) is filled on every wet, non-vanished cell (harmless — the column reads it only where it has ice). But fb is filled ONLY where BOTH the cell is wet+non-vanished AND it carries ice (sum_cat m_ice > ICE_RHO_ICE*H_VANISHED, exactly ice_thermo_columns’ own per-cat ice threshold, summed): no ice base ⇒ no basal flux. This keeps the coupler and the column in lockstep on which cells exchange, so an ice-free warm ocean cell (SST > T_f, the normal open-ocean state) reports fb = 0 and the melt-side reduce kernel’s -fb term is a harmless subtraction of zero (rather than a spurious ocean-cooling Q_heat = -fb).

Ice specific enthalpy (J/kg) from temperature + bulk salinity — SIS2 enth_from_TS (SIS2_ice_thm.F90:1647). NOTE: here (and only here) the freezing point uses max(0, s), per SIS2. The SIS2 else branch at :1667-1671 (Cp_brine /= Cp_ice, the T_fr*log(T_fr/T) form) is deliberately NOT ported — with ICE_CP_BRINE == ICE_CP_ICE it is unreachable, and dropping it keeps the map closed-form-invertible (see module docstring).

Enthalpy (J/kg) of liquid water at temperature t — SIS2 enthalpy_liquid (SIS2_ice_thm.F90:1691). s is unused in this linear form; the argument is kept for SIS2 call-site parity (enthalpy_liquid(T, S, ITV)).

Enthalpy (J/kg) of liquid water at the freezing point for salinity s — SIS2 enthalpy_liquid_freeze (SIS2_ice_thm.F90:1679).

One outer (slow) EVP call: evp_sub_steps subcycles advancing ui/vi/str_d/str_t/str_s, plus the subcycle-averaged ice->ocean stress fxoc/fyoc.

Flat-impl core of ice_evp_dynamics: the EVP subcycle body with the persistent scratch passed as EXPLICIT-SHAPE dummies (memory: never assumed-shape into a do concurrent feeder — NVHPC would emit descriptor-walk memcpys per launch). The *_w scratch names mirror the retired module workspace 1:1, so the body below is unchanged from the pre-slot version.

mi_ratio_A_q at a single corner (SIS2 :926-964), FULL form — all four branches (interior / corner-coast / straight-coast / land). weak_coast_stress=.false. hardwired (SIS2 default): sum_area is the MASKED area sum of the 4 surrounding T-cells. Factored out of the fill kernel so a unit test can pin it directly (SPEC §7 gate 8).

Small constructor — build once from &ocean_ice_nml config. 11 args (> the style guide’s 6): pre-existing deviation, sanctioned by the derived-type-grouping escape hatch (FORTRAN_STYLE.md §Public procedure arguments) — the whole point of ice_evp_params_t is to be this constructor’s one-shot host. Positional (not optional): a knob threaded to the config/schema but not to this constructor reads from the namelist, validates, and does nothing — the dead-knob class the audit indicts. The compiler catches the omission; optional would not.

Gathers mis/mice/ci from the category state (mode-branched, mirrors PR 4b’s IST->CAS dispatch), pulls the one-step-lagged ocean surface velocity, and calls ice_evp_dynamics on ice%u_ice/v_ice/str_d/str_t/str_s/fxoc/fyoc. No-op when the ice slot is not live or dynamics is off (defence-in-depth; the driver already gates this call on ice%dynamics).

Fill fb_part_sum(i,j) = Σ_c part_size(i,j,c) restricted to categories passing the column’s OWN entry gate (m_ice(i,j,c) > ICE_RHO_ICE*H_VANISHED) — MUST run BEFORE ice_thermo_columns mutates m_ice (module docstring). Device kernel over PHYSICAL cells, inner if/serial do cat — never a masked do concurrent header.

Outer shim (outer-shim + flat-impl pattern): dereference the tracer registry (ms%tracers(idx)%hTr) on the HOST and forward bare arrays to the device kernel — NVHPC stdpar cannot follow the array-of-derived-types indirection inside a do-concurrent body. No-op when either S or T is unregistered.

Device kernel over PHYSICAL surface cells (ghosts excluded — wall ghosts are inert and seam ghosts are rebuilt by the wrap / exchange, so clamping them would double-count the bank in any area integral). Per wet cell at k = nz:

Outer shim (outer-shim + flat-impl pattern): dereference the tracer registry (ms%tracers(idx)%hTr) on the HOST and forward bare arrays to the device kernel — NVHPC stdpar cannot follow the array-of-derived-types indirection inside a do-concurrent body (same rule as ice_frazil_accumulate). No-op when either S or T is unregistered.

Per-(cell,category-1) column worker: gather + flip (bottom-up -> top-down, TRAP #2 discipline), spend the whole frazil bank evenly over the nk layers (SIS2:1402), rebalance, scatter + flip back. Port of add_frazil_SIS2 (SIS2_ice_thm.F90:1342-1444), bulk-salinity mode (salin_freeze = ICE_BULK_SALINITY, SIS_slow_thermo.F90: 1207-1209) — the ice_rel_salin mode is NOT ported.

Device kernel over PHYSICAL cells (ghosts excluded — same physical-cells-only contract as ice_frazil_accumulate_impl and ice_thermo_columns). Per wet, non-vanished, banked cell: sample SST/SSS at k = nz, compute the seawater freezing point, spend the WHOLE bank on category-1’s column via ice_frazil_uptake_column, and reset frazil_heat to 0 (fully spent). The per-cell diags (m_frozen_diag, salt_flux_diag) are zeroed UNCONDITIONALLY first, then overwritten under the gate — an unbanked/dry/land/vanished cell reports zero, not a stale value from a prior window.

ncat>1 SIS2 ITD-mode frazil spend. Port of SIS2 SIS_slow_thermo.F90:1121-1146 + 1181-1186 (non-filling mode — SIS2’s default SIS2_FILLING_FRAZIL=.true. thin-category fill is DEFERRED, see module docstring). Per banked cell (same wet/non-vanished/bank>0 gate as ice_frazil_uptake_impl): 1. k_merge scan (SIS2:1124-1129): first category c with part(0) + part(c) > 0.01; falls back to k_merge = 1 if no category qualifies (SIS2’s k_merge default). 2. Open-water annexation (SIS2:1131-1145): if part(0) > 0, dilute category k_merge’s thickness at CONSTANT MASS — m_ice(k_merge)/m_snow(k_merge) scale by part(k_merge)/(part(k_merge)+part(0)), part(k_merge) absorbs all of part(0), part(0) -> 0. enth/sal are per-MASS intensive — untouched by an area-only dilution. 3. Per-ice-area spend (SIS2:1181-1186): frazil_col = frazil_heat/part(k_merge) (J per m² of category area; the denominator is > 0 by construction — either an occupied category was found, or step 2 just grew part(k_merge) from the Σpart=1 invariant), then the SAME UNCHANGED ice_frazil_uptake_column spends it on category k_merge’s column. 4. Diags (per cell, part-weighted back to CELL-area units to match the ncat==1 diag convention that the couplers and rdb_ice_thermo_driver consume): m_frozen_diag = part(k_merge)*m_frozen_pt; salt_flux_diag = part(k_merge)*(m_frozen_pt*s_surf - salt_to_ice_pt)/dt_therm; frazil_heat reset to 0 (fully spent). Both diags are zeroed UNCONDITIONALLY at loop top, same contract as ice_frazil_uptake_impl.

Gather pass (race-free): tr_flux_x_work(I,c) = val(donor of I,c), donor by the sign of the FLUX. uh is exactly 0.0 at the array-edge faces I=1 and I=nx+1 (the flux kernel zeros them unconditionally), so the I=1 donor-by-sign branch (uh>=0) would read the out-of-bounds val(0,j,c) — guarded explicitly below (I==1/I==nx+1 fall back to the IN-BOUNDS neighbour; the value is never actually consumed by the update kernel there since uh==0 at both those faces makes them inert, but the gather must still avoid the invalid index). Reads val ONLY at a donor cell — never writes val — so this kernel has no race with any other iteration.

Layer-indexed twin of ice_gather_flux_x_impl for enth_ice/sal_ice/enth_snow (shape (nx,ny,ncat,nk)).

Meridional twin of ice_gather_flux_x_impl.