continuity_t Derived Type

   type :: continuity_t
      logical :: is_init = .false.
         !! True between `init` and `destroy`.  Prefer this to
         !! `allocated(...)` — tracks GPU device attachment too.
      real(wp) :: angstrom_h = 0.0_wp
         !! Phase-1 Lagrangian minimum-thickness floor (m). Set from
         !! cfg%ocean%isopycnal%angstrom_h at setup, but only PASSED to the
         !! h-update kernels when the active vcoord is VCOORD_LAGRANGIAN
         !! (gated at the dyn call site). 0 ⇒ off ⇒ bit-identical.
         !! R7: floor lifts h but NOT hTr, so Tr=hTr/h shifts on a floored
         !! layer. Harmless for the adiabatic isopycnal config; thermo-on
         !! isopycnal correctness is OUT OF SCOPE for v1.
      logical :: conservative_floor = .false.
         !! When `.true.` (set from cfg%ocean%isopycnal%conservative_floor at
         !! setup; requires angstrom_h > 0 + VCOORD_LAGRANGIAN) the injecting
         !! `max(h, angstrom_h)` floor is skipped in the h-update and replaced
         !! by the conservative per-column borrow in `rdb_ocean_min_thickness`,
         !! invoked from the dyn continuity site. Uses `mt_h_new` as scratch.
         !! Default `.false.` ⇒ legacy injecting floor ⇒ bit-identical.
      real(wp) :: h_min = 1.0e-6_wp
         !! Lower clip on cell-centred thickness during the update.
         !! Also used as the floor for `ppm_limit_pos` when
         !! `use_ppm_limit_pos = .true.` — same semantic role as
         !! MOM6's `GV%Angstrom_H` for the vanishing-layer limiter.
      real(wp) :: cfl_max = 0.5_wp
         !! Soft cap on per-face CFL before falling back to upwind.
      logical :: renorm_legacy_single_step = .false.
         !! Use the pre-Newton single-linear-step uhbt renormalisation
         !! (donors picked at the UNCORRECTED velocity, no CFL bracket,
         !! no iteration) — wet/dry composition, set from
         !! `cfg%ocean%wetdry%enable` at setup.  The Newton form's donor
         !! RE-PICK + CFL bracket are what stabilise Lagrangian grounding
         !! (LAGRANGIAN_PGF_BUG.md §0), but under wet/dry they interact
         !! with the drying-front face gating and drive a drying column's
         !! `h_layer` negative (`test_ocean_wetdry_driver`); wet/dry owns
         !! its positivity via `ppm_limit_pos` + the BT limiter and ran
         !! validated on the single-step form, so it keeps it
         !! (bit-identical there).
      logical :: renorm_consistent_flux = .true.
         !! `&ocean_continuity_nml renorm_consistent_flux`.  When `.true.`
         !! the `uhbt`/`vhbt` renormalisation evaluates a layer whose upwind
         !! donor FLIPS under the correction as `(u0 + du)·h_face(new
         !! donor)`, i.e. as the flux of its corrected velocity, and
         !! brackets the Newton solve with bisection (MOM6
         !! `zonal_flux_adjust`).  The historical model
         !! `flux0 + du·h_face(new donor)` keeps the OLD donor's `u0·h_old`
         !! and is DISCONTINUOUS at the flip, by `u0·(h_new − h_old)·w`:
         !! whenever `uhbt` falls in that gap (a face where the corrected
         !! velocity must change sign across a thickness jump — a sigma
         !! layer over a bathymetric step, where `h_old ≠ h_new` after the
         !! PPM limiter flattens both edges), Newton has NO root, cycles
         !! for `RENORM_MAXIT` iterations and hands continuity a layer
         !! transport of the wrong SIGN.  The layer `η` then departs from
         !! the barotropic `η_end` by O(η) at the step every such step —
         !! a spurious η dipole that pumps the (undamped) barotropic
         !! grid-scale mode.  Default `.true.` (MOM6 behaviour); faces
         !! where no donor flips are bit-identical to the historical
         !! model, which `.false.` restores.
      logical :: use_ppm_limit_pos = .false.
         !! MOM6 `PPM_limit_pos` analogue.  When `.true.`, the PPM
         !! face-thickness reconstruction in continuity adds a
         !! positivity-preserving limiter that shrinks h_left /
         !! h_right toward h_centre whenever the parabolic fit would
         !! produce an interior minimum below `h_min`.  At the
         !! `h_centre ≤ h_min` limit the reconstruction collapses
         !! to a constant (= upwind for that cell), bounding the
         !! mass flux by the actual layer thickness.  Default off
         !! keeps the pre-knob behaviour bit-identical.  Driver writes
         !! from `cfg%ocean%continuity%ppm_limit_pos`.
      logical :: vol_cfl = .false.
         !! MOM6 `vol_CFL` analogue.  When `.false.` (default) the
         !! continuity flux uses the PPM downwind EDGE value (the
         !! CFL→0 limit), bit-identical to the pre-knob behaviour.
         !! When `.true.` the donor-side face thickness is the
         !! swept-volume integral of the reconstructed parabola
         !! (`volcfl_face`), adding the missing O(CFL) term.  Fixes
         !! the dt-sensitive near-bed residual at steep shelf breaks.
         !! Driver writes from `cfg%ocean%continuity%vol_cfl`.
      logical :: positive_definite = .false.
         !! Positive-definite split continuity master switch.  When `.true.`
         !! the split layer continuity (`continuity_tracer_step_split` via
         !! `continuity_zonal_flux`/`continuity_meridional_flux`) keeps every
         !! layer `>= h_lim`: P1 floors the PPM reconstruction edges at
         !! `2·h_lim` (MOM6's positive-definite reconstruction); P2 (later)
         !! scales down per-donor outfluxes so no mass is created.  Set from
         !! `cfg%ocean%continuity%positive_definite` at setup.  Host-side
         !! control scalar — read into a local before each DC loop (never
         !! dereferenced inside a device loop).  Default `.false.` ⇒ untaken
         !! branches only ⇒ bit-identical.
      real(wp) :: h_lim = 0.0_wp
         !! Positive-definite thickness floor (m).  Derived at setup:
         !! `cfg%ocean%isopycnal%angstrom_h` on VCOORD_LAGRANGIAN, else 0.
         !! Only consumed under `positive_definite = .true.`; at `h_lim = 0`
         !! (every non-Lagrangian coord) the P1 `max(edge, 2·h_lim)` floor is
         !! inert (edges are already `>= 0`), but the P2 outflux limiter still
         !! engages to keep every layer `>= 0` (avail = max(h, 0)).
      logical :: hTr_holds_conc = .false.
         !! True while the frozen tracer content `hTr` has been re-weighted
         !! onto the CURRENT `h_layer` so that the concentration
         !! `T = hTr/h_layer` is the (invariant) window-start value.  Set by
         !! `continuity_tracer_step_split` in TR_MODE_ACCUMULATE, cleared by
         !! `continuity_tracer_drain` once it has converted `hTr` back onto
         !! the reconstructed window-start thickness `hprev`.
         !!
         !! It exists so the drain's input contract stays EXPLICIT: `.false.`
         !! means "hTr is content on the window-start thickness" (the
         !! white-box unit tests that seed the drain directly, and the
         !! historical behaviour); `.true.` means "hTr is content on the
         !! current thickness".  Host-only — never read inside a kernel.
      logical :: windowed_advection = .false.
         !! Gates the ALLOCATION of the Phase-2 windowed-advection state (13
         !! 3D arrays: `uhtr`/`vhtr` accumulators + the (6b) drain
         !! workspace, ~4.2 GB at 1000x800x50) — consumed only when
         !! `dt_tracer_advect_ratio > 1`.  Latched from cfg BEFORE
         !! `init(grid)` by `ocean_state_init_from_config` (the same
         !! conditional-allocation contract as the default-off closures);
         !! the `enter_data`/`exit_data`/drain paths already guard on
         !! `allocated()`.
         !!
         !! Default `.false.`: the workspace is opt-in, so a `ct%init(...)`
         !! on any path that does NOT run through
         !! `ocean_state_init_from_config` costs nothing.  (It used to
         !! default `.true.` for the convenience of direct test call sites —
         !! which meant every such path silently paid the multi-GB
         !! allocation whether or not the drain could ever run.  Those call
         !! sites now set the flag explicitly before `init`.)
      integer :: n_limited_step = 0
         !! P2 diagnostic counter: number of interior faces whose mass flux
         !! the positive-definite limiter scaled (θ_donor < 1) this outer-step
         !! call, summed over the zonal + meridional passes.  Host-side scalar
         !! (the `!$acc parallel loop reduction` returns to the host); zeroed
         !! at the top of every `continuity_tracer_step_split` call.  P3 drains
         !! it to the console stats line.  0 ⇒ no limiting fired (the healthy
         !! case; MOM6-style graceful degradation must be loud).
      integer(int64) :: n_limited_total = 0_int64
         !! P3 running total of `n_limited_step` across every
         !! `continuity_tracer_step_split` call (accumulated once per call, at
         !! the end — one add per split call, mirroring `dyn%ntrunc_total`).
         !! The driver drains its per-report DELTA to the console next to the
         !! CFL-truncation line.  Grows ONLY when `positive_definite = .true.`
         !! (else the passes are skipped), so a nonzero total is itself the
         !! signal the limiter is active.  Loud-by-design caveat: a
         !! `uniform_z` isopycnal (VCOORD_LAGRANGIAN) stack shows a permanently
         !! LARGE, steadily-growing count because ~most layers sit AT the floor
         !! and θ=0 correctly freezes their (massless) outflow every call —
         !! that is not a pathology.  A healthy zstar/sigma run sits at 0.
         !! int64 (not int32 like `dyn%ntrunc_total`) precisely because that
         !! floored-stack case can accumulate > 2·10⁹ over a long run.
      integer :: tracer_recon = TRACER_RECON_PPM
         !! Face-reconstruction scheme for the WINDOWED tracer-advection
         !! drain (Q6; `dt_tracer_advect_ratio > 1` path only).  0 = CW-PPM
         !! (default, bit-identical), 1/2/3 = WENO5/7/9-Z swept-average
         !! (rung-adaptive, degrades near land/walls).  Set at configure
         !! from `&ocean_vmix_nml tracer_recon`.  Host-side control knob:
         !! the drain reads it to pick a face kernel — never dereferenced
         !! inside a device loop, so it rides the struct copyin.
         !! The every-step (ratio = 1) advect path is unaffected (stays
         !! CW-PPM) — WENO is a drain-only scheme in this wave.

      ! ---- Face-reconstruction workspace ----
      ! Phase 2b sizes these for the barotropic case (n3 = 1).  When
      ! Phase 5 adds layers, the multilayer kernel re-`init`s them
      ! at n3 = nz_ml.  Naming convention:
      !
      !   h_face_left_x(i, j, k)  — value AT east face i extrapolated
      !                             from the LEFT-side cell (i-1, j, k),
      !                             i.e. h_R of cell i-1.
      !   h_face_right_x(i, j, k) — value AT east face i extrapolated
      !                             from the RIGHT-side cell (i, j, k),
      !                             i.e. h_L of cell i.
      !
      ! Upwind: the kernel picks left if u >= 0 (left cell donates),
      ! right if u < 0.
      type(scratch_3d_buffer_t) :: h_face_left_x
         !! Left-state thickness at east faces.
      type(scratch_3d_buffer_t) :: h_face_right_x
         !! Right-state thickness at east faces.
      type(scratch_3d_buffer_t) :: h_face_left_y
         !! Left-state thickness at north faces.
      type(scratch_3d_buffer_t) :: h_face_right_y
         !! Right-state thickness at north faces.
      type(scratch_3d_buffer_t) :: mt_h_new
         !! Cell-centred `(nx,ny,nz)` scratch for the conservative
         !! minimum-thickness borrow (`conservative_floor`).  Holds the
         !! floor-only target thickness field between the h-update and the
         !! h_layer overwrite.  Unused (but allocated + mapped) when the knob
         !! is off — the DC kernels never touch it in that case.
      type(scratch_3d_buffer_t) :: mt_grounded
         !! `(nx,ny,1)` grounded-column mask (1.0 = any layer below floor)
         !! for the borrow's early-exit restructure: built in one coalesced
         !! pass, consumed by every borrow kernel in place of per-face
         !! column re-scans.  Same lifetime/mapping as `mt_h_new`.
      type(scratch_3d_buffer_t) :: pd_theta
         !! `(nx,ny,nz)` cell-centred per-donor availability factor θ(i,j,k)
         !! for the P2 positive-definite outflux limiter.  Built once per
         !! direction pass (over the FULL range incl. ghosts, so any cell that
         !! can donate to a swept face has a current θ), then each interior
         !! face is scaled by its upwind donor's θ.  Allocated + device-mapped
         !! UNCONDITIONALLY (like `mt_h_new`/`mt_grounded`) — the enter_data
         !! contract has no config visibility, and at `nx·ny·nz` reals the
         !! footprint matches `mt_h_new` already carried; unused (but present)
         !! when `positive_definite = .false.`.

      ! ---- Phase 2 flux-accumulator slots (DT_TRACER_ADVECT) ----
      ! Persistent device-resident accumulators of the per-stage C-grid
      ! face transports (`mass_flux_x/y_layer`).  When
      ! `dt_tracer_advect_ratio > 1` the windowed horizontal tracer
      ! advect (Phase 6b) drains these over the accumulation window;
      ! reconstruction `hprev = areaT·h_end + div(uhtr)` closes
      ! continuity by construction.  Allocated-but-unused at the default
      ! ratio = 1 (the every-step path bypasses them) — they are still
      ! mapped in enter_data so a ratio>1 run finds them present.
      ! Shapes mirror mass_flux_x/y_layer: east-face (nx+1,ny,nz),
      ! north-face (nx,ny+1,nz).
      real(wp), allocatable :: uhtr(:, :, :)
         !! Accumulated zonal face transport (m^3, area-weighted ·dt).
      real(wp), allocatable :: vhtr(:, :, :)
         !! Accumulated meridional face transport (m^3, area-weighted ·dt).
      real(wp) :: t_dyn_rel_adv = 0.0_wp
         !! Elapsed dynamics time (s) accumulated since the last tracer
         !! advect / accumulator reset.  Adds `dt` once per outer step.

      ! ---- Phase 2 (6b) windowed-drain workspace ----
      ! Persistent device-resident scratch for the fixed-budget PPM drain
      ! that spends `uhtr/vhtr` at the DT_TRACER_ADVECT boundary.  All
      ! allocated-but-unused at ratio = 1 (the every-step bypass never
      ! touches them); mapped in enter_data so a ratio > 1 run finds them
      ! present.  Sized at the same C-grid face / centre shapes as the
      ! prognostics.  Per-pass re-reconstruction (V2) writes Tr_work /
      ! pal / par / pa6 each sub-cycle pass (Reichl & Hallberg / MOM6
      ! ADVECT_PPM pattern); hprev_work carries the evolving thickness.
      real(wp), allocatable :: hprev_work(:, :, :)
      real(wp), allocatable :: h_win_start(:, :, :)
         !! Layer thickness captured when an accumulation window OPENS.
         !! Paired with `hTr_holds_conc`: the per-stage concentration hold
         !! re-weights `hTr` onto the evolving `h_layer`, and the drain undoes
         !! it with THIS array, which returns `hTr` bit-for-bit to the frozen
         !! window-start content the drain has always consumed.  (Undoing via
         !! the reconstructed `hprev` instead would be equivalent only where
         !! `hprev == h_win_start` exactly, and would silently convert any
         !! reconstruction residual — e.g. from a thin-layer `h_min` clip —
         !! into a tracer-mass drift.)
         !! Evolving (drained) layer thickness during the sub-cycle (m).
      real(wp), allocatable :: uhr_x(:, :, :)
         !! Remaining unspent zonal transport this window (m^3).
      real(wp), allocatable :: uhr_y(:, :, :)
         !! Remaining unspent meridional transport this window (m^3).
      real(wp), allocatable :: uhh_x(:, :, :)
         !! Per-pass limited zonal transport portion (m^3).
      real(wp), allocatable :: uhh_y(:, :, :)
         !! Per-pass limited meridional transport portion (m^3).
      real(wp), allocatable :: tr_flux_x(:, :, :)
         !! Per-pass zonal tracer flux F (m^3 · concentration).
      real(wp), allocatable :: tr_flux_y(:, :, :)
         !! Per-pass meridional tracer flux F (m^3 · concentration).
      real(wp), allocatable :: tr_work(:, :, :)
         !! Current concentration Tr = hTr/hprev_work (rebuilt per pass).
      real(wp), allocatable :: pal(:, :, :)
         !! CW parabola left-edge value per cell (rebuilt per pass).
      real(wp), allocatable :: par(:, :, :)
         !! CW parabola right-edge value per cell (rebuilt per pass).
      real(wp), allocatable :: pa6(:, :, :)
         !! CW parabola curvature a6 = 6·Tr − 3·(aL+aR) (rebuilt per pass).
   contains
      procedure, non_overridable :: init => continuity_init
      procedure, non_overridable :: destroy => continuity_destroy
      procedure, non_overridable :: enter_data => continuity_enter_data
      procedure, non_overridable :: exit_data => continuity_exit_data
      procedure, non_overridable :: bytes => continuity_bytes
   end type continuity_t