ocean_dyn_step_split Subroutine

   subroutine ocean_dyn_step_split(grid, metrics, dyn, eos, cor, ct, pgf, hv, bd, ss, &
                                   va, hd, vd, vmix, ms, dt, n_inner, sf, geo, vcoord, bc, sp, t, &
                                   lateral_mix, epbl, kshear, mle, slopes, gm, varmix, wavespeed, &
                                   redi, meke, vmix_tidal, tides, psurf, td, cav)
      !! Split-explicit SSP-RK2 outer step on the multilayer state.
      !! Parallel to `ocean_dyn_step` (the unsplit driver
      !! still ships for tests + reference).  Phase 4b-MVP scope:
      !! gravity-wave-stability fix on momentum only.  ALE-aware
      !! `h_layer` redistribution and full nonlinear-bt corrections
      !! are deferred to a follow-up branch.
      !!
      !! Algorithm per SSP-RK2 stage:
      !!   1. Derive (η^n, u_bt^n, v_bt^n) from current `ms` and stash
      !!      u_bt^n, v_bt^n in `dyn%bt_work%ubt_at_n` / `dyn%bt_work%vbt_at_n`.
      !!   2. Run all slow computes (writes the per-kernel scratch
      !!      buffers).
      !!   3. Sum the per-face slow tendency from those scratches
      !!      into `dyn%bt_work%F_slow_u` / `dyn%bt_work%F_slow_v`.
      !!   4. Depth-mean → `dyn%bt_work%F_bt_u` / `dyn%bt_work%F_bt_v`.  (Weighted by
      !!      face thickness.)
      !!   5. Existing applies for momentum + continuity (uses the
      !!      same per-kernel scratches — produces u^*, v^*, h^* the
      !!      unsplit driver would).
      !!   6. Run the barotropic substep with F_bt as constant forcing for
      !!      `n_inner` substeps at dt_inner = dt / n_inner — fills
      !!      `dyn%bt_work%bt_eta` / `dyn%bt_work%bt_ubt` / `dyn%bt_work%bt_vbt` with the
      !!      time-mean.
      !!   7. Correction step: add (⟨u_bt⟩ - u_bt^n - dt·F_bt_u) to
      !!      every layer's u, v.  This replaces the unsplit bt mode
      !!      (which is FE-amplified on the gravity wave) with the
      !!      barotropic-substep's resolved bt mode (stable under FBE for
      !!      dt_inner·c·k < 2).
      !!
      !! Two RK2 stages, then RK2 average.  Per-tracer save/average
      !! is identical to the unsplit driver.
      !!
      !! `dyn%bt_work%bt_H_ref` must be set BEFORE the first split step (the
      !! caller initialises it from the time-mean total H of the
      !! initial multilayer state; for Eulerian-z with a flat
      !! bathymetry it's just the total column depth).
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
         !! Curvilinear horizontal metrics — forwarded to the PGF (and,
         !! in slice 2, the other geometry-aware kernels).
      type(ocean_dyn_t), intent(inout) :: dyn
      type(eos_t), intent(in) :: eos
      type(coriolis_adv_t), intent(inout) :: cor
      type(continuity_t), intent(inout) :: ct
      type(ocean_pressure_force_t), intent(inout) :: pgf
      type(ocean_horizontal_viscosity_t), intent(inout) :: hv
      type(ocean_bottom_drag_t), intent(inout) :: bd
      type(ocean_top_drag_t), intent(inout), optional :: td
         !! Ice-shelf TOP-drag slot (`&ocean_tdrag_nml`).  OPTIONAL so the
         !! many direct `ocean_dyn_step*` / `run_stage*` call sites in the
         !! test suite need no churn; the production driver always passes
         !! it.  Absent, or present and disabled, => no kernel launch and
         !! a bit-identical step.
      type(ocean_cavity_flux_t), intent(inout), optional :: cav
         !! Ice-shelf basal-melt slot (`&ocean_cavity_melt_nml`).
         !! OPTIONAL for the same reason `td` is: the direct
         !! `ocean_dyn_step*` / `run_stage*` call sites in the test suite
         !! need no churn, and the production driver always passes it.
         !! Absent, disabled, or `freshwater="virtual"` => no kernel
         !! launch and a bit-identical step.  It is threaded down here
         !! rather than acted on in `engine_step_finalize` because the
         !! real-freshwater volume must be spent in the SAME stage, at
         !! the SAME stage weight and from the SAME `melt` value as the
         !! salt and heat halves the surface-flux apply spends -- see
         !! `ocean_cavity_mass_step`'s docstring.
      type(ocean_surface_stress_t), intent(inout) :: ss
      type(ocean_vertical_advection_t), intent(inout) :: va
      type(ocean_hdiff_tracer_t), intent(inout) :: hd
      type(ocean_vdiff_t), intent(inout) :: vd
      type(ocean_vmix_t), intent(inout) :: vmix
      type(multilayer_state_t), intent(inout) :: ms
      real(wp), intent(in) :: dt
      integer, intent(in) :: n_inner
      type(ocean_surface_flux_t), intent(in), optional :: sf
      type(ocean_geothermal_t), intent(in), optional :: geo
         !! Geothermal bottom-heat-flux slot.  Absent or `enable=.false.`
         !! preserves the historical no-geothermal path bit-identically.
      type(ocean_vcoord_t), intent(inout), optional :: vcoord
         !! Vertical-coordinate state.  When present and
         !! `vcoord%coord_type /= VCOORD_EULERIAN_Z`, the orchestrator
         !! runs the ALE remap step after the RK2 average (Lagrangian-
         !! then-remap pattern, MOM6-style).  Absent or EULERIAN_Z
         !! preserves the pre-existing Eulerian-z behaviour bit-
         !! identically.
      type(ocean_bc_state_t), intent(inout), optional :: bc
         !! Open-boundary config.  Absent or all-OBC_WALL preserves
         !! the closed-wall behaviour bit-identically (the barotropic substep's
         !! tag dispatch falls through to hard-zero).  When any edge
         !! is non-WALL, the corresponding BC variant fires inside
         !! the barotropic substep's per-substep wall closure.
      type(ocean_sponge_t), intent(in), optional :: sp
         !! Map-driven sponge slot (PR-23).  Absent or `enable=.false.`
         !! (default) preserves the legacy `bc`-band sponge path bit-
         !! identically; `enable=.true.` supersedes it (exactly one of the
         !! two runs — see `run_stage_split`'s dispatch).
      real(wp), intent(in), optional :: t
         !! Wall-clock time at the start of this outer step (s),
         !! used to evaluate the OBC_TIDAL constituent table.
      type(ocean_lateral_mix_t), intent(inout), optional :: lateral_mix
         !! Flow-aware lateral-viscosity closure (Leith / Smagorinsky).
         !! Absent or `closure = LMIX_NONE` falls through to the
         !! scalar `nu_h` in `hv`, bit-identically.
      type(ocean_epbl_t), intent(inout), optional :: epbl
         !! Energetics-based PBL slot.  Absent or `enable=.false.`
         !! preserves the historical PP81/KPP path bit-identically.
      type(ocean_kappa_shear_t), intent(inout), optional :: kshear
         !! Kappa-shear interior closure slot.  Absent or
         !! `enable=.false.` preserves the historical path bit-identically.
      type(ocean_tidal_mixing_t), intent(inout), optional :: vmix_tidal
         !! Tidal-mixing interior closure slot.  Absent or
         !! `enable=.false.` preserves the historical path bit-identically.
      type(ocean_mle_t), intent(inout), optional :: mle
         !! Fox-Kemper MLE slot (B5).  Absent or `enable=.false.`
         !! preserves bit-identity.  Transports are computed once per
         !! outer step at THERMO cadence (from the prior step's
         !! `epbl%mld`) and folded into the continuity mass fluxes in
         !! both RK2 stages.
      type(ocean_slopes_t), intent(inout), optional :: slopes
         !! Isopycnal-slope diagnostics slot.  Refreshed at THERMO cadence
         !! (the split driver otherwise never calls `ocean_slopes_compute`)
         !! so the GM slot has a fresh slope to consume.  Absent or
         !! disabled ⇒ no-op.
      type(ocean_gm_t), intent(inout), optional :: gm
         !! Gent-McWilliams thickness-diffusion slot (capability [2]).
         !! Absent or `enable=.false.` preserves bit-identity.  Its slopes /
         !! VarMix / MEKE inputs refresh at THERMO cadence at the top of the
         !! step; the bolus transport itself is computed from the
         !! post-dynamics thickness and applied as its own sequential
         !! operator EVERY outer step after the stage loop (`run_gm_step`,
         !! MOM6 `thickness_diffuse` after `step_MOM_dyn_split_RK2`).
      type(ocean_varmix_t), intent(inout), optional :: varmix
         !! VarMix slot (capability [4]): spatially-varying GM/Redi
         !! coefficients.  When present + `enable=.true.` (and wavespeed
         !! present) `varmix_compute` fills the pre-CFL base `khth_u/v` at
         !! THERMO cadence BEFORE GM, and GM consumes them as its external
         !! base.  Absent or disabled ⇒ GM uses its scalar `khth` ⇒
         !! bit-identity.
      type(ocean_wave_speed_t), intent(inout), optional :: wavespeed
         !! Wave-speed slot supplying `cg1` to the VarMix resolution
         !! function.  Only read when `varmix` is active.
      type(ocean_redi_t), intent(inout), optional :: redi
         !! Redi continuous neutral-diffusion slot (capability [3]).  Absent
         !! or `enable=.false.` preserves bit-identity.  Phase-A neutral-
         !! surface coefficients are computed ONCE per outer step at THERMO
         !! cadence here (tracer-independent geometry); Phase B applies the
         !! rotated flux per tracer inside `run_stage_split` after the
         !! along-coordinate `tracer_hdiff`.
      type(ocean_meke_t), intent(inout), optional :: meke
         !! MEKE prognostic eddy-energy slot (capability [5]).  Stepped once
         !! per outer step at THERMO cadence right after `varmix_compute`:
         !! it reads `gm%gm_src` from the PREVIOUS outer step's GM operator
         !! (one-step lag; `gm_src` is restart-registered) and feeds the
         !! geom-mean of its derived `kh` into `varmix%khth_u/v` (+ khtr),
         !! which this step's GM operator CFL-clamps after the dynamics.
         !! Absent or `enable=.false.` ⇒ no-op (bit-identical).
      type(ocean_tides_t), intent(inout), optional :: tides
         !! Equilibrium body-force tide slot (C1) + scalar SAL (C2).  When
         !! present and `enable`, `eta_eq` is refreshed ONCE per outer step
         !! (held static across the inner substep loop); scalar SAL then
         !! folds `beta_sal*eta` (lagged barotropic SSH) into the combined
         !! `eta_forcing = eta_eq + eta_sal` forwarded to the barotropic
         !! PGF.  Absent / disabled / `use_sal=.false.` ⇒ bit-identical.
      type(ocean_p_surf_t), intent(inout), optional :: psurf
         !! Atmospheric surface-pressure loading / inverse barometer slot
         !! (PR-17).  When present and `enable`, the assembled surface
         !! pressure `sf%p_surf` (read only; Pa) is converted to
         !! `eta_ib = -p_surf/(rho0 g_bt)` ONCE per outer step and folded
         !! into `eta_seam = eta_ib [+ tide eta_forcing]`, which then feeds
         !! the barotropic PGF in place of the tide-only seam.  Requires the
         !! PR-12 component set (`sf%p_surf`) — guarded at configure.
         !! Absent / disabled ⇒ bit-identical.

      integer :: it, stage
      integer :: i, j, nx_ptop, ny_ptop
         !! Loop indices + extents for the E3 `ms%p_top` refresh below.
      logical :: tide_on, psurf_on, p_top_live
      real(wp) :: t_now
         !! Model time (s) for `ocean_ideal_age_young_val`; `t` fallback (PR-7).
      integer :: nan_i, nan_j, nan_k
      logical :: nan_is_u
      real(wp) :: nan_dx, nan_visc_cfl
         !! Location + local diagnostics for the FIRST non-finite face the
         !! velocity-truncation NaN-catch finds this outer step (see the
         !! "[nan-catch] outer step" report below) — always requested from
         !! `apply_velocity_truncation` (cheap: only actually searched
         !! inside that call when `n_nan>0`).

      ! ---- Ghost-band poison (debug knob) ----
      ! Sentinel-NaN every exchange-covered ghost band BEFORE any exchange
      ! or kernel.  Any kernel that reads a ghost that was not overwritten by
      ! the subsequent exchange will encounter a NaN, failing the run loudly at
      ! the offending step.  Default off (dyn%poison_ghosts = .false.) = the
      ! branch is never taken and the run is bit-identical to the unpoisoned path.
      ! Requires bc (edge topology); bc absent means no topology to gate on.
      if (dyn%poison_ghosts .and. present(bc)) then
         call ocean_poison_ghost_bands( &
            grid, ms, dyn%bt_work, ss, &
            (.not. bc%has_west) .or. bc%periodic_x, &
            (.not. bc%has_east) .or. bc%periodic_x, &
            (.not. bc%has_south) .or. bc%periodic_y, &
            (.not. bc%has_north) .or. bc%periodic_y)
      end if

      tide_on = .false.
      if (present(tides)) then
         if (tides%enable) tide_on = .true.
      end if
      ! Refresh the equilibrium-tide elevation ONCE per outer step, before
      ! the RK2 stages (tidal periods are O(hr) >> inner dt).
      if (tide_on) then
         if (present(t)) then
            call tides_update_eta_eq(tides, t)
         else
            call tides_update_eta_eq(tides, 0.0_wp)
         end if
         ! Scalar SAL (C2): fold `beta_sal*eta` into the combined seam field.
         ! Uses the LAGGED barotropic surface elevation `bt_work%bt_eta`
         ! (the SSH the barotropic substep re-derives each stage; here it
         ! still holds the previous outer step's value — standard
         ! one-step-lagged scalar SAL).  `use_sal=.false.` ⇒
         ! `eta_forcing == eta_eq`, bit-identical to C1.
         call tides_update_eta_sal(tides, dyn%bt_work%bt_eta)
      end if

      ! Atmospheric surface-pressure loading / inverse barometer (PR-17).
      ! Refresh the combined seam field ONCE per outer step (held static
      ! across the inner substep loop, exactly like the tide): fold
      ! `eta_ib = -p_surf/(rho0 g_bt)` into `eta_seam = eta_ib [+ tide
      ! eta_forcing]`.  `g_bt` is the barotropic substep gravity (what the
      ! seam consumer runs on).  `sf%p_surf` is the assembled total (read
      ! only here; seeded from p_surf_atm at configure, and — once PR-18
      ! lands — overwritten with the ice load via the ice path's own inout
      ! access before this step).  The seam is differenced in the ghost
      ! band, so a halo exchange follows the fill (single-rank no-op;
      ! correct multi-rank).  `sf%p_surf` is mapped under the PR-12
      ! component set; `validate_config` requires it (and `present(sf)`)
      ! when psurf is enabled, so the read is always device-present.
      psurf_on = .false.
      if (present(psurf) .and. present(sf)) then
         if (psurf%enable) psurf_on = .true.
      end if
      if (psurf_on) then
         if (tide_on) then
            call p_surf_update_seam(psurf, sf%p_surf, dyn%bt_work%g_bt, &
                                    eta_tide=tides%eta_forcing)
         else
            call p_surf_update_seam(psurf, sf%p_surf, dyn%bt_work%g_bt)
         end if
         call ocean_halo_centre(psurf%eta_seam, device_resident=.true.)
      end if

      ! E3: the SAME assembled `sf%p_surf` also becomes the top-of-column
      ! pressure the EOS's IN-SITU builders measure down from (`ms%p_top`,
      ! Pa) when `&ocean_psurf_nml in_eos` is set — the ice-shelf-cavity
      ! seam, where 1e6-2e7 Pa of overburden makes the historical
      ! "in-situ pressure starts at 0 at the free surface" a systematic
      ! ~4-5 kg/m^3 density error under a nonlinear EOS.  It does NOT
      ! touch the POTENTIAL density `ms%rho_layer`, which stays at the
      ! horizontally uniform `eos%p_ref` by design.  Refreshed here, once
      ! per outer step, before the PGF of the step and held static across
      ! the stages.
      !
      ! P5.0: the SAME refresh also covers `&ocean_pgf_nml p_top_in_bc`,
      ! the second consumer of `ms%p_top` — it puts the load in the
      ! FV_MOM6 pressure-stack surface BC (`pa(nz+1)`), independently of
      ! whether it also reaches the EOS arguments.  The gate is the
      ! DISJUNCTION so neither consumer can ever read a p_top that the
      ! configure-time seed left behind while `sf%p_surf` moved on.
      !
      ! P5.2 completes the PARTITION: the assembled top-of-column load is
      !
      !     ms%p_top = metrics%p_ice_ref + sf%p_surf
      !
      ! — the STATIC isostatic ice load (absorbed into the barotropic
      ! datum `bt_H_ref = b - z_draft`, and therefore deliberately NOT a
      ! component of `sf%p_surf`, which the `eta_ib` seam is built from)
      ! plus whatever atmospheric / anomaly load the psurf seam carries.
      ! Without the cavity term this refresh would OVERWRITE the
      ! configure-time ice load with `p_surf` alone on step 1, silently
      ! unloading the column for every consumer of `p_top`.
      !
      ! A cavity WITHOUT the psurf seam needs no refresh at all: the
      ! draft is static, so the configure-time seed in
      ! `configure_ocean_cavity` is already the final value and this
      ! whole block stays switched off (`psurf_on = .false.`).  That is
      ! why the gate below is still the psurf gate.
      !
      ! Written INLINE as a `do concurrent` rather than as a call: a
      ! host-gated call handing a state array to an EXTERNAL subroutine
      ! makes nvfortran treat that array as escaping and pessimises EVERY
      ! `do concurrent` in this routine, whether or not the branch is
      ! taken (CLAUDE.md, measured at +4.8 % on `ocean_continuity`).
      ! The copy spans the WHOLE array, ghosts included, so `p_top`
      ! inherits exactly the halo validity `p_surf` has and needs no
      ! exchange of its own (`p_ice_ref` is ghost-filled at configure,
      ! from a `z_draft` that went through the bathymetry's own re-wrap
      ! and halo exchange).
      p_top_live = .false.
      if (psurf_on) then
         if (psurf%in_eos .or. pgf%p_top_in_bc) p_top_live = .true.
      end if
      if (p_top_live) then
         nx_ptop = size(ms%p_top, 1)
         ny_ptop = size(ms%p_top, 2)
         ! Two inline loops rather than one with a branch on `use_cavity`:
         ! without a cavity `p_ice_ref` is a `(1,1)` PLACEHOLDER, so the
         ! cavity spelling may not even be written in a form the compiler
         ! could speculate an index out of.
         if (metrics%use_cavity) then
            do concurrent(j=1:ny_ptop, i=1:nx_ptop)
               ms%p_top(i, j) = metrics%p_ice_ref(i, j) + sf%p_surf(i, j)
            end do
         else
            do concurrent(j=1:ny_ptop, i=1:nx_ptop)
               ms%p_top(i, j) = sf%p_surf(i, j)
            end do
         end if
      end if

      call probe_dS(grid, ms, "outer step entry", 0, dyn%outer_step_count + 1)

      ! Fox-Kemper MLE restratification (B5): compute the ML-confined
      ! overturning transports ONCE per outer step at THERMO cadence (it
      ! is a slow buoyancy-driven process; per-stage recompute doubles
      ! cost for no benefit).  Uses the MLD diagnosed by EPBL at the end
      ! of the previous step.  `run_stage_split` then folds the same
      ! uhml/vhml into the mass fluxes in both stages.  No-op when
      ! absent / disabled.
      if (present(mle) .and. present(epbl)) then
         if (dyn%enable_thermodynamics .and. dyn%is_thermo_step()) then
            call profiler_start("ocean_foxkemper")
            ! Pass the thermo window so the per-layer FK availability cap
            ! keeps the windowed tracer-drain hprev reconstruction positive
            ! in thin (EPBL-MLD) surface layers (dt_tracer_advect_ratio > 1).
            ! At dt_therm_ratio = 1 this is dt — a tighter-but-harmless cap.
            if (present(bc)) then
               call mle_compute_transports(grid, metrics, mle, ms, epbl, ss=ss, &
                                           dt_limit=dyn%therm_dt(dt), bc=bc)
            else
               call mle_compute_transports(grid, metrics, mle, ms, epbl, ss=ss, &
                                           dt_limit=dyn%therm_dt(dt))
            end if
            call profiler_stop("ocean_foxkemper")
         end if
      end if

      ! Wave speed (B1): refresh the first-baroclinic gravity-wave speed
      ! `cg1` + the Rossby deformation radius `rd`/`rd_over_dx` ONCE per
      ! outer step at THERMO cadence (further gated by `n_wavespeed`,
      ! `mod(outer_step_count, n_wavespeed) == 0`; effective cadence is
      ! lcm(dt_therm_ratio, n_wavespeed) when both are > 1).  Placed at the
      ! TOP of the outer step — before both `varmix_compute` call sites
      ! and `run_meke_step` — so `cg1` is fresh for all three consumers
      ! this step; `rho_layer` is whatever the PREVIOUS outer step (or,
      ! at `outer_step_count == 0`, the initial condition) left, the same
      ! one-step lag GM/MEKE already accept and the same placement MOM6's
      ! `calc_resoln_function` uses relative to `step_MOM`.  At step 0
      ! this means `cg1 = 0` (if `rho_layer` is uniform-rho0 IC) so
      ! `Res_fn = 1` for one step — a documented decision, not a bug.
      ! No-op when absent / disabled ⇒ bit-identical.  `wavespeed_compute`
      ! is `pure` (called from `do concurrent`-adjacent code), so the
      ! profiler wrapper lives here, not inside it.
      if (present(wavespeed)) then
         if (wavespeed%is_init .and. wavespeed%enable) then
            if (dyn%enable_thermodynamics .and. dyn%is_thermo_step() .and. &
                mod(dyn%outer_step_count, wavespeed%n_wavespeed) == 0) then
               call profiler_start("ocean_wavespeed")
               call wavespeed_compute(grid, metrics, wavespeed, ms)
               call profiler_stop("ocean_wavespeed")
            end if
         end if
      end if

      ! Resolution-scaled momentum viscosity (Gap 1): the lateral-mix
      ! kernels read the VarMix resolution function `res_fn_u/v`.  When GM
      ! is active it owns the VarMix refresh (below), so this standalone
      ! pass runs ONLY when GM is absent / disabled but `resoln_scaled_visc`
      ! still needs a fresh `res_fn` — fired once per outer step at THERMO
      ! cadence (`res_fn` is slowly varying; it stays device-resident
      ! between thermo steps).  No-op when VarMix is disabled or the knob
      ! is off ⇒ bit-identical.  `res_fn` depends only on the static grid
      ! terms + cg1 (not on the slopes GM refreshes), so a no-GM run still
      ! yields the correct resolution function.
      if (present(varmix) .and. present(wavespeed) .and. &
          present(lateral_mix) .and. present(slopes)) then
         if (varmix%enable .and. lateral_mix%resoln_scaled_visc) then
            if (dyn%enable_thermodynamics .and. dyn%is_thermo_step()) then
               if (.not. gm_refreshes_varmix(gm, slopes, dyn)) then
                  ! `slopes` is guaranteed present here; `varmix_compute`
                  ! self-gates on `slopes%is_init` (the configure invariant
                  ! ensures VarMix-enabled ⇒ slopes allocated) and uses the
                  ! slope only for the (here-zero) Eady SN term — the
                  ! resolution function it fills needs only cg1.
                  call varmix_compute(grid, metrics, varmix, slopes, &
                                      wavespeed, ms)
               end if
            end if
         end if
      end if

      ! Gent-McWilliams thickness diffusion (capability [2]) — INPUTS: refresh
      ! the isopycnal slopes, the VarMix base KhTh and MEKE once per outer
      ! step at THERMO cadence (a slow buoyancy-driven process).  The bolus
      ! transport itself is NOT computed here: it is computed from the
      ! thickness the dynamics LEAVES and applied as its own operator after
      ! the stage loop (`run_gm_step`, below) — MOM6's `thickness_diffuse`
      ! after `step_MOM_dyn_split_RK2`.  The split driver does not otherwise
      ! call `ocean_slopes_compute`, so GM owns the slope refresh.  No-op
      ! when absent / disabled.
      if (present(gm)) then
         if (dyn%enable_thermodynamics .and. dyn%is_thermo_step()) then
            if (gm%enable .and. present(slopes)) then
               call profiler_start("ocean_gm")
               call ocean_slopes_compute(grid, metrics, eos, slopes, ms, dyn%therm_dt(dt))
               ! VarMix (capability [4]): fill the spatially-varying pre-CFL
               ! base KhTh face field from the fresh slopes + cg1; the GM
               ! operator consumes it as its external base.  When VarMix is
               ! absent / disabled GM falls back to its scalar khth.
               if (present(varmix) .and. present(wavespeed)) then
                  if (varmix%enable) then
                     call varmix_compute(grid, metrics, varmix, slopes, &
                                         wavespeed, ms)
                  end if
               end if
               ! MEKE (capability [5]): step the prognostic eddy-energy field
               ! AFTER VarMix so it reads the previous step's gm_src
               ! (one-step lag) and adds the geom-mean kh into varmix%khth
               ! before the GM operator's CFL clamp.  VarMix off ⇒ MEKE
               ! still evolves E (feedback inert).  No-op when absent /
               ! disabled.
               call run_meke_step(grid, metrics, gm, varmix, wavespeed, hv, &
                                  ms, dyn%therm_dt(dt), meke)
               call profiler_stop("ocean_gm")
            end if
         end if
      end if

      ! Redi neutral diffusion (capability [3]): build the tracer-
      ! INDEPENDENT neutral-surface coefficients ONCE per outer step at
      ! THERMO cadence (a slow geometry — recomputing per stage/tracer is
      ! wasteful).  `run_stage_split` then applies the rotated flux to each
      ! tracer (Phase B) after `tracer_hdiff`.  No-op when absent / disabled.
      if (present(redi)) then
         if (dyn%enable_thermodynamics .and. dyn%is_thermo_step()) then
            if (redi%enable) then
               call profiler_start("ocean_redi")
               call redi_calc_coeffs(grid, metrics, eos, redi, ms)
               call profiler_stop("ocean_redi")
            end if
         end if
      end if

      ! SPEC S1 seed (pred_corr): before the very first step the time-mean
      ! family must hold the initial state — the first predictor's
      ! CorAd/hvisc read u_av/h_av, and the allocation default (0) would
      ! hand them a zero-thickness field.
      !
      ! The seed MUST obey the land contract (`mask_time_mean_velocities`;
      ! the prognostic itself is left to the stage-end
      ! `mask_layer_velocities`, as before).  Nothing ever rewrites `u_av`
      ! at a masked face afterwards:
      ! its only writer is the transport renormaliser's `u_cor`, which
      ! skips the physical walls and has `Σ h·dy_cu = 0` (no write) on
      ! every land face.  A non-zero initial velocity there therefore
      ! lived in `u_av` for the whole run — and the three Coriolis forms
      ! read it differently: the enstrophy form's `v_at_u` sees it, the
      ! energy / HK transport forms do not (their `vh` rides the masked
      ! `dx_cv = 0`), and the fast-loop reference (`set_cor_ref_velocity`
      ! → `subtract_fast_cor_ref`) does.  Under `sadourny_energy` /
      ! `sadourny_hk` the reference then removed a Coriolis term the slow
      ! forcing never contained, a constant `−(f/4)·(v̄_av(i−1)+v̄_av(i))`
      ! on every barotropic substep of the wall-adjacent rows — measured
      ! ×126 in KE+PE over 4.2 d on `test_ocean_cor_ref_seiche`'s basin,
      ! its time-integrated work matching the gain to 0.1 %.  A masked IC
      ! (every configured run: `ocean_state_seed_land_cells`) makes the
      ! mask idempotent, i.e. byte-identical.
      if (dyn%split_scheme == SPLIT_SCHEME_PRED_CORR .and. dyn%outer_step_count == 0) then
         call copy_field_3d(ms%u_face_x_layer, ms%u_av_layer, &
                            size(ms%u_face_x_layer, 1), size(ms%u_face_x_layer, 2), &
                            size(ms%u_face_x_layer, 3))
         call copy_field_3d(ms%v_face_y_layer, ms%v_av_layer, &
                            size(ms%v_face_y_layer, 1), size(ms%v_face_y_layer, 2), &
                            size(ms%v_face_y_layer, 3))
         call mask_time_mean_velocities(metrics, ms)
         call copy_field_3d(ms%h_layer, ms%h_av_layer, &
                            size(ms%h_layer, 1), size(ms%h_layer, 2), &
                            size(ms%h_layer, 3))
      end if

      call save_state(ms)
      if (allocated(ms%tracers)) then
         do it = 1, size(ms%tracers)
            call copy_field_3d(ms%tracers(it)%hTr, ms%tracers(it)%hTr0, &
                               size(ms%tracers(it)%hTr, 1), &
                               size(ms%tracers(it)%hTr, 2), &
                               size(ms%tracers(it)%hTr, 3))
         end do
      end if

      ! Forward all optionals straight through to `run_stage_split` —
      ! Fortran 2008+ propagates `present(...)` correctly across
      ! optional dummy arguments, so an absent optional passed as the
      ! actual arg stays absent in the callee.  This collapses the
      ! previous 4-branch cartesian-product dispatch (vcoord × bc)
      ! into a single stage loop.
      !
      ! MOM6 `set_viscous_BBL` (`step_MOM_dynamics`, once per step before
      ! the predictor): the per-face bottom boundary layer — `kv_bbl`,
      ! `bbl_thick` — that the vdiff glue (and the visc_rem producer) read
      ! in both stages / both schemes.  From the start-of-step state.  No-op
      ! unless the per-face BBL glue is configured.
      call vdiff_set_viscous_bbl(grid, vd, ms, eos, cor%f_corner)
      do stage = 1, 2
         if (psurf_on .and. pgf%p_top_in_bc) then
            ! The load reaches the slow PGF too (`ms%p_top` in the FV_MOM6
            ! top BC), so its depth mean is already in `F_bt`: hand the
            ! stage `eta_ib` so the MOM6-split forcing counts it once
            ! (`set_fast_forcing_eta_pf`).
            call run_stage_split(grid, metrics, dyn, eos, cor, ct, pgf, hv, bd, ss, &
                                 va, hd, vd, vmix, ms, dt, n_inner, &
                                 sf=sf, geo=geo, stage=stage, vcoord=vcoord, bc=bc, sp=sp, t=t, &
                                 lateral_mix=lateral_mix, epbl=epbl, kshear=kshear, mle=mle, &
                                 redi=redi, varmix=varmix, vmix_tidal=vmix_tidal, meke=meke, &
                                 eta_forcing=psurf%eta_seam, td=td, cav=cav, &
                                 eta_pf_seam=psurf%eta_ib)
         else if (psurf_on) then
            ! `eta_seam` already includes the tide when it is on (folded in
            ! p_surf_update_seam above), so this single branch subsumes the
            ! tide-on case; the two branches below are the pre-PR-17 code.
            call run_stage_split(grid, metrics, dyn, eos, cor, ct, pgf, hv, bd, ss, &
                                 va, hd, vd, vmix, ms, dt, n_inner, &
                                 sf=sf, geo=geo, stage=stage, vcoord=vcoord, bc=bc, sp=sp, t=t, &
                                 lateral_mix=lateral_mix, epbl=epbl, kshear=kshear, mle=mle, &
                                 redi=redi, varmix=varmix, vmix_tidal=vmix_tidal, meke=meke, &
                                 eta_forcing=psurf%eta_seam, td=td, cav=cav)
         else if (tide_on) then
            call run_stage_split(grid, metrics, dyn, eos, cor, ct, pgf, hv, bd, ss, &
                                 va, hd, vd, vmix, ms, dt, n_inner, &
                                 sf=sf, geo=geo, stage=stage, vcoord=vcoord, bc=bc, sp=sp, t=t, &
                                 lateral_mix=lateral_mix, epbl=epbl, kshear=kshear, mle=mle, &
                                 redi=redi, varmix=varmix, vmix_tidal=vmix_tidal, meke=meke, &
                                 eta_forcing=tides%eta_forcing, td=td, cav=cav)
         else
            call run_stage_split(grid, metrics, dyn, eos, cor, ct, pgf, hv, bd, ss, &
                                 va, hd, vd, vmix, ms, dt, n_inner, &
                                 sf=sf, geo=geo, stage=stage, vcoord=vcoord, bc=bc, sp=sp, t=t, &
                                 lateral_mix=lateral_mix, epbl=epbl, kshear=kshear, mle=mle, &
                                 redi=redi, varmix=varmix, vmix_tidal=vmix_tidal, meke=meke, td=td, &
                                 cav=cav)
         end if
         ! pred_corr between-stage reset (SPEC §2): the predictor's
         ! provisional up/vp/hp are discarded — only u_av/v_av/h_av carry
         ! forward — and the corrector advances from u^n / h^n.  The
         ! stage-2 entry then re-derives the barotropic state from the
         ! restored layers (derive_bt_from_layers), so the corrector's
         ! btstep also starts from the step-entry state, as MOM6's does.
         if (dyn%split_scheme == SPLIT_SCHEME_PRED_CORR .and. stage == 1) then
            call restore_state(ms)
         end if
      end do

      if (dyn%split_scheme == SPLIT_SCHEME_PRED_CORR) then
         ! No SSP average: the corrector's single full-dt update IS the
         ! step (SPEC §1 fact 1, §5 trap 4 — averaging FB stages
         ! annihilates the internal wave, |R| = 0.089 at θ = 1.35).
      else
         call rk2_average(ms)
         if (dyn%check_h_positive) then
            call check_h_positive_or_die(grid, ms, "after rk2_average", 3, &
                                         dyn%outer_step_count + 1, check_layers=.true.)
         end if
         ! Land-face velocity reset on the averaged state (C4 / R4b.2).
         call mask_layer_velocities(grid, metrics, ms, bt_work=dyn%bt_work)
         ! Phase-2 vanished-layer velocity reset on the RK2-averaged state.
         ! Only for VCOORD_LAGRANGIAN; no-op otherwise (bit-identical).
         if (dyn%reset_vanished_u .and. present(vcoord)) then
            if (vcoord%coord_type == VCOORD_LAGRANGIAN) then
               call reset_vanished_layer_velocities(ms, &
                                                    isopycnal_vanish_tol(dyn%angstrom_h, pd_floor=ct%positive_definite))
            end if
         end if
         if (allocated(ms%tracers)) then
            do it = 1, size(ms%tracers)
               call rk2_average_field_3d(ms%tracers(it)%hTr0, ms%tracers(it)%hTr, &
                                         size(ms%tracers(it)%hTr, 1), &
                                         size(ms%tracers(it)%hTr, 2), &
                                         size(ms%tracers(it)%hTr, 3))
            end do
         end if
      end if
      call probe_dS(grid, ms, "after rk2_average", 3, dyn%outer_step_count + 1)

      ! ---- Gent-McWilliams thickness diffusion: its OWN sequential operator
      ! on the CURRENT thickness, after the dynamics (MOM6 `thickness_diffuse`),
      ! every outer step.  See `run_gm_step`.
      call run_gm_step(grid, metrics, dyn, ct, va, ms, dt, gm, slopes, varmix, &
                       wavespeed, vcoord, bc)

      ! Phase 2 (6b): advance the windowed-advect clock once per OUTER step
      ! (not per RK2 stage) when accumulating.  Reset to 0 by the drain.
      if (dyn%dt_tracer_advect_ratio > 1) ct%t_dyn_rel_adv = ct%t_dyn_rel_adv + dt

      ! Velocity housekeeping (E7): advective-CFL truncation then the
      ! absolute maxvel cap.  Applied after the RK2 average but before
      ! the ALE remap — so the remap sees a bounded velocity field.
      ! Both no-op when their knob <= 0.
      ! Phase-3: pass vanish_tol when cfl_ignore_vanished + VCOORD_LAGRANGIAN.
      if (dyn%cfl_ignore_vanished .and. present(vcoord)) then
         if (vcoord%coord_type == VCOORD_LAGRANGIAN) then
            call apply_velocity_truncation(ms, metrics, dt, dyn%cfl_trunc, dyn%maxvel, &
                                           dyn%ntrunc_step, &
                                           vanish_tol=isopycnal_vanish_tol(dyn%angstrom_h, &
                                                                           pd_floor=ct%positive_definite), &
                                           n_nanzero=dyn%n_nanzero_step, &
                                           clip_cell_metric=.true., &
                                           nan_i=nan_i, nan_j=nan_j, nan_k=nan_k, nan_is_u=nan_is_u, &
                                           nan_dx=nan_dx, nan_visc_cfl=nan_visc_cfl, nu_h=hv%nu_h)
         else
            call apply_velocity_truncation(ms, metrics, dt, dyn%cfl_trunc, dyn%maxvel, &
                                           dyn%ntrunc_step, n_nanzero=dyn%n_nanzero_step, &
                                           clip_cell_metric=.true., &
                                           nan_i=nan_i, nan_j=nan_j, nan_k=nan_k, nan_is_u=nan_is_u, &
                                           nan_dx=nan_dx, nan_visc_cfl=nan_visc_cfl, nu_h=hv%nu_h)
         end if
      else
         call apply_velocity_truncation(ms, metrics, dt, dyn%cfl_trunc, dyn%maxvel, &
                                        dyn%ntrunc_step, n_nanzero=dyn%n_nanzero_step, &
                                        clip_cell_metric=.true., &
                                        nan_i=nan_i, nan_j=nan_j, nan_k=nan_k, nan_is_u=nan_is_u, &
                                        nan_dx=nan_dx, nan_visc_cfl=nan_visc_cfl, nu_h=hv%nu_h)
      end if
      dyn%ntrunc_total = dyn%ntrunc_total + dyn%ntrunc_step
      ! Loud NaN-catch accounting: any non-zero is a producer 0/0 upstream.
      ! Actionable (FINDINGS.md, the global-tripolar-aquaplanet debugging
      ! session that had nothing but this bare count to go on): name the
      ! first non-finite face's location, its local cell size, and — since
      ! `hv%nu_h` was cheap to thread through — the local viscous CFL a
      ! plain scalar-nu_h Laplacian would see there, so an under-resolved
      ! viscous CFL (the actual root cause that day) is visible immediately
      ! instead of requiring a separate bisection.
      if (dyn%n_nanzero_step > 0) then
         dyn%n_nanzero_total = dyn%n_nanzero_total + int(dyn%n_nanzero_step, int64)
         write (output_unit, '("[nan-catch] outer step ", i0, ": zeroed ", i0, &
            &" non-finite face velocities; first at ", a, "-face (i=", i0, ", j=", i0, &
            &", k=", i0, "), local dx=", es10.3, " m, local nu_h*dt/dx^2=", es10.3, &
            &" (nu_h=", es10.3, " m2/s)")') &
            dyn%outer_step_count + 1, dyn%n_nanzero_step, merge("u", "v", nan_is_u), &
            nan_i, nan_j, nan_k, nan_dx, nan_visc_cfl, hv%nu_h
         flush (output_unit)
      end if
      ! Outer-step seam (stage 0): the RK2-averaged, truncation-bounded
      ! state the ALE remap + next step will consume.
      call chksum_state(grid, ms, dyn%chksum_probe, "post_trunc", 0, dyn%outer_step_count + 1)

      ! Phase 2 (6b) windowed horizontal tracer-advect DRAIN.  Only active
      ! when `dt_tracer_advect_ratio > 1`; the every-step path (ratio = 1)
      ! never accumulated, so this is a hard no-op there (bit-identical).
      ! Runs AFTER the RK2 average (hTr is the window-start frozen mass,
      ! h_layer is the RK2-averaged window-end thickness) and BEFORE the
      ! ALE remap (spec §(c) ordering: advect-then-reset → remap).
      !
      ! Fires when the window just filled (`is_tracer_advect_step`), OR as
      ! a mandatory FLUSH before any ALE remap if the window is non-empty
      ! (`t_dyn_rel_adv > 0`) — so no Lagrangian-then-remapped state is ever
      ! emitted with un-drained accumulated transport.  The integer-multiple
      ! config constraint makes the natural window-close coincide with the
      ! thermo step, so the flush path is the safety net, not the norm.
      ! TODO(flush): wire the same drain into the output/restart cadence in
      ! the driver for the end-of-run segment (here it is remap-aligned).
      if (dyn%dt_tracer_advect_ratio > 1) then
         if (dyn%is_tracer_advect_step() .or. &
             (present(vcoord) .and. dyn%is_thermo_step() .and. ct%t_dyn_rel_adv > 0.0_wp)) then
            call profiler_start("ocean_tracer_drain")
            if (present(bc)) then
               call continuity_tracer_drain(grid, metrics, ct, ms, dyn%dt_tracer_advect_ratio, bc=bc)
            else
               call continuity_tracer_drain(grid, metrics, ct, ms, dyn%dt_tracer_advect_ratio)
            end if
            call profiler_stop("ocean_tracer_drain")
            call probe_dS(grid, ms, "after tracer drain", 3, dyn%outer_step_count + 1)
         end if
      end if

      ! Lagrangian-then-remap: dynamics has advanced on the existing
      ! grid; now relayer onto vcoord%target_h.  Orchestrator no-ops
      ! for VCOORD_EULERIAN_Z, so the call is bit-safe when vcoord
      ! is in its default configuration.
      !
      ! Gated on the THERMO cadence (`is_thermo_step()`): with the default
      ! `dt_therm_ratio = 1` every outer step is a thermo step, so the
      ! remap fires every step exactly as before (bit-identical).  With
      ! `dt_therm_ratio > 1` the layers run Lagrangian (h_layer is
      ! prognostic; continuity advances it every step) and relayer onto
      ! `target_h` once per thermo interval — the MOM6 DT_THERM design
      ! (Adcroft & Hallberg 2006).  The Lagrangian state between remaps is
      ! self-consistent: diag output remaps to fixed-z independent of the
      ! layer grid, and restart checkpoints carry `outer_step_count` so a
      ! mid-Lagrangian resume re-aligns the cadence bit-exactly.
      if (present(vcoord) .and. dyn%is_thermo_step()) then
         call profiler_start("ocean_ale_remap")
         ! The grid time-filter relaxes once per thermo interval, so it
         ! must see the aggregate thermo dt (= dt at dt_therm_ratio=1, the
         ! bit-identity default).
         call ocean_apply_ale_remap_step(grid, vcoord, ms, dyn%bt_work%bt_eta, dyn%bt_work%bt_H_ref, &
                                         method=vcoord%remap_method, eos=eos, dt=dyn%therm_dt(dt))
         call check_remap_preconditions_or_die(grid, vcoord, ms%nz_ml, &
                                               dyn%outer_step_count + 1)
         call profiler_stop("ocean_ale_remap")
         call probe_dS(grid, ms, "after ALE remap", 3, dyn%outer_step_count + 1)
         ! z-level closed faces: the ALE remap is the LAST velocity writer
         ! of the outer step, and it is a COLUMN operator -- it redistributes
         ! momentum along a face column without consulting any horizontal
         ! mask.  Its own `min(h_L,h_R)` face column (see
         ! `remap_x_face_velocity`) already gives a closed layer an
         ! exactly-zero target so nothing is poured IN, but the layers
         ! ABOVE and BELOW it still shift, and a PPM reconstruction whose
         ! stencil straddles the gap can leave a non-zero value in the
         ! zero-thickness cell.  Re-assert the wall here: a closed face
         ! carries exactly zero normal velocity at the END of the step, not
         ! merely at the end of the last stage.  Gated inside
         ! `mask_layer_velocities`, so this is a no-op with the knob off.
         if (metrics%use_closed_faces) then
            call mask_layer_velocities(grid, metrics, ms, bt_work=dyn%bt_work)
         end if
      end if

      ! Ideal-age surface reset (PR-7): the Dirichlet BC `age = young_val`
      ! on k=nz must be the LAST operator to touch the top layer this outer
      ! step — applied per-RK2-stage it is halved by `rk2_average_field_3d`
      ! above (a Dirichlet condition inside an averaged sub-step is not a
      ! Dirichlet condition); applied before the ALE remap it is overwritten
      ! by the remap's vertical redistribution of subsurface age into the
      ! new top cell. So: after rk2_average, after continuity_tracer_drain,
      ! after the ALE remap, before the ghost-SSH refill (SSH-only, doesn't
      ! touch tracers). Thermo-cadence gated to compose with the remap;
      ! self-gates internally on `ms%idx_age <= 0`.
      if (dyn%is_thermo_step()) then
         t_now = 0.0_wp
         if (present(t)) t_now = t
         call ocean_ideal_age_reset_surface(grid, ms, &
                                            ocean_ideal_age_young_val(dyn%ideal_age_young_val, &
                                                                      dyn%ideal_age_sfc_growth_rate, t_now))
      end if

      ! Re-establish the zero-gradient free surface in the open-edge GHOST
      ! columns before the diagnostic manager reads the state.  The slow
      ! continuity drifts the ghost h_layer (array-edge fluxes) and the
      ! conservative remap preserves that drift, so `SSH = Σh_layer − b`
      ! blows up in the halo (worst at open corners) while the physical
      ! interior is healthy.  Physics-neutral: next step's stage-1
      ! `ocean_obc_fill_ghosts` re-fills these ghosts before any kernel
      ! reads them.  No-op for WALL/PERIODIC edges (bit-identical).
      if (present(bc)) call ocean_obc_refill_ghost_ssh(grid, bc, ms, dyn%bt_work%bt_H_ref)

      ! ---- Invariant I1′: `h <= H_VANISHED ⇒ hTr = h·c_live` ---------------
      ! THE enforcement point.  Every tracer writer above (surface flux, melt,
      ! sponge, hdiff, vdiff, vertical advection, the OBC ghost fills, the
      ! windowed drain) deposits content into whatever layer it was handed;
      ! only the ALE remap checked `h` on the way in.  Rather than ask forty
      ! kernel authors to remember the rule, establish it ONCE here, at the end
      ! of the outer step, over the whole tracer registry: every filler is
      ! pooled with its donor live layer and the pool mixed to one
      ! concentration, so the filler carries its donor's `c_live`.
      !
      ! Content is moved WITHIN the column (filler ↔ donor), so the column
      ! integral is unchanged and NOTHING is recorded in any budget — a
      ! contributor that always sums to zero is noise in the one instrument
      ! that detects real leaks.
      !
      ! A column with no sub-threshold layer is a textual no-op, so every
      ! sigma / z*-lite / eulerian_z configuration in the tree is bit-identical
      ! and pays only the sweep.
      call profiler_start("ocean_vanished_i1")
      call ms%enforce_vanished_content(grid%nx_total, grid%ny_total)
      call profiler_stop("ocean_vanished_i1")
      call check_vanished_invariant_or_die(grid, vcoord, ms, dyn%outer_step_count + 1)

      dyn%outer_step_count = dyn%outer_step_count + 1
   end subroutine ocean_dyn_step_split