multilayer_state_t Derived Type

   type :: multilayer_state_t
      !! Per-layer multilayer C-grid state.

      logical :: is_init = .false.
         !! True between `init` and `destroy`.  Prefer this to
         !! `allocated(...)` — tracks GPU device attachment too.

      ! Conservation-budget accumulators (host scalars; not device-mapped).
      real(wp) :: mass_out = 0.0_wp
         !! Cumulative mass (kg) that has left the domain through its open
         !! boundaries since t=0 (positive = outflow), accumulated per RK2
         !! stage from the continuity divergence (`flux_h_layer`) so the
         !! console mass `Error` closes to round-off even with open BCs.
      logical :: mass_out_efp_on = .false.
         !! Also accumulate `mass_out` as an order-invariant extended-fixed-
         !! point sum (`mass_out_efp`), set from `&ocean_diag_nml
         !! reproducing_sums`: the FP running sum's last digits depend on the
         !! decomposition (it is a telescoping sum of large cancelling terms),
         !! the EFP one does not, so the console `out` column is the same on
         !! every rank count.
      integer(int64) :: mass_out_efp(EFP_DIGITS) = 0_int64
         !! This rank's `mass_out` as EFP bins (`rdb_efp` layout).
      integer(int64) :: mass_out_efp_poison = 0_int64
         !! Non-finite counter for `mass_out_efp`, mirroring `efp_t%poison`
         !! (this accumulator is a raw bin array, not an `efp_t`, since it
         !! is a standalone module-level running total rather than a
         !! collective-combined value -- see `efp_t`'s docstring in
         !! `rdb_efp`). `ocean_accumulate_mass_out` adds a slab's poison
         !! count here the same way it adds the slab's carried bins;
         !! `ocean_console_stats_report` folds it into
         !! `efp_local(IX_MOUT)%poison` so a NaN/Inf flux_h_layer poisons
         !! the console's `Mass out` column instead of laundering into a
         !! plausible finite number.
      logical :: mass_out_tracked = .false.
         !! Set once the dyn step has accumulated `mass_out`, so the console
         !! only activates the mass budget on a path that feeds it.
      real(wp) :: mass_src = 0.0_wp
         !! Cumulative mass (kg) ADDED to the domain since t=0 by a
         !! tracked volume SOURCE (positive = added), the mass twin of
         !! `salt_budget_surface` / `heat_budget_surface`.  Accumulated
         !! with the same per-stage weight as `mass_out`, so the console
         !! residual `(M - M0) + mass_out - mass_src` stays at round-off
         !! while the total legitimately grows.
         !!
         !! Fed today by the ice-shelf real-freshwater path
         !! (`&ocean_cavity_melt_nml freshwater="mass"`) and by its
         !! `volume_compensation` sink (which enters NEGATIVE).  Zero on
         !! every other path ⇒ the printed budget is unchanged.
         !!
         !! Scaled by `RHO_WATER`, not by the configured `rho_0`, because
         !! that is the density the console's own `total_mass =
         !! sum(h*areaT)*RHO_WATER` uses: the accumulator has to measure
         !! the same mass the total does.  The VOLUME it came from was
         !! converted from a kg/m^2/s flux with `rho_0` (Boussinesq
         !! volume conservation) — see `ocean_cavity_mass_step`.
         !!
         !! Host scalar, not device-mapped, and NOT restart-registered —
         !! the same policy `mass_out` follows (the registry carries
         !! device-mapped 2-D/3-D fields; these two cumulative host
         !! scalars restart at zero together with the console's own
         !! `mass0` reference latch, so the residual is measured over the
         !! resumed window rather than across the gap).

      integer :: nz_ml = 0
         !! Number of multilayer levels (k=1 bed, k=nz_ml surface).

      ! ---- Cell-centred per-layer thickness ----
      real(wp), allocatable :: h_layer(:, :, :)
         !! Layer thickness at cell centres (m), shape (nx, ny, nz_ml).
      real(wp), allocatable :: h_layer0(:, :, :)
         !! RK2 save of h_layer at start of outer step.

      ! ---- Face-located layer velocities + momenta ----
      ! Same face-indexing convention as the barotropic C-grid state.
      real(wp), allocatable :: u_face_x_layer(:, :, :)
         !! x-velocity at the WEST face of cell (i,j,k), shape
         !! (nx+1, ny, nz_ml): face i sits between cells i-1 and i, so a
         !! cell's divergence reads faces (i, i+1) — the convention every
         !! consumer (continuity `flux(i+1)-flux(i)`, `metrics%wet_u`,
         !! the kappa-shear centre average) actually uses.  ("east face"
         !! here previously was a stale docstring.)
      real(wp), allocatable :: hu_face_x_layer(:, :, :)
         !! h*u at the same west-face stagger (m^2/s).
      real(wp), allocatable :: v_face_y_layer(:, :, :)
         !! y-velocity at the SOUTH face of cell (i,j,k), shape
         !! (nx, ny+1, nz_ml): face j sits between cells j-1 and j
         !! (same stagger rule as `u_face_x_layer`).
      real(wp), allocatable :: hv_face_y_layer(:, :, :)
         !! h*v at the same south-face stagger.

      ! ---- MOM6 split-RK2 time-mean fields (docs/MOM6_SPLIT_RK2_SPEC.md §1) ----
      ! MOM6 carries an INSTANTANEOUS prognostic velocity AND a step time-mean
      ! (`CS%u_av`, `CS%h_av`).  Every slow tendency (CorAdCalc, horizontal
      ! viscosity) is evaluated on the TIME-MEAN, never on the prognostic —
      ! that is where the predictor-corrector gets its second-order character
      ! without an SSP stage average.  `u_av` is produced by continuity's
      ! `u_cor` (the transport-matched velocity satisfying Sum_k u*h = uhbt)
      ! and MUST NOT be written back into the prognostic.
      real(wp), allocatable :: u_av_layer(:, :, :)
         !! Step time-mean x face velocity, shape (nx+1, ny, nz_ml). MOM6 `u_av`.
      real(wp), allocatable :: v_av_layer(:, :, :)
         !! Step time-mean y face velocity, shape (nx, ny+1, nz_ml). MOM6 `v_av`.
      real(wp), allocatable :: h_av_layer(:, :, :)
         !! Step time-mean layer thickness, shape (nx, ny, nz_ml). MOM6 `h_av`.

      ! ---- RK2 face-velocity saves ----
      real(wp), allocatable :: u_face_x_layer0(:, :, :)
      real(wp), allocatable :: v_face_y_layer0(:, :, :)

      ! ---- Per-face per-layer mass fluxes (continuity-PPM output) ----
      real(wp), allocatable :: mass_flux_x_layer(:, :, :)
      real(wp), allocatable :: mass_flux_y_layer(:, :, :)

      ! ---- Cell-centred per-layer flux divergence ----
      ! continuity-PPM writes here; the apply step reads it for the
      ! forward-Euler h_layer update.  Shape (nx, ny, nz_ml).
      real(wp), allocatable :: flux_h_layer(:, :, :)

      ! ---- Cell-centred density ----
      ! EOS writes here from T, S; the PGF kernel reads it.  Shape
      ! (nx, ny, nz_ml), matches h_layer.
      real(wp), allocatable :: rho_layer(:, :, :)

      ! ---- Top-of-column pressure for the IN-SITU EOS builders ----
      ! (nx, ny) Pa, cell-centred INCLUDING ghosts, `>= 0`.  The pressure
      ! standing on the free surface that an IN-SITU pressure argument is
      ! measured DOWN from: a hydrostatic builder that used to start its
      ! column at 0 Pa starts it at `p_top(i,j)` instead.  Zero on every
      ! shipped configuration — filled only when `&ocean_psurf_nml
      ! in_eos = .true.`, from the assembled total
      ! `ocean_surface_flux_t%p_surf` (atmospheric load today, ice-shelf
      ! / sea-ice mass load later), refreshed once per outer step.
      ! ALWAYS allocated and zero-filled so every kernel has ONE code
      ! path: no optional dummy, no host-gated call handing a state array
      ! to an external subroutine, and `p_top(i,j) + p` is bit-exact `p`
      ! under IEEE-754 when the knob is off.
      !
      ! Three things it is NOT:
      !   * NOT the POTENTIAL-density reference.  `rho_layer` is evaluated
      !     at the scalar `eos%p_ref` and must stay that way — it is
      !     differenced along layers, so a per-column reference pressure
      !     would fabricate an along-layer density gradient.
      !   * NOT the barotropic `eta_forcing` seam — that is a separate,
      !     gradient-only, metres-valued field on `ocean_p_surf_t`, and it
      !     already carries the depth-uniform load gradient.
      !   * NOT the PGF top boundary condition (`pa(nz+1)`, still
      !     `rho_ref*g*eta`; `p_edge(nz+1)`, still 0).
      real(wp), allocatable :: p_top(:, :)

      ! ---- Vertical (cross-layer) velocity at layer interfaces ----
      ! Diagnostic in z*; reads as residual w under ALE.  Shape
      ! (nx, ny, nz_ml+1) with k=1 the bed (0) and k=nz_ml+1 the surface.
      real(wp), allocatable :: w_interface(:, :, :)

      ! ---- First LIVE layer, counting down from the top ----
      ! `k_top(i,j)` is the index of the shallowest layer that carries
      ! mass — the largest `k` with `h_layer(i,j,k) > H_VANISHED` — with
      ! a fallback of `nz_ml` when the column has no live layer at all
      ! (land, or a fully collapsed column).  It exists because under a
      ! quasi-geopotential coordinate beneath an ice shelf
      ! (`vcoord_type = "z_fixed"`, `vcoord%z_top > 0`) the layers whose
      ! nominal geopotential range lies INSIDE the ice are inert fillers
      ! at `zstar_h_min`, so on an ice-covered column `k = nz_ml` is NOT
      ! the ice-adjacent layer.  Every top-side consumer that used to
      ! spell `nz` literally reads this instead.
      !
      ! **Fallback `nz_ml` is what makes the indirection free.**  On
      ! sigma, z*-lite, and every other family shipped today no wet
      ! column ever vanishes its top layer, so `k_top ≡ nz_ml`, every
      ! rewritten consumer reads the same memory with the same
      ! arithmetic, and the answer is bit-identical.  A land column also
      ! reads `nz_ml` (all its layers sit AT the marker under the
      ! land-state contract), matching what those consumers do today.
      !
      ! **Static, and deliberately so.**  It is filled ONCE at configure
      ! (`configure_ocean_k_top`) from the `z_fixed` target at `eta = 0`
      ! — the same kernel and the same input the closed-face mask uses,
      ! so there is exactly one definition of "live".  Under `z_fixed`
      ! `eta` is absorbed by the first LIVE layer (the partial cell) and
      ! a filler's target is `zstar_h_min` whatever `eta` does, so the
      ! live/filler pattern does not move and there is nothing to
      ! recompute per stage.  `tests/test_ocean_ktop.F90` pins the
      ! static claim against the mask builder's own live pattern.
      integer, allocatable :: k_top(:, :)
         !! Shallowest live layer at cell centres, shape (nx, ny).
      integer, allocatable :: k_top_u(:, :)
         !! u-face twin, shape (nx+1, ny).  `min` of the two bounding
         !! columns, NOT `max`: a velocity face carries water in layer
         !! `k` only where BOTH abutting columns are live there — that is
         !! the same statement `metrics%open_u` makes — so the shallowest
         !! layer the FACE has is the DEEPER of the two column tops, i.e.
         !! the smaller index.  Taking `max` would put the ice-ocean drag
         !! and the implicit stress fold on a row that is a filler on one
         !! side.
      integer, allocatable :: k_top_v(:, :)
         !! v-face twin, shape (nx, ny+1).  Same `min` rule.

      ! ---- First LIVE layer, counting UP from the bed ----
      ! `k_bot(i,j)` is the bed-side mirror of `k_top`: the index of the
      ! deepest layer that carries mass — the smallest `k` with
      ! `h_layer(i,j,k) > H_VANISHED` — with a fallback of `1` when the
      ! column has no live layer at all.  Under `vcoord_type = "z_fixed"`
      ! the layers whose nominal geopotential range lies INSIDE the bed
      ! are inert fillers at `zstar_h_min`, so on every column shallower
      ! than `z_fixed_h_ref` `k = 1` is NOT the bed-adjacent layer.
      ! Every bed-side consumer (bottom drag + its HBBL band and the
      ! implicit-fold rate, the vdiff bed row, geothermal, the tidal-
      ! mixing bed anchor, the MEKE bed speed, the bed-reaching shortwave
      ! residual) reads this instead of spelling `1`.
      !
      ! **Fallback `1` is what makes the indirection free.**  No other
      ! family has a STATIC bed filler: sigma / z*-lite / zsigma never
      ! vanish a layer, and the dynamic bed pinch of `zstar_full` is not
      ! a configure-time pattern (consumers that care keep their own
      ! `h`-gated scan above `k_bot`), so `k_bot ≡ 1` there and every
      ! rewritten consumer reads the same memory with the same
      ! arithmetic.  A land column also reads `1`.
      !
      ! **Static, and deliberately so** — exactly `k_top`'s argument:
      ! filled ONCE at configure (`configure_ocean_k_bot`) from the
      ! `z_fixed` target at `eta = 0`, the same kernel and input the
      ! closed-face mask and `k_top` use; the bed is static and `eta` is
      ! absorbed by the first live layer at the TOP, so the bed-side
      ! live/filler pattern does not move.  Derived from bathymetry +
      ! the vcoord config, so it is rebuilt on every start and is NOT
      ! restart state.  `tests/test_ocean_zfixed_k_bot.F90` pins it.
      integer, allocatable :: k_bot(:, :)
         !! Deepest live layer at cell centres, shape (nx, ny).
      integer, allocatable :: k_bot_u(:, :)
         !! u-face twin, shape (nx+1, ny).  `max` of the two bounding
         !! columns — the mirror of `k_top_u`'s `min`: a face carries
         !! water in layer `k` only where BOTH columns are live there
         !! (`metrics%open_u`), so the deepest layer the FACE has is the
         !! SHALLOWER of the two column bottoms, i.e. the larger index.
         !! `min` would put the bottom drag and the implicit-drag fold on
         !! a row that is a filler on one side (a closed face).
      integer, allocatable :: k_bot_v(:, :)
         !! v-face twin, shape (nx, ny+1).  Same `max` rule.

      ! ---- Land / ocean mask (surface-forcing mask) ----
      ! 2D wet-cell indicator at cell centres: 1.0 = ocean, 0.0 = land.
      ! Populated at IC time from the bathymetry threshold; consumed by
      ! the three surface-forcing kernels (`ocean_surface_stress`,
      ! `ocean_bottom_drag`, `ocean_surface_flux`) so forcing doesn't
      ! drive fictitious land-column currents.  Default 1.0 everywhere
      ! is bit-identical.  Not yet plumbed through the prognostic
      ! kernels (continuity, PGF, advection, viscosity, vdiff).
      real(wp), allocatable :: wet_mask(:, :)

      ! ---- Tracer registry ----
      ! Reuses the coastal `tracer_t` verbatim, cell-centred shape
      ! `(nx, ny, nz)` (stagger only affects velocities).  Salinity +
      ! temperature first; passive tracers append.  Special-shape
      ! kernels (EOS, surface flux) locate S/T via the named indices.
      type(tracer_t), allocatable :: tracers(:)
         !! Registered prognostic tracers.
      integer :: idx_salinity = 0
         !! Index into tracers(:) for salinity. 0 = not registered.
      integer :: idx_temperature = 0
         !! Index into tracers(:) for temperature. 0 = not registered.
      integer :: idx_age = 0
         !! Index into tracers(:) for the ideal-age tracer.  0 = not
         !! registered.  When > 0, the dyn step ages at 1 s/s and zeros
         !! the surface layer (k = nz_ml) every step.  No EOS coupling.
      integer :: idx_pseudo_salt = 0
         !! Index into tracers(:) for the pseudo-salt verification tracer.
         !! 0 = not registered.
      logical :: registry_locked = .false.
         !! Set by `enter_data`, cleared by `exit_data`.  While locked,
         !! `register_passive_tracer` REFUSES: the device map snapshots
         !! `tracers(:)` element-by-element, so a slot appended after the
         !! map has no device `hTr` and the first kernel touching it faults
         !! on the `mem:separate` build.

      ! ---- Budget contributor slots ----
      ! Per-kernel mass / heat / salt accumulators: each physics kernel
      ! writes its per-step `dt · delta` here; the budget module's
      ! `drain_contributors` spatially integrates + zeroes them at each
      ! eval cadence.
      real(wp), allocatable :: mass_budget_continuity(:, :, :)
         !! Mass change per cell per step attributed to continuity-PPM
         !! divergence (m·dt units; the budget integral over volume
         !! recovers m³).  Shape (nx, ny, nz_ml).  Zero in a closed
         !! basin (perfect telescope).
      real(wp), allocatable :: heat_budget_surface(:, :, :)
         !! hTr (K·m) change per cell per step attributed to the surface
         !! heat-flux kernel.  Populated only at `k_top(i,j)`, the first
         !! LIVE layer — which is `nz_ml` on every column with no
         !! top-side filler, i.e. everywhere but an ice-covered column
         !! under a quasi-geopotential coordinate; zero elsewhere.  Sign
         !! convention: positive = source into the ocean.
      real(wp), allocatable :: salt_budget_surface(:, :, :)
         !! hTr (PSU·m) change per cell per step attributed to the
         !! surface salt-flux kernel.  Same shape + indexing convention
         !! as `heat_budget_surface` (`k_top`, not a literal `nz_ml`).
      real(wp), allocatable :: heat_budget_geothermal(:, :, :)
         !! hTr (K·m) change per cell per step attributed to the
         !! geothermal bottom-heat-flux kernel.  Populated only at the
         !! lowest massive layer (k=1 in the common case); zero
         !! elsewhere.  Sign convention: positive = source into the
         !! ocean from below.
      real(wp), allocatable :: heat_budget_sponge(:, :, :)
         !! hTr (K·m) change per cell per layer per step attributed to
         !! the map-driven sponge's tracer relaxation
         !! (`rdb_ocean_sponge::relax_tracer_budget_impl`).  Zero unless
         !! `&ocean_sponge_nml enable=.true., relax_tracers=.true.`.
         !! Sign convention: positive = source into the ocean (relaxing
         !! toward a warmer reference).  Not drained by the sponge.
      real(wp), allocatable :: salt_budget_sponge(:, :, :)
         !! hTr (PSU·m) change per cell per layer per step attributed to
         !! the map-driven sponge's tracer relaxation.  Same shape +
         !! indexing convention as `heat_budget_sponge`.
      real(wp), allocatable :: heat_budget_vert_adv(:, :, :)
         !! hTr (K·m) change per cell per layer per step attributed to
         !! the first-order-upwind vertical-advection kernel for
         !! temperature.  Closed-BC kernel: column sum telescopes to
         !! zero, so the spatial integral is zero to FP.
      real(wp), allocatable :: salt_budget_vert_adv(:, :, :)
         !! hTr (PSU·m) change per cell per layer per step for
         !! salinity.  Same closed-BC telescope property.
      real(wp), allocatable :: heat_budget_vdiff(:, :, :)
         !! hTr (K·m) change per cell per layer per step attributed to
         !! the backward-Euler vertical-diffusion tridiag solve for
         !! temperature.  Closed BCs (no flux through bed or surface)
         !! ⇒ column sum telescopes to zero.
      real(wp), allocatable :: salt_budget_vdiff(:, :, :)
         !! Salinity analogue of `heat_budget_vdiff`.
      real(wp), allocatable :: heat_budget_hdiff(:, :, :)
         !! hTr (K·m) change per cell per layer per step attributed to
         !! the Laplacian horizontal-diffusion kernel for temperature.
         !! Closed walls (wall faces forced to zero) ⇒ spatial integral
         !! telescopes to zero.
      real(wp), allocatable :: salt_budget_hdiff(:, :, :)
         !! Salinity analogue of `heat_budget_hdiff`.
      real(wp), allocatable :: heat_budget_horiz_adv(:, :, :)
         !! hTr (°C·m) change per cell per layer, accumulated (summed)
         !! across BOTH RK2 stages from t=0 (`+=` each stage, never zeroed
         !! mid-step; no 0.5 weight here — the console applies it),
         !! attributed to the continuity-PPM HORIZONTAL tracer advection
         !! (zonal + meridional divergence of the tracer mass flux).
         !! Interior sum telescopes to the net advective flux across the
         !! open boundaries; zero in a closed basin.  Only the fused
         !! (`dt_tracer_advect_ratio=1`) path fills it — see the console
         !! reporter for the >1 fallback.
      real(wp), allocatable :: salt_budget_horiz_adv(:, :, :)
         !! Salinity analogue (PSU·m) of `heat_budget_horiz_adv`.
      real(wp), allocatable :: mass_budget_remap(:, :, :)
         !! Per-cell `h_layer_new − h_layer_old` from the ALE remap
         !! step.  Column-conservative ⇒ Σ_k = 0 per (i, j).  Non-zero
         !! residual flags a remap conservation leak.
      real(wp), allocatable :: heat_budget_remap(:, :, :)
         !! Per-cell `hTr_new − hTr_old` from the ALE remap step for
         !! temperature.  Same column-telescope property.
      real(wp), allocatable :: salt_budget_remap(:, :, :)
         !! Salinity analogue.

   contains
      procedure, non_overridable :: init => multilayer_state_init
      procedure, non_overridable :: destroy => multilayer_state_destroy
      procedure, non_overridable :: enter_data => multilayer_state_enter_data
      procedure, non_overridable :: exit_data => multilayer_state_exit_data
      procedure, non_overridable :: bytes => multilayer_state_bytes
      procedure, non_overridable :: register_passive_tracer => &
         multilayer_register_passive_tracer
      procedure, non_overridable :: enforce_vanished_content => &
         multilayer_enforce_vanished_content
      procedure, non_overridable :: enforce_vanished_content_host => &
         multilayer_enforce_vanished_content_host
      procedure, non_overridable :: scan_vanished_content => &
         multilayer_scan_vanished_content
   end type multilayer_state_t