barotropic_workstate_t Derived Type

   type :: barotropic_workstate_t
      logical :: is_init = .false.
         !! True between `init` and `destroy`.  Tracks GPU device
         !! attachment too — prefer to `allocated(...)` which only
         !! sees the host pointer.

      logical :: bt_substep_drag = .false.
         !! When `.true.`, driver fills `bt_rem_u/v` with the
         !! multiplicative BT-substep drag factor that damps `bt_ubt`
         !! / `bt_vbt` each inner step.  Default off ⇒ `bt_rem` ≡ 1,
         !! multiplication is a no-op (bit-identical).

      logical :: bt_correction_visc_rem = .false.
         !! RETIRED-by-D1 (2026-10, PR-3 follow-up): when `.true.`,
         !! `apply_bt_correction` weights the per-layer barotropic
         !! increment by `visc_rem_*(k)/⟨visc_rem⟩_h` instead of
         !! uniformly, biasing the Δu distribution toward layers LESS
         !! damped by vertical viscosity (`&ocean_bt_nml
         !! correction_visc_rem`, now fail-loud at configure outside
         !! direct test construction of `cfg`). MOM6's split-explicit
         !! barotropic solver (Hallberg 1997; Hallberg & Adcroft 2009)
         !! gives every layer the SAME `u_accel_bt` via `accel_layer_u`
         !! (plus only the depth-mean-zero `pbce` baroclinic
         !! term) — NO `visc_rem` weight — so this fold has no MOM6
         !! counterpart, and on a 1-degree Southern Ocean z* OPEN-step
         !! probe it is the mechanism that NaNs at step ~40 under
         !! `bbl_glue` (the weight ratio is unbounded when a column's
         !! glue damping is uneven across layers; MOM6 never risks this
         !! because it never weights the fold at all).  `visc_rem_chain`
         !! does NOT set this field — see `bt_visc_rem_producer` below for
         !! the (now decoupled) producer gate.
      logical :: bt_visc_rem_producer = .false.
         !! D1 follow-up: gates the visc_rem PRODUCER fused into
         !! `vmix_apply_in_stage`'s momentum vdiff solve (MOM6
         !! `vertvisc_remnant`, sharing
         !! `vertvisc_coef`'s SAME coupling coefficients `a_u` — includes
         !! `kv_bbl`/the BBL glue and the Rayleigh/bed piston whenever the
         !! glue or `implicit_drag` folds them into the matrix; the
         !! producer itself does NOT require `implicit_drag` — see
         !! `vdiff_apply_momentum`'s `do_remnant`).  Decoupled from
         !! `bt_correction_visc_rem`: the producer must run whenever ANY
         !! consumer needs `visc_rem_u/v` fresh — `bt_forcing_visc_rem`
         !! (wt_u), `bt_renorm_visc_rem` (continuity u_cor), or
         !! `bt_rem_from_visc_rem` (av_rem/bt_rem) — not only the
         !! (retired) weighted BT-correction fold.  Set to the OR of all
         !! four (including the legacy `bt_correction_visc_rem`, so a
         !! test that constructs `cfg` directly and sets it still gets a
         !! live producer) by `configure_ocean_bt`.
      logical :: bt_forcing_visc_rem = .false.
         !! MOM6 `wt_u` parity for the BT forcing assembly: weight the
         !! `F_bt_u/v` depth-mean (and the PGF-projection subtraction) by
         !! `h_face·visc_rem(k)` (`&ocean_bt_nml forcing_visc_rem`).
      logical :: bt_renorm_visc_rem = .false.
         !! MOM6 continuity-inversion parity (SPEC S2b): γ-weighted
         !! transport-matching renormaliser (`&ocean_bt_nml
         !! renorm_visc_rem`) — `visc_rem_u/v` forwarded into the slow
         !! continuity so `u_cor = u + du·γ_k` and the fluxes carry the
         !! same weights.
         !!
         !! **PR-1 call-point / dt mapping to MOM6's three
         !! `vertvisc_remnant` call sites** (all at the OUTER step's
         !! `dt` — `VISC_REM_TIMESTEP_
         !! BUG` defaults `.false.`, so none of them use `dt_pred`):
         !!   * The pre-predictor call (`dt`) maps to
         !!     `visc_rem_precompute`'s pre-substep refresh in
         !!     `run_stage_split`, which always runs at the stage's `dt`
         !!     (gated `is_pc .or. bt_forcing_visc_rem .or.
         !!     bt_renorm_visc_rem .or. bt_rem_from_visc_rem`, i.e.
         !!     unconditionally once per `pred_corr` stage).
         !!   * The post-predictor call (full `dt`, NOT `dt_pred`) maps
         !!     to `vmix_apply_in_stage`'s stage-end producer called with
         !!     the `dt_remnant=dt` argument at the PREDICTOR stage —
         !!     decoupled from the predictor's own velocity-apply `dt_vel
         !!     = pc_be·dt` so the remnant matrix is built at the full
         !!     step, matching MOM6's default (non-buggy) behaviour.  Both
         !!     this and the pre-predictor call build the SAME linear
         !!     matrix (visc_rem does not depend on velocity, only on
         !!     dt/h/kv/drag), so they agree exactly — mirroring MOM6,
         !!     where both calls reuse the SAME `vertvisc_coef` output and
         !!     so are identical by construction.
         !!   * The corrector call (`dt`) maps to the stage-end producer's
         !!     existing fused call at the CORRECTOR stage, where `dt_vel
         !!     ≡ dt` already (no predictor off-centring) — unchanged.
         !! `ssp_rk2` (no predictor/corrector split): one refresh per
         !! stage, before that stage's barotropic step — the pre-substep
         !! `visc_rem_precompute` call (gated on `bt_forcing_visc_rem
         !! .or. bt_renorm_visc_rem .or. bt_rem_from_visc_rem`) plus the
         !! stage-end producer when `bt_visc_rem_producer` is on; `dt_vel
         !! ≡ dt` on both ssp_rk2 stages, so the split `dt_remnant` path
         !! is never taken there.

      logical :: bt_rem_from_visc_rem = .false.
         !! PR-2 (bt-rem-from-av-rem, `&ocean_bt_nml
         !! bt_rem_from_visc_rem`): `compute_bt_rem_from_visc_rem` builds
         !! `bt_rem_u/v` from `av_rem_u/v` (`:= Σ_k frhat_k·visc_rem_k`,
         !! `frhat_k` the face layer fraction `derive_bt_from_layers`
         !! uses) instead of the linear-piston law — MOM6's barotropic
         !! solver.  Mutually exclusive with
         !! `bt_substep_drag` (double-counted bed drag) and `bt_halo > 0`
         !! (checked in `validate_config`).
      logical :: bt_strong_drag = .false.
         !! MOM6 `BT_STRONG_DRAG`: the rational-
         !! approximation `bt_rem` form.  Requires `bt_rem_from_visc_rem`.
      logical :: bt_rescale_strong_drag = .false.
         !! MOM6 `RESCALE_STRONG_DRAG`.  Requires
         !! `bt_strong_drag`.

      integer :: frhat_scheme = FRHAT_ARITHMETIC
         !! `&ocean_bt_nml frhat_scheme` (`rdb_constants::FRHAT_*`,
         !! `parse_frhat_scheme`; default `FRHAT_ARITHMETIC`, bit-identical
         !! to the pre-port tree).  `FRHAT_HYBRID` ports MOM6's `btcalc`
         !! HVEL_SCHEME=HYBRID face-thickness closure
         !! into every barotropic depth
         !! mean that reads a layer's face thickness: `derive_bt_from_
         !! layers`, `face_depth_mean_u/v`, `face_depth_mean_rem_u/v`,
         !! `apply_bt_correction`'s open/visc_rem folds, and (through
         !! `face_depth_mean_u/v`) `compute_bt_rem_from_visc_rem`'s
         !! `av_rem`. Read on the host before each `do concurrent` (a plain
         !! scalar dispatch flag, like `bt_strong_drag` — never dereferenced
         !! from `bt_work` INSIDE a kernel, so it needs no `enter_data` map
         !! of its own). See `rdb_barotropic_coupling::frhat_h_face_step`
         !! (`src/shared_module_utilities/rdb_frhat_face.inc`) for the
         !! bottom-up port of MOM6's `e_u`/`D_shallow_u` sweep, and
         !! `derive_bt_from_layers`'s docstring for the one site
         !! (`use_upstream_h_face`) that has no MOM6 frhat counterpart and
         !! is deliberately left untouched.

      logical :: bt_correction_bc_pgf = .false.
         !! When `.true.`, `apply_bt_correction` adds the per-layer
         !! baroclinic-PGF retro-correction on top of the uniform /
         !! visc_rem-weighted Δu.  Requires `ocean_pgf_form = "fv_mom6"` and
         !! the `pbce` / `gtot_*` / `e_anom` / `eta_PF` fields filled.

      logical :: bt_bc_pgf_forcing = .true.
         !! `&ocean_bt_nml bc_pgf_forcing` (default `.true.`, MOM6
         !! parity).  When `.true.` the fast forcing is
         !! `F_bt_*_fast = F_bt + g_pf·∇(η_PF − η_seam)` — the depth
         !! mean of the FULL slow PGF stays in the forcing and only the
         !! free-surface term the slow PGF itself carries (`g_pf`,
         !! `pgf_free_surface_gravity`: 0 for MONT/FV_LITE/FV_WRIGHT) at
         !! the stage-entry η it was evaluated on is removed.  `.false.` =
         !! legacy `F_bt − ⟨PGF⟩`, which discarded the depth-mean
         !! baroclinic PGF.  See `set_fast_forcing_eta_pf`.

      real(wp) :: g_bt = 9.81_wp
         !! Acceleration in the barotropic-substep η-gradient PGF
         !! (`a = -g_bt · ∇η`), m/s².  Defaults to full gravity
         !! (bit-identical for Mont/FV/coastal).  For
         !! `OPGF_VARIANT_GPRIME` the driver overrides with
         !! `pgf%gprime_gfs` so the BT mode runs at the reduced-gravity
         !! speed `sqrt(g_FS · H)`; else the gprime knobs net out to
         !! full `g`.

      real(wp) :: bebt = 0.1_wp
         !! MOM6 `BEBT` (default 0.1, as MOM6).  Continuity flux uses
         !! `ubt_trans = (1+bebt)·ubt^n − bebt·ubt^{n-1}` — the
         !! `BT_PROJECT_VELOCITY` spelling.  Because η is updated BEFORE
         !! the velocity in this loop, it is the same scheme as MOM6's
         !! default (`BT_PROJECT_VELOCITY = .false.`: predictor η, then
         !! `(1−bebt)·ubt^n + bebt·ubt^{n+1}` transport): this loop's η is
         !! MOM6's `eta_pred` and the velocity sequence is identical, so
         !! the per-substep damping `|λ|² = 1 − bebt·a²` and the stability
         !! limit `a ≤ 2/√(1+2·bebt)` are MOM6's.  `bebt = 0` ⇒
         !! `ubt_trans = ubt^n` (neutral forward-backward Euler).  Set via
         !! `&ocean_bt_nml bebt`.

      logical :: substep_zeta_ke = .true.
         !! Live `(ζ_bt+f)·v − ∇KE` in the fast loop (default,
         !! bit-identical).  `.false.` = MOM6-parity planetary-only
         !! substeps; the `subtract_fast_cor_ref` reference reduces to
         !! `f·v̄` to match.  See `&ocean_bt_nml substep_zeta_ke`.

      ! ---- Primary barotropic-substep state ----
      real(wp), allocatable :: bt_eta(:, :)
         !! Barotropic SSH at cell centres, shape (nx, ny).  Defined
         !! as `sum_k(h_layer) - bt_H_ref`.
      real(wp), allocatable :: bt_ubt(:, :)
         !! Depth-mean u at east faces, shape (nx+1, ny).
      real(wp), allocatable :: bt_vbt(:, :)
         !! Depth-mean v at north faces, shape (nx, ny+1).
      real(wp), allocatable :: bt_H_ref(:, :)
         !! Reference column thickness (m) at cell centres, shape
         !! (nx, ny).  Constant for Eulerian-z; set by the driver
         !! from the initial multilayer state.  Used as `H` in the
         !! linearized continuity `∂η/∂t = -div(H · u_bt)`.

      ! ---- Fast-mode time-averaging accumulators ----
      real(wp), allocatable :: ubt_sum(:, :)
      real(wp), allocatable :: vbt_sum(:, :)
      real(wp), allocatable :: eta_sum(:, :)

      ! ---- Time-mean depth-integrated transport ----
      ! Accumulated each fast substep as `(H_ref + η_inst) · u_bt_inst`
      ! at faces, time-averaged at loop end; consumed by slow continuity
      ! as a barotropic constraint so per-layer mass fluxes vertically
      ! sum to the substep transport (makes `sum_k(h_layer) = H_ref +
      ! η_end`, so `apply_bt_correction`'s h-rescale a no-op).
      real(wp), allocatable :: uhbt_sum(:, :)
         !! East-face transport accumulator, shape `(nx+1, ny)`.
      real(wp), allocatable :: vhbt_sum(:, :)
         !! North-face transport accumulator, shape `(nx, ny+1)`.
      real(wp), allocatable :: bt_uhbt(:, :)
         !! Time-mean east-face transport (m²/s), shape `(nx+1, ny)`.
         !! `uhbt_sum / n_inner`.
      real(wp), allocatable :: bt_vhbt(:, :)
         !! Time-mean north-face transport (m²/s), shape `(nx, ny+1)`.

      ! ---- End-of-barotropic-substep snapshot ----
      ! Saved BEFORE the time-mean overwrite at fast-loop end.  Used by
      ! `apply_bt_correction` so the recombined per-layer momentum AND
      ! the h_layer rescale both see the end-of-step barotropic mode
      ! (Hallberg 2009 split-explicit convention); mixing end-step
      ! velocity with time-mean SSH corrupts gravity-wave dispersion.
      real(wp), allocatable :: bt_ubt_end(:, :)
      real(wp), allocatable :: bt_vbt_end(:, :)
      real(wp), allocatable :: bt_eta_end(:, :)

      ! ---- Nonlinear barotropic-substep scratch ----
      real(wp), allocatable :: bt_zeta_corner(:, :)
         !! Barotropic relative vorticity at corners, shape (nx+1, ny+1).
      real(wp), allocatable :: bt_ke_centre(:, :)
         !! Barotropic kinetic energy at cell centres, shape (nx, ny).
      real(wp), allocatable :: bt_eta_new(:, :)
         !! Per-substep η^{n+1} scratch, shape (nx, ny).  Jacobi
         !! buffer so the face-thickness reads in Pass 1 stay race-
         !! free under the Fortran 2018 `do concurrent` semantics.

      ! ---- Coriolis/advection reference velocity (MOM6 `ubt_Cor`) ----
      ! The barotropic velocity at which `subtract_fast_cor_ref`
      ! evaluates the reference `(ζ+f)·v − ∇KE` it removes from the
      ! substep forcing.  It MUST be the depth mean of the SAME layer
      ! velocity the slow `cor%pv_flux_*` in `F_bt` was evaluated on,
      ! or the uncancelled residual `f × (v̄_ref − v̄_slow)` is injected
      ! into every barotropic substep as a near-constant forcing and
      ! pumps the basin's gravest Poincaré seiche exponentially.
      !   * `ssp_rk2` — slow tendencies are evaluated on the prognostic
      !     `u^n`, so this is a copy of the stage-entry `bt_ubt/bt_vbt`
      !     (bit-identical to reading `bt_ubt/bt_vbt` directly).
      !   * `pred_corr` — slow tendencies are evaluated on `u_av/v_av`,
      !     so this is the depth mean of `u_av/v_av` under the SAME
      !     weights the forcing depth-mean used (h, or h·visc_rem when
      !     `&ocean_bt_nml forcing_visc_rem`).  MOM6 `MOM_barotropic`
      !     builds `Cor_ref_u` from `ubt_Cor = Σ_k wt_u·U_Cor` with
      !     `U_Cor = u_av`, i.e. the same velocity `CorAdCalc` used.
      real(wp), allocatable :: cor_ref_u(:, :)
         !! u-face Coriolis/advection reference velocity, shape (nx+1, ny).
      real(wp), allocatable :: cor_ref_v(:, :)
         !! v-face counterpart, shape (nx, ny+1).

      ! ---- BEBT projection: previous-substep velocity snapshots ----
      ! Hold u^{n-1} / v^{n-1} between substeps for the η-update
      ! extrapolation `(1+bebt)·u^n − bebt·u^{n-1}`.  Initialised to
      ! `bt_ubt / bt_vbt` at the top of each outer step so the first
      ! substep has zero extrapolation (bit-identical when `bebt = 0`).
      real(wp), allocatable :: bt_ubt_prev(:, :)
         !! u-face velocity from the previous substep, shape (nx+1, ny).
      real(wp), allocatable :: bt_vbt_prev(:, :)
         !! v-face counterpart, shape (nx, ny+1).

      ! ---- Split-driver slow-tendency accumulators ----
      ! Only allocated when `init` is called with the optional `nz_ml`.
      real(wp), allocatable :: F_slow_u(:, :, :)
         !! Per-face slow u-acceleration sum, shape (nx+1, ny, nz_ml).
      real(wp), allocatable :: F_slow_v(:, :, :)
         !! Per-face slow v-acceleration sum, shape (nx, ny+1, nz_ml).
      real(wp), allocatable :: F_bt_u(:, :)
         !! Depth-mean of `F_slow_u`, shape (nx+1, ny).  Used in the
         !! `apply_bt_correction` subtraction to clean the slow-apply's
         !! bt projection out of every layer before installing the
         !! barotropic-substep bt mode.
      real(wp), allocatable :: F_bt_v(:, :)
         !! Depth-mean of `F_slow_v`, shape (nx, ny+1).
      real(wp), allocatable :: F_bt_u_fast(:, :)
         !! `F_bt_u` minus the bt projection of the PGF, shape
         !! (nx+1, ny).  Passed as the substep's `force_u` so the
         !! substep's own `-G·∂η/∂x` is the only bt PGF on the bt mode
         !! (else the slow + internal PGF stack to `-2G·∂η/∂x`,
         !! doubling the effective gravity-wave speed).
      real(wp), allocatable :: F_bt_v_fast(:, :)
         !! v counterpart, shape (nx, ny+1).
      real(wp), allocatable :: ubt_at_n(:, :)
         !! Depth-mean u at the start of the outer step, shape
         !! (nx+1, ny).  Used in the recombine step:
         !!   u^{n+1}(k) = u^*(k) + (⟨u_bt⟩ - ubt_at_n - dt·F_bt_u)
      real(wp), allocatable :: vbt_at_n(:, :)
         !! v counterpart, shape (nx, ny+1).

      ! ---- bc-PGF per-layer correction ----
      ! Per-layer baroclinic-PGF retro-correction for the η change
      ! during the BT substep: the slow PGF used `eta_PF`, but the BT
      ! substep evolves η to `bt_eta_end`, so each layer would feel a
      ! PGF based on stale `eta_PF`.  Adds back the layer-dependent
      ! response, scaled by `pbce(k)` minus column-mean `gtot_face`.
      ! Fills only when `ocean_bt_correction_bc_pgf = .true.`; else
      ! holds init zeros and the corrector skips it (bit-identical).
      real(wp), allocatable :: pbce(:, :, :)
         !! Per-layer pressure-anomaly gravity coefficient (m/s²),
         !! shape (nx, ny, nz_ml).  Column-mean equals `gtot_face`.
         !! Built by `compute_pbce`; requires the FV_MOM6 PGF.
      real(wp), allocatable :: gtot_E(:, :)
         !! Depth-weighted column average of `pbce` evaluated at the
         !! east face of cell `(i, j)`, shape (nx, ny).  Built by
         !! `compute_gtot_faces`.
      real(wp), allocatable :: gtot_W(:, :)
         !! West-face counterpart, shape (nx, ny).
      real(wp), allocatable :: gtot_N(:, :)
         !! North-face counterpart, shape (nx, ny).
      real(wp), allocatable :: gtot_S(:, :)
         !! South-face counterpart, shape (nx, ny).
      real(wp), allocatable :: e_anom(:, :)
         !! SSH anomaly relative to `eta_PF`, shape (nx, ny).
         !! `e_anom = 0.5·(bt_eta_end + bt_eta) − eta_PF`.
      real(wp), allocatable :: eta_PF(:, :)
         !! Snapshot of `bt_eta` taken just before the slow PGF is
         !! computed (= η the PGF "saw").  Updated each RK2 stage.
         !! Shape (nx, ny).

      ! ---- BT corrector visc_rem weights ----
      ! Per-face per-layer fraction of velocity remaining after viscous
      ! damping over one outer step; biases the h-weighted corrector's
      ! per-layer Δu toward LESS-damped layers.  Default 1.0 ⇒ h-only
      ! path bit-identically.
      real(wp), allocatable :: visc_rem_u(:, :, :)
         !! u-face per-layer visc_rem, shape (nx+1, ny, nz_ml).
      real(wp), allocatable :: visc_rem_v(:, :, :)
         !! v-face per-layer visc_rem, shape (nx, ny+1, nz_ml).

      ! ---- bt-substep multiplicative drag damping ----
      ! Per-face damping factor in [0,1], applied inside the BT substep
      ! loop after each u/v update.  Computed once per stage from the
      ! linear-drag coefficient + face Htot:
      !   bt_rem_face = Htot / (Htot + r · HBBL · dt_inner)
      ! Cumulative over substeps ⇒ exp(-r·HBBL·t/Htot)-style damping.
      ! `bt_substep_drag` off ⇒ arrays hold 1, no-op (bit-identical).
      real(wp), allocatable :: bt_rem_u(:, :)
         !! u-face damping factor, shape (nx+1, ny).
      real(wp), allocatable :: bt_rem_v(:, :)
         !! v-face counterpart, shape (nx, ny+1).

      ! ---- PR-2 (bt-rem-from-av-rem): visc_rem depth mean ----
      ! `av_rem_u/v := Σ_k frhat_k·visc_rem_k`, the frhat-weighted (plain
      ! face-thickness-fraction, NOT visc_rem-weighted — see
      ! `compute_av_rem`'s docstring for the distinction from
      ! `face_depth_mean_rem_u`'s `wt_u`) depth mean of the viscous
      ! remnant, built once per barotropic call by
      ! `compute_bt_rem_from_visc_rem` when `bt_rem_from_visc_rem` is on
      ! (MOM6's barotropic solver).  Allocated unconditionally
      ! (cheap, 2D, same class as `bt_rem_u/v`); default 1.0 so an
      ! unused array is still well-defined if ever read.
      real(wp), allocatable :: av_rem_u(:, :)
         !! u-face visc_rem depth mean, shape (nx+1, ny).
      real(wp), allocatable :: av_rem_v(:, :)
         !! v-face counterpart, shape (nx, ny+1).

      ! ---- Barotropic linear (Rayleigh) wave drag (Egbert & Ray 2001;
      ! Jayne & St Laurent 2001) ----
      ! Static per-face piston velocity `r_H` [m/s] representing the
      ! barotropic-to-internal-tide energy sink.  MULTIPLIED into
      ! `bt_rem_u/v` by `compute_bt_rem_wave_drag` exactly as MOM6
      ! composes `lin_drag_u` with the viscous remnant.
      ! `lwd_enable = .false.` (default)
      ! ⇒ arrays stay unallocated and every BT path is bit-identical —
      ! same "arrays stay unallocated" contract as `use_bt_cont_type`
      ! below.  `lwd_` (not `wd_`) because `wd_` is taken by wet/dry.
      logical :: lwd_enable = .false.
         !! Driver writes from `&ocean_bt_nml wave_drag`.
      real(wp), allocatable :: lwd_drag_u(:, :)
         !! u-face piston velocity `r_H` [m/s], shape (nx+1, ny).
         !! Static after `configure_ocean_wave_drag`; allocated only
         !! when `lwd_enable = .true.`.
      real(wp), allocatable :: lwd_drag_v(:, :)
         !! v-face counterpart, shape (nx, ny+1).

      ! ---- BT_cont_type flux-bounded continuity ----
      ! Per-face piecewise-cubic flux closure replacing the naive
      ! `uh = u · h_face` (consumed by `find_uhbt`).  When
      ! `use_bt_cont_type = .false.` the arrays stay unallocated and
      ! every BT path is bit-identical.
      logical :: use_bt_cont_type = .false.
         !! Driver writes from the `ocean_use_bt_cont_type` namelist.
      type(local_BT_cont_u_type), allocatable :: BTCL_u(:, :)
         !! Per-u-face flux-closure coefficients, shape (nx+1, ny).
         !! Allocated only when `use_bt_cont_type = .true.`.
      type(local_BT_cont_v_type), allocatable :: BTCL_v(:, :)
         !! Per-v-face counterpart, shape (nx, ny+1).

      ! ---- Upstream-PPM face thickness for the BT chain ----
      ! Per-face column-sum of `h_layer` at the upstream face side,
      ! built once per outer step from `u/v_face_*_layer` signs and
      ! consumed across the BT chain (`derive_bt_from_layers`, substep,
      ! `apply_bt_correction`) so it shares the per-layer PPM
      ! face-thickness convention — kills phantom bed-layer velocity at
      ! slopes.  `use_upstream_h_face = .false.` ⇒ unallocated, legacy
      ! centred-h (bit-identical).
      logical :: use_upstream_h_face = .false.
         !! Driver writes from the `ocean_bt_upstream_h_face` namelist.
      real(wp), allocatable :: h_face_up_x(:, :)
         !! East-face upstream-h column sum (m), shape (nx+1, ny).
         !! `Σ_k h_layer(upstream_cell, j, k)` where the per-layer
         !! upstream selection follows `u_face_x_layer` sign.
         !! Allocated only when `use_upstream_h_face = .true.`.
      real(wp), allocatable :: h_face_up_y(:, :)
         !! North-face counterpart, shape (nx, ny+1).  Built from
         !! `v_face_y_layer` sign.

      ! ---- Dynamic wetting/drying (docs/ocean_wetdry_plan.md) ----
      ! Knobs + workspaces for the positive-definite wet/dry BT substep
      ! branch.  `wetdry_enable = .false.` (default) leaves every wd_*
      ! array unallocated and the substep on its unmodified centred-face
      ! path — byte-identical.  Allocated by `configure_ocean_wetdry`
      ! (rdb_ocean_setup) when `&ocean_wetdry_nml enable` is on.
      logical :: wetdry_enable = .false.
         !! Driver writes from `&ocean_wetdry_nml enable`.
      real(wp) :: wd_dry_depth = 0.05_wp
         !! Total-depth dry threshold (m); cells with `D < wd_dry_depth`
         !! are dynamically dry.  From `&ocean_wetdry_nml dry_depth`.
      real(wp) :: wd_rewet_depth = 0.10_wp
         !! Hysteresis re-wet threshold (m), > `wd_dry_depth`.  From
         !! `&ocean_wetdry_nml rewet_depth`.
      real(wp), allocatable :: wd_wet_dyn(:, :)
         !! Dynamic cell wet mask (1 = wet, 0 = dry), shape (nx, ny).
         !! PERSISTENT hysteresis state: updated every BT substep from
         !! `D = bt_H_ref + bt_eta` (wet above `wd_rewet_depth`, dry
         !! below `wd_dry_depth`, held in between).  Seeded from the
         !! initial D at configure.  Composes multiplicatively ON TOP of
         !! the static land masks (`metrics%wet_u/v` + zeroed metrics) —
         !! a static-land face can never be dynamically opened.
      real(wp), allocatable :: wd_theta(:, :)
         !! Per-cell positive-definite outflow limiter factor in [0, 1],
         !! shape (nx, ny).  `theta = min(1, available_volume /
         !! substep_outflow_volume)`; each face flux is scaled by
         !! `min(theta_L, theta_R)`, which guarantees `D >= 0` every
         !! substep (a cell drains at most what it holds).  Deep water
         !! ⇒ theta ≡ 1 ⇒ the limiter is exactly inert.
      real(wp), allocatable :: wd_flux_x(:, :)
         !! East-face provisional-then-limited volume flux (m³/s), shape
         !! (nx+1, ny).  Pass A fills `h_up · ubt_trans · dy_cu` with
         !! UPWIND face thickness + the FROUDE_CAP thin-face velocity
         !! guard; Pass C scales it by `min(theta_L, theta_R)` in place.
      real(wp), allocatable :: wd_flux_y(:, :)
         !! North-face counterpart, shape (nx, ny+1).
      real(wp), allocatable :: wd_open_u(:, :)
         !! Dynamic u-face open mask (1 = open, 0 = blocked), shape
         !! (nx+1, ny).  Bed-blocking (C-grid analogue of the coastal
         !! hydrostatic reconstruction): a face into a dry cell is a
         !! wall unless the wet side's surface stands above the dry
         !! side's bed elevation + dry_depth.  Rewetting starts from
         !! rest at the face.  Consumed by the substep Pass 2 gate and
         !! by the driver's layer-velocity masking.
      real(wp), allocatable :: wd_open_v(:, :)
         !! v-face counterpart, shape (nx, ny+1).
   contains
      procedure, non_overridable :: init => barotropic_workstate_init
      procedure, non_overridable :: destroy => barotropic_workstate_destroy
      procedure, non_overridable :: enter_data => barotropic_workstate_enter_data
      procedure, non_overridable :: exit_data => barotropic_workstate_exit_data
      procedure, non_overridable :: bytes => barotropic_workstate_bytes
   end type barotropic_workstate_t