Modules

Depth-integrated state on an Arakawa C-grid for the ocean dynamical core. C-grid layout: * Scalars (h, b, flux divergence) at cell centres, (nx_total, ny_total). * x-velocity/momentum at east faces, (nx_total+1, ny_total). * y-velocity/momentum at north faces, (nx_total, ny_total+1). Lives at ocean_state%barotropic for sim_type='ocean'.

Transient barotropic state the split driver reads/writes between the slow baroclinic step and the fast barotropic substep loop. Transient (not restartable), unlike the prognostic barotropic_state_t. C-grid stagger: scalars at centres, u at east faces, v at north faces, ζ at corners. F_slow_*, F_bt_*, ubt_at_n, vbt_at_n, etc. are allocated only when init is passed nz_ml (barotropic-only unit tests skip them).

find_uhbt/find_duhbt_du are pure + !$acc routine seq (callable from the BT substep do concurrent); uhbt_to_ubt keeps a Newton iteration and is host/single-thread only (OBC setup).

Reads simulation parameters from a namelist input file

The shared print layout for the periodic conservation + stability console line, used by BOTH the ocean and coastal drivers so the two regimes report identically. Each regime computes its own totals (area-weighted device reductions over its own state), then hands the scalars here; this module owns the reference-snapshot latching, the relative-drift arithmetic, the exact line format, and the NaN / CFL panic guards.

Holds scheme-variant flags and reusable workspace for the C-grid continuity-PPM (Lin & Rood) thickness-flux kernel, plus the public step routines that the split-explicit driver calls per barotropic substep. The kernel is the primary mass-flux producer for the ocean path: it consumes face velocities + cell-centred thickness and emits per-face mass fluxes that the tracer-advection kernels then reuse under CWC for free.

Pure Fortran, no MPI dependency. Each rank calls decomp_init with its process-grid coordinates to learn its local subdomain size and offsets.

Order-invariant real-to-integer summation (Hallberg & Adcroft, 2014, An Order-invariant Real-to-Integer Conversion Sum, Parallel Computing 40(5-6), doi:10.1016/j.parco.2014.04.007).

P0/P0.1 gave every setup procedure an optional, intent(out) :: ierr returning a coarse STAGE code (rdb_ocean_status) while the specific reason went only to global_logger%error — reachable from a Fortran log, not from a C-ABI caller. RuntimeError: error 3 is exactly the failure this module exists to prevent.

Exchanges ghost-cell strips between neighbouring MPI ranks

Fills ocean_sea_ice_t’s coupleable seam (atm_sf0/atm_dsfdt/atm_sw_dn/atm_fprec) from the &ocean_ice_nml scalar restoring knobs. SF(T) = lambda(T - T_air) => sf_0 = -lambdaT_air (upward-positive SEB intercept), dsf_dt = lambda (thermostat slope), sw_dn = sw_down (uniform), atm_fprec = snowfall (uniform, PR 26). The seam is identical whether this filler, a bulk formula, or a live coupler fills it – that is what makes the model coupleable (PLAN_ICE_PR3c S”coupleable seam”; PLAN_PR26_snowfall.md S13 fixes atm_fprec’s shape/units/lifecycle for PR-55 to fill from data).

Port of ice_temp_SIS2 + laytemp_SIS2 + update_lay_enth (SIS2_ice_thm.F90:169-945, Apache-2.0) under the ICE_CP_BRINE == ICE_CP_ICE simplification (rdb_ice_enthalpy module docstring): every per-layer implicit solve and T<->E inversion is a closed-form quadratic — no Newton / false-position iteration anywhere. Ported line-by-line from the validated stdlib-only Python prototype tmp_local_artifacts/ice_pr3a_prototype/ sis2_column.py; when any formula here and SIS2 itself seem to disagree, the prototype (which reproduces SIS2’s own commented-out col_check energy-closure diagnostic to ~1e-14 fractional) is the tiebreaker.

Mechanical transliteration of the validated Python prototype tmp_local_artifacts/proto_evp_core.py (Channel1D) + tmp_local_artifacts/proto_evp_validate.py (Box2D), grounded against SIS_dyn_cgrid.F90 (SIS_C_dynamics :603-1614, limit_stresses :1619-1741, SIS_C_dyn_init :209-270) per SPEC_ice-pr5-evp.md. Field names mirror SIS2 1:1 (str_d, str_t, str_s, sh_Dd, sh_Dt, sh_Ds, zeta, del_sh, mi_ratio_A_q, Tdamp, EC, …).

Port of add_frazil_SIS2 (SIS2_ice_thm.F90:1342-1444, Apache-2.0), bulk-salinity mode only. Spends the WHOLE per-cell frazil bank (ocean_sea_ice_t%frazil_heat, banked by rdb_ice_frazil’s surface-freezing clamp) once per thermo step by depositing new ice mass, evenly split across the nk_ice layers of category 1 (cat = 1 — the only category this v1 model fills; part_size is untouched, same convention as PR 3a).

Port of SIS2 adjust_ice_categories (SIS_transport.F90:611-891, Apache-2.0). SIS2’s shipped algorithm is NOT the Lipscomb (2001) linear-profile remap — that is a SIS_transport.F90:737 TODO comment, never code. What SIS2 actually runs (and what this module ports) is the “for now move all of it” WHOLE-CATEGORY shift: when a category’s mean ice mass crosses a bin boundary, its ENTIRE area+mass moves to the adjacent category, merged into that category’s existing area+mass via an area-weighted destination- thickness average (mass AND area conserving, SIS2:730-741) with upwind mass-weighted tracer mixing (equivalent in effect to SIS2’s deferred advect_tracers_thicker, SIS_tracer_advect.F90:1791-1850 — each category boundary is visited once per pass in a fixed order with running masses, so a single closed-form mass-weighted mix at the transfer site reproduces the same result; verified by the validated Python prototype, tmp_local_artifacts/proto_itd_adjust.py).

Bridges ocean_sea_ice_t%salt_flux_diag (filled by rdb_ice_frazil_uptake%ice_frazil_uptake) into ocean_surface_flux_t%Q_salt — the field the REAL production kernel ocean_surface_flux_apply_tracers (rdb_ocean_surface_flux.F90) reads every thermo step.

Complements rdb_ice_mass%ice_snow_accumulate (the on-ice snow add): the share of &ocean_ice_nml snowfall that lands where there is no ice (open water at ncat>1, ice-free cells at ncat==1, and any category that fails the column’s own entry gate, m_ice(i,j,c) > ICE_RHO_ICE*H_VANISHED) does not accumulate — it melts on contact with seawater. This module delivers that mass’s latent heat + virtual freshening to the ocean through the EXISTING heat_flux_diag/salt_flux_diag contributor seam (rdb_ice_thermo_driver.F90’s ordering contract), the same one ice_frazil_uptake and the melt-side reduce kernels already use. Without this, snowfall > 0 over open water would silently annihilate mass and energy (the sw_thru failure mode, PLAN_PR26_snowfall.md §2 / roadmap §5 trap 15).

Bridges ice_thermo_columns (PR 3a) into the ocean-coupling seam: runs the column, then reduces h2o_ocn_to_ice/h2o_ice_to_ocn/heat_to_ocn over categories into the per-cell heat_flux_diag/m_melt_diag and the net-melt salt contribution added into salt_flux_diag. Q_heat = sum_cat(heat_to_ocn)/dt_therm - fb [W/m^2, +down] m_net = sum_cat(h2o_ocn_to_ice - h2o_ice_to_ocn) [kg/m^2, +freeze] salt += m_net*(s_surf - ICE_BULK_SALINITY)/dt_therm (ADD, not set) sw_thru_diag = sum_cat(sw_thru) [W/m^2, +down] PR 31: sw_thru (the shortwave penetrating the ice to the water) is now reduced to the per-cell sw_thru_diag and coupled to the ocean by ice_ocean_sw_flux (rdb_ice_ocean_coupler). It is DELIBERATELY kept OUT of the heat_flux_diag sum above — heat_flux_diag carries only the NON-shortwave heat; the coupler delivers the shortwave via a separate path (the q_sw component, or a direct Q_heat add) so the energy reaches Q_heat exactly once in either component mode.

Provides NetCDF read/write primitives and an error-checking helper. Used by bathymetry loader and output writer.

Free procedures that estimate, log, and verify the device memory footprint of the simulation state before the OpenACC enter_data mapping is attempted — so an over-budget run reports a readable diagnostic instead of dying opaquely inside the CUDA runtime (the motivating failure: a 10M-cell ocean config that exhausted device memory at enter_data with no setup-time hint).

A host-only, schema-driven namelist parser and validator. Validates a namelist file against a registered schema of groups and keys, collecting ALL errors (with path:line: prefixes, did-you-mean suggestions, range/enum/length checks) and rejecting unsupported syntax (repeat-counts, indexed assignment, derived-type refs).

Oceananigans-style geometry construction (docs/ocean_python_api_plan.md S5b, 06_python_surface_design.md D6.2) hands bathymetry to the library as an interior-sized array instead of a topo_config formula or a NetCDF file. Two pieces live here, deliberately in a NetCDF-FREE module (unlike rdb_bathymetry, which requires RDB_ENABLE_NETCDF=ON purely because it also happens to contain the NetCDF reader): sign normalisation/validation, and ghost fill. Array injection is explicitly meant to need no filesystem/NetCDF at all — a caller building a grid in memory and handing it straight to create() must not be forced onto a NetCDF build just to get its ghosts filled.

Type taxonomy + per-edge BC config + the composed slot on ocean_state_t. Six BC types implemented end-to-end (WALL, OPEN, TIDAL, CLAMPED, SPONGE, CHAPMAN); INFLOW/DISCHARGE/NESTED tags are declared for cross-backend alignment but error stop if encountered. Per-edge granularity: each of the four outer edges carries one ocean_bc_face_tag_t, read independently by the dispatch helpers in rdb_ocean_boundary.

Shadow state for the wide-halo barotropic march-in. When bt_halo > 0, the BT fast loop runs on WIDE arrays with ghost width ng_wide = nghost + bt_halo. Ghost cells outside the valid band evolve stale data that creeps INWARD at 2 cells/substep; one grouped exchange every bt_halo/2 substeps keeps the physical interior clean. The mid-substep u exchange is absorbed (within the 2-cell/substep stencil budget).

Static ice-shelf cavity GEOMETRY: the prescribed ice draft z_draft(i,j) and the arithmetic that turns it into a water column, a grounding decision and an isostatic load.

Periodic conservation + stability scalars printed to stdout at the driver’s status_interval cadence (separate from the diag manager’s NetCDF). Mirrors MOM6’s “MOM Day N:” status line:

Generic, decomposition-aware, device-resident reader for time-varying NetCDF fields field(x, y[, z], t). Owns the file handles, the time axis, the bracket bookkeeping, the H->D motion (only when the bracket advances), and the per-step linear-in-time blend kernel. Owns NO field names and NO state-slot mapping: a consumer registers its own (file, variable, destination) triple via ocean_data_input_register_2d/_3d (or the edge-segment variants for OBC files) and gets back an opaque id.

Driver-level state + step routines for the outer (baroclinic) + inner (barotropic) split-explicit RK2 integrator of the ocean dynamical core. Carries the substep ratio n_inner, the fast-mode time-averaging accumulators, and energy/CFL diagnostic slots.

Hosts ocean_eos_compute, the outer shim that pulls the registered S, T tracer arrays off multilayer_state_t and forwards them as bare 3D arrays to the regime-agnostic eos_compute_arrays dispatch in rdb_eos.

Discrete tripolar north-fold exchange (Murray 1996), LOCAL to the tile that holds the whole fold row. Pure seam operators only (state orchestration lives in rdb_ocean_fold_apply). Free procedures, explicit-shape dummies, j-outer / i-inner do concurrent.

Owner-routed north-fold exchange for an east-west split north rank row, plus the px-dispatching single-field fold entry points (ocean_fold_north_*): px = 1 keeps the local kernels of rdb_ocean_fold verbatim, px > 1 always uses the exchange (the self-mirror tile through a local copy).

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

Roundabout has (at the time of writing) six independent minimum-nghost rules spread across five files, each enforced at its own call site in rdb_config.F90/rdb_ocean_boundary_types.F90/rdb_ocean_halo.F90/ rdb_ocean_setup.F90. A Python caller building a grid in memory (docs/ocean_python_api_plan.md S5b) needs to be able to ask “how wide must my halo be” WITHOUT reproducing that scatter — the model Oceananigans uses is a per-scheme trait folded by max() over every tendency term (inflate_halo_size, automatic_halo_sizing.jl:73-81; 06_python_surface_design.md B2.4). required_halo below is that fold, reusing the EXISTING per-rule functions (pv_adv_required_nghost, tracer_recon_required_nghost) rather than re-deriving their numbers, so this module and the enforcement sites it does not replace can never silently disagree.

Carries the closure parameters and per-step tendency workspace for the Laplacian horizontal-momentum-viscosity kernel. Sits alongside rdb_ocean_lateral_mix (which holds the variable-coefficient state for Leith/Smag closures): this module is the kernel, that one is the closure. Phase Tier-1 ships the constant-nu_h variant; later phases will read coefficients out of ocean_lateral_mix_t%ah_face_* and replace the scalar.

Neutral-density slope S = -∇ρ/∂_zρ and interface stratification N² at C-grid layer interfaces (u-faces → slope_x/n2_u, v-faces → slope_y/n2_v). Purely diagnostic — consumed by GM/Redi/VarMix/ MLE; no flux consumer here. Harmonic-thickness-weighted FD form (Griffies 1998).

Rebuild of the transient per-term KE attribution meter used to localise the split-corrector anti-damping term (LAGRANGIAN_PGF_BUG.md §R). Samples the h-weighted layer kinetic energy between the split stage’s tendency applies and prints KE_ATTR rows: the dKE between consecutive samples IS the energy injected/removed by the segment just executed, so one instrumented run names the anti-damping term directly.

Mesoscale eddy kinetic energy parameterization. Carries a single 2D vertically-averaged eddy-energy field meke(nx,ny) [m^2/s^2], sourced by the GM potential-energy release, damped by an implicit (backward-Euler) bottom drag, transported laterally (harmonic-mass Laplacian + optional biharmonic + advection), and fed back as a thickness/tracer diffusivity kh = khcoeff·sqrt(2·gamma_t²·E)·Lmix added (geom mean) into VarMix’s per-face KhTh before GM’s CFL clamp — closing the GM↔eddy-energy loop. Updated via a Strang split each thermo step.

Conservative minimum-layer-thickness adjustment for the isopycnal (VCOORD_LAGRANGIAN) ocean path.

Parent-state ingestion + interpolation state for nesting a regional ocean run inside a global/basin parent. Distinct from coastal BC_NESTED: parent state on a coarser mesh (ingest + remap onto our vertical grid), tides composed in via rdb_ocean_tides, mandatory T/S clamp, and a Flow-Relaxation-Scheme (FRS) zone of nrelax cells rather than a hard clamp. Scaffold (kernels not yet implemented).

Surface-pressure loading for the ocean C-grid dyn-core (PR-17). An atmospheric surface pressure p_surf(x,y) (Pa) adds a depth-uniform acceleration -(1/rho0) grad(p_surf) to the horizontal momentum (Wunsch & Stammer 1997). Because the term is depth-independent it is a purely BAROTROPIC forcing: it drives the free surface and, at rest, integrates to the equilibrium inverse-barometer response eta_ib = -p_surf/(rho0 g) + const (~1 cm depression per hPa of atmospheric high; Ponte 2006 for the closed-domain mean removal).

The 5-point Boole-quadrature density integrals, through the generic eos_t handle, that feed the finite-volume pressure-gradient force (FV_MOM6 variant of rdb_ocean_pressure_force): the in-layer vertical rule over a PCM / PLM / PPM sub-layer T/S profile (boole_dpa_intz_layer) and the cross-face rules (boole_dpa_face, boole_dpa_face_pcm). They are the REFERENCE the per-EOS fast paths are tested against, and the rule the kernel runs for an EOS without a twin (the linear EOS).

Carries scheme-variant flags + reusable workspace for the layered hydrostatic pressure-force kernels on C-grid face arrays: Montgomery potential, finite-volume (lite / Wright / MOM6) and the NK=2 reduced-gravity form.

Redi neutral (along-isopycnal) tracer diffusion. Continuous (non-iterative) neutral-surface geometry of the Griffies rotated-diffusion tensor + the two-phase GPU flux kernel. Clean-room from Redi (1982), Griffies et al. (1998), Griffies (2004).

Restart cadence + per-slot checkpoint registry for the ocean dyn-core. Each prognostic-owning slot registers its arrays at birth; the manager walks the registry to do the I/O, so new state-carrying slots join by registering rather than editing here.

Applies the per-concern &ocean_* namelist knobs onto an ocean_state_t’s kernel slots (bottom drag, vmix, lateral viscosity, vertical coord, PGF, barotropic split). Pure cfg → slot wiring + rank-0 logging — no I/O, no NetCDF — so it is shared by the production driver and the NetCDF-free benchmark.

Motivating failure (tmp_local_artifacts/global_run/FINDINGS.md, 2026-09-11): a global tripolar aquaplanet NaN’d at outer step 7. The ONLY diagnostic on offer was a bare non-finite-face count — “producer 0/0 upstream, investigate”. A human needed several runs + a bisection to find the actual cause: nu_h = 2.0e4 with dt = 900 s at the ~3 km polar cells gives a viscous-diffusion number of 1.65 against an explicit-Laplacian bound of 0.125 — 13x over — and bound_kh (the per-cell runtime clamp that would have protected against exactly this) was never enabled.

A library cannot abort its host process. Every procedure on the solver-creation path (config parse/validate, the configure_ocean_* setup chain, IC seeding) that used to error stop now takes an optional, intent(out) :: ierr — present ⇒ return one of these codes instead of aborting; absent ⇒ unchanged legacy error stop behaviour (every existing caller keeps aborting exactly as it does today).

Clean-room astronomy for the equilibrium (astronomical) body-force tide. Computes the four mean longitudes (moon s, sun h, lunar perigee p, ascending node N) from Schureman’s polynomials, the per-constituent equilibrium argument V_c (with the load-bearing ±pi/2 diurnal signs), and the slowly-varying nodal amplitude/phase corrections f_c(N), u_c(N). Host-side scalar generator — called once at init (to bake the ref-date offset into phase0) and once per outer step (to advance the running argument); NOT a device kernel.

Holds the vertical-coordinate configuration and the per-step target grid that the ALE remap step relamps multilayer.h_layer + every multilayer.tracers(t)%hTr onto. The coastal path has VCOORD_SIGMA, VCOORD_ZSIGMA, VCOORD_ZSTAR, VCOORD_ZSTAR_SIGMA, VCOORD_ZSTAR_FULL plus rdb_remap_column; this module is the C-grid counterpart for the ocean dynamical core.

Per-column first-baroclinic internal gravity-wave speed cg1 (m/s) and first-mode Rossby radius Rd (m), plus Rd/dx (the GM/Redi/MEKE resolution ratio). Solves the rigid-lid Sturm-Liouville eigenproblem discretised from layer thicknesses and per-interface reduced gravities gprime = (g/rho0)*max(0,drho) (rho_layer authoritative); largest c^2 via a fixed-budget Sturm-count bisection (no early exit -> warp-divergence-free). Reference: Chelton et al. (1998).

Pure per-column conservative vertical remap, shared by both backends. Reconstruction methods: PCM — piecewise constant (donor cell) PLM — piecewise linear, minmod limiter PPM — piecewise parabolic (Colella & Woodward 1984) PPM_H4 — PPM with non-uniform 4th-order edge values (White & Adcroft 2008) PQM — piecewise quartic (White & Adcroft 2008)

Single owner for (n1, n2, n3) scratch arrays. Every column- local kernel slot (continuity, coriolis_adv, pressure_force, lateral-mix vorticity scratch, KPP solve scratch) declares its per-step workspace as a scratch_3d_buffer_t rather than as a bare real(wp), allocatable :: foo(:, :, :).

Per-tracer state and config for the layered ocean solver.