ocean_vcoord_t Derived Type

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

      ! ---- Active coordinate ----
      integer :: coord_type = VCOORD_EULERIAN_Z
         !! Selected vertical-coordinate variant.

      ! ---- Per-layer fractional thickness ----
      ! Sums to 1.0 across the column.  For `VCOORD_SIGMA` this is the
      ! target σ stencil (and `VCOORD_Z_FIXED` / `VCOORD_ZSTAR` fall back
      ! to it when no nominal profile was resolved, `z_fixed_h_ref = 0`):
      !   target_h(i,j,k) = (H(i,j) + eta(i,j)) * dsig(k)
      ! For `VCOORD_ZSTAR_FULL` dsig is a fallback used when the per-
      ! column `z_ref` table is not populated.
      real(wp), allocatable :: dsig(:)
         !! Per-layer σ-fraction.  Sums to 1.0; size `nz_ml`.

      ! ---- Global z-level reference profile ----
      ! Reference z-interfaces in metres (positive-down), `z_ref_global(0)
      ! = 0` is the surface, `z_ref_global(nz_ml)` is the deepest
      ! reference interface.  Drives the z-level branch of `VCOORD_ZSIGMA`
      ! and the z*-lite branch of `VCOORD_ZSTAR_SIGMA`.  Default at init
      ! is uniform 0..1 (normalised) — the namelist parser populates it
      ! with absolute depths when those cases are activated.
      real(wp), allocatable :: z_ref_global(:)
         !! Global reference z-interfaces (m, positive-down), shape
         !! `0:nz_ml`.  Used by ZSIGMA / ZSTAR_SIGMA (NOT by ZSTAR, whose
         !! nominal profile is the `z_fixed` one, `z_fixed_zi`).

      ! ---- Per-column target thickness ----
      ! Recomputed every outer step from (H, eta) per the coord_type
      ! case.  Consumed by the remap kernel (Layer 3) to advance
      ! `multilayer.h_layer` and every `tracers(t)%hTr`.
      real(wp), allocatable :: target_h(:, :, :)
         !! Target layer thickness (m), shape `(nx, ny, nz_ml)`.

      ! ---- Per-column z* reference profile ----
      ! Anchored to local bathymetry for `VCOORD_ZSTAR_FULL`.  Indexed
      ! from `k=0` (surface) to `k=nz_ml` (bed) — opposite of the
      ! bottom-up state convention so the surface anchor is at index 0
      ! (matches the coastal convention in `vcoord_target_dz_column_zstar_full`).
      ! Populated by `build_zref_full(h_bed)`; consumed by the
      ! ZSTAR_FULL branch of `compute_target_h`.
      real(wp), allocatable :: z_ref(:, :, :)
         !! Per-column z* reference (m), shape `(nx, ny, 0:nz_ml)`.

      ! ---- Geopotential depth of the column TOP ----
      ! The rigid-lid seam of the z-like families.  `z_top(i,j)` is the
      ! depth (m, positive down, `>= 0`) of the top of the WATER column
      ! below the `z = 0` datum, i.e. `metrics%z_draft` under an
      ! ice-shelf cavity and identically `0` everywhere else (open
      ! ocean, and every non-cavity run).  Filled once at configure
      ! (`configure_ocean_cavity`) — the draft is static — and never
      ! touched again, so the remap driver's signature is unchanged and
      ! no second static 2-D array is threaded through the `pure` call
      ! chain.
      !
      ! Allocated UNCONDITIONALLY with `source = 0.0_wp`, so it is always
      ! safe to hand to an explicit-shape device dummy: unlike
      ! `metrics%z_draft` (a `(1,1)` placeholder when no cavity is
      ! configured) this is always `(nx_total, ny_total)`.
      !
      ! Consumed by the `VCOORD_Z_FIXED` branch of `compute_target_h`,
      ! which measures its nominal interface depths from `z = 0` and
      ! clips the stack against `z_top`.  `z_top ≡ 0` ⇒ every geometric
      ! branch reproduces its pre-cavity arithmetic bit-for-bit.
      real(wp), allocatable :: z_top(:, :)
         !! Geopotential depth of the column top (m, positive down),
         !! shape `(nx, ny)`.  `0` = the free surface at `z = 0`.

      ! ---- Isopycnal (VCOORD_RHO) target densities ----
      ! Monotone-increasing nominal interface potential densities
      ! (kg/m³) referenced to `rho_ref_pressure`.  Indexed `0:nz_ml`:
      ! `rho_target(0)` is the lightest (surface, maps to the k=nz
      ! interface in bottom-up state); `rho_target(nz_ml)` the densest
      ! (bed, k=1 interface).  The `compute_target_h_rho` inversion
      ! places interior interfaces where the reconstructed column
      ! density equals each interior target.  Sized + populated only
      ! when `coord_type == VCOORD_RHO`; ignored otherwise.
      real(wp), allocatable :: rho_target(:)
         !! Target interface potential densities (kg/m³), shape `0:nz_ml`.

      ! ---- ALE remap workspaces ----
      ! Allocated once at init + enter_data'd to device.  Previously
      ! these were allocated per call inside `ocean_apply_ale_remap_*`
      ! which generated thousands of host-allocated buffers per day
      ! that the device-side DCs in the remap step had to implicitly
      ! transfer back and forth.  Persistent device-mapped scratch
      ! eliminates that overhead entirely.
      real(wp), allocatable :: remap_total_h(:, :)
         !! Column-total h_layer scratch, shape `(nx, ny)`.
      real(wp), allocatable :: remap_h_ref(:, :)
         !! H reference (total_h − bt_eta) scratch, shape `(nx, ny)`.
      real(wp), allocatable :: remap_h_old(:, :, :)
         !! Snapshot of `h_layer` before the remap, shape `(nx, ny, nz)`.
      real(wp), allocatable :: remap_conc_t(:, :, :)
         !! Layer-mean T concentration scratch for the `VCOORD_RHO`
         !! density inversion, shape `(nx, ny, nz)`.  Built from
         !! `hTr / remap_h_old` (vanishing-layer-guarded) once per remap.
      real(wp), allocatable :: remap_conc_s(:, :, :)
         !! Layer-mean S concentration scratch for `VCOORD_RHO`, shape
         !! `(nx, ny, nz)`.

      ! ---- Tuning knobs ----
      integer :: remap_method = REMAP_PPM
         !! ALE remap reconstruction order (REMAP_PCM/PLM/PPM/PPM_H4/PQM).
         !! PQM falls back to PPM for nz < 5 (see `remap_column_pqm`).
      real(wp) :: zstar_h_surf_target = 5.0_wp
         !! Surface-layer thickness anchor for `VCOORD_ZSTAR_FULL` (m).
      real(wp) :: zstar_h_min = 1.0e-4_wp
         !! Bed-side vanishing-layer floor (m).  **Two contracts, picked by
         !! the coordinate family — see `rdb_vcoord :: vcoord_h_min_role`.**
         !!
         !! On the GEOMETRIC families (`VCOORD_ZSTAR_FULL`, `VCOORD_Z_FIXED`)
         !! this is the thickness handed to filler layers that lie BELOW the
         !! local bed.  They hold no water; the floor exists ONLY so
         !! `target_h` is never exactly zero and no h-dividing kernel can
         !! 1/0.  They are MEANT to be classified vanished downstream, so the
         !! default sits deliberately BELOW the D4 skip/merge marker
         !! `H_VANISHED = 1.5e-4` — not by accident, and not a floor in the
         !! `angstrom_h` sense (the D4 taxonomy forbids using `H_VANISHED`
         !! as a positivity floor).  Thinner is also better physics here:
         !! each filler interface carries the full topographic slope, so the
         !! spurious rest PGF transport it drives scales WITH the floor (the
         !! same argument that took `seed_h_layer_uniform_z_impl` off
         !! `2*H_VANISHED`).  `validate_config` warns on a value above
         !! `H_VANISHED` under these families, and refuses a non-positive one.
         !!
         !! On the DENSITY families (`VCOORD_RHO`, `VCOORD_HYCOM`) the
         !! collapsed layers are real layers the inversion squeezed shut
         !! anywhere in the column; they carry tracer mass, so
         !! `compute_target_h_rho_impl` inflates them to
         !! `max(zstar_h_min, 2*H_VANISHED)` to keep them above the remap
         !! drain.  There `zstar_h_min` is additionally the pre-compaction
         !! strip threshold, so a large value is meaningful rather than wrong.
      integer  :: zstar_n_surf = 0
         !! Number of fine near-surface layers for ZSTAR_FULL.  ≤ 0 =
         !! auto-pick (max(1, nz_ml/3)).
      integer  :: zstar_stretching = STRETCH_UNIFORM
         !! Stretching mode for ZSTAR_FULL.  `STRETCH_UNIFORM` (default)
         !! or `STRETCH_LOG` for a geometric near-surface fine zone.
      real(wp) :: zsigma_depth_transition = 200.0_wp
         !! Sigma → z* transition depth (m) for `VCOORD_ZSTAR_SIGMA`.
      real(wp) :: zsigma_blend_width = 100.0_wp
         !! Smoothstep blend width (m) above the transition depth.
      real(wp) :: rho_ref_pressure = 2.0e7_wp
         !! Reference pressure (Pa, default 2e7 = 2000 dbar) for the
         !! potential density that defines the `VCOORD_RHO` coordinate.
         !! A rdb convention (not MOM6-inherited).
      real(wp) :: z_fixed_h_ref = 0.0_wp
         !! Total reference depth (m) for `VCOORD_Z_FIXED` and
         !! `VCOORD_ZSTAR` (the MOM6 z* nominal profile is the `z_fixed`
         !! one, dilated per column — `ocean_vcoord_zstar_target`).  Layer
         !! interfaces sit at `z = k · h_ref / nz_ml` from the surface,
         !! same as MOM6's `COORD_CONFIG = "gprime"` with `MAXIMUM_DEPTH
         !! = h_ref`.  Driver writes from `cfg%ocean%topo%max_depth` at init.
         !! When 0 (default) the `compute_target_h` Z_FIXED branch falls
         !! back to a uniform `H · dsig(k)` target so the path stays
         !! sane in tests that don't explicitly set this knob.
         !! Under a stretched profile (`z_fixed_use_profile`) it is the
         !! profile's total depth, `z_fixed_zi(0)`, and the nominal
         !! interfaces come from `z_fixed_zi` instead of `h_ref/nz`.
      logical :: z_fixed_use_profile = .false.
         !! `&vcoord_nml z_fixed_profile /= "uniform"`: the `VCOORD_Z_FIXED`
         !! (and `VCOORD_ZSTAR`) nominal interfaces come from `z_fixed_zi` /
         !! `z_fixed_dz` (set by `ocean_vcoord_set_z_fixed_profile`) rather than from the uniform
         !! `z_fixed_h_ref/nz_ml`.  Scalar, rides `copyin(this)`.  Default
         !! `.false.` ⇒ the uniform arithmetic, byte-identical.
      real(wp), allocatable :: z_fixed_zi(:)
         !! `VCOORD_Z_FIXED` nominal interface depths (m, positive down,
         !! below `z = 0`), shape `0:nz_ml`, BOTTOM-UP like the state:
         !! `z_fixed_zi(k)` is the TOP interface of layer `k`, so
         !! `z_fixed_zi(nz_ml) = 0` (the surface) and `z_fixed_zi(0)` is the
         !! profile's total depth.  Allocated at init (zeros), so it is never
         !! a placeholder; read only when `z_fixed_use_profile`.
      real(wp), allocatable :: z_fixed_dz(:)
         !! `VCOORD_Z_FIXED` nominal layer thicknesses (m), shape `nz_ml`,
         !! bottom-up: `z_fixed_zi(k-1) - z_fixed_zi(k)`, stored separately
         !! so the partial-top-cell threshold uses the exact namelist
         !! value.  Read only when `z_fixed_use_profile`.
      real(wp) :: regrid_time_scale = 0.0_wp
         !! Grid time-filter timescale τ (s) for the ALE regrid.  After
         !! `compute_target_h` builds the new target grid, the remap step
         !! relaxes the coordinate a fraction `dt/(τ+dt)` toward that
         !! target each outer step rather than jumping to it — damping the
         !! per-step grid-motion shock that drives the σ/z* PGE
         !! (White & Adcroft 2008, the grid time-filter).  Scalar on the
         !! type, reaches the device through the existing `copyin(this)`;
         !! no new device array.  Default `0.0` ⇒ `wtd = 1` ⇒ jump to
         !! target ⇒ bit-identical to the no-filter remap.
      logical :: remap_boundary_extrap = .false.
         !! Close the ALE remap's reconstruction at the two boundary cells
         !! (`k=1`, `k=nz`) with the linear-exact one-sided edge pair
         !! instead of the PCM flatten (MOM6 `BOUNDARY_EXTRAPOLATION`).
         !!
         !! The default closure makes PLM/PPM/PPM_H4/PQM first-order in
         !! exactly the two cells adjacent to the bed and the surface, so
         !! a column whose tracer is linear in z is remapped with an O(h)
         !! error there every thermo step.  Under a terrain-following
         !! coordinate over a slope that error differs between neighbouring
         !! columns, which is a horizontal density gradient, which is a
         !! spurious pressure-gradient force — and with rotation it feeds a
         !! growing grid mode trapped in those same layers (see
         !! `docs/CAPABILITIES_AND_LIMITATIONS.md`).  Scalar on the type,
         !! reaches the device through the existing `copyin(this)`; no new
         !! device array.  Default `.false.` ⇒ bit-identical.
      logical :: remap_nonuniform_weights = .false.
         !! Use the non-uniform-grid reconstruction weights in the ALE
         !! remap's PLM slope and PPM edge estimate — Colella & Woodward
         !! (1984) eqs (1.6)-(1.8) — instead of their equal-thickness
         !! specialisations (`0.5·minmod` and `(7/12, -1/12)`).
         !!
         !! The shipped formulae are linear-exact only when the SOURCE
         !! column is uniform, which under every geometric family but
         !! `sigma`-on-flat-bed it is not: a stretched column carries an
         !! O(Δh/h) reconstruction error on a profile linear in z, in the
         !! whole interior rather than only at the two boundary cells
         !! `remap_boundary_extrap` addresses.  The two knobs are
         !! complementary — the interior needs this one, the outermost two
         !! cells need that one, and a column is exact only with BOTH.
         !! PPM_H4 and PQM carry thickness-weighted stencils already and are
         !! unaffected (their small-`nz` fallbacks excepted).  Scalar on the
         !! type, reaches the device through the existing `copyin(this)`; no
         !! new device array.  Default `.false.` ⇒ bit-identical.
      logical :: remap_check_preconditions = .false.
         !! Assert the ALE remap's column preconditions once per remap and
         !! fail loud on a violation (audit findings V5, V6).
         !!
         !! The overlap sweep every reconstruction shares assumes both
         !! `dz >= 0` (a negative source thickness makes the cumulative
         !! interface stack NON-MONOTONE, and the sweep then integrates the
         !! reversed interval twice — creating mass with no NaN and no bounds
         !! hit) and `sum(dz_old) == sum(dz_new)` (a short target silently
         !! deletes the non-overlapping tail; a long one integrates it as
         !! `q = 0`).  Neither has ever been checked, and the target builders
         !! break the second one on degenerate columns.  Diagnostic — a
         !! per-column reduction at the THERMO cadence, two scalars back to
         !! the host.  Default `.false.` ⇒ the check never runs.
      logical :: remap_vel_conserve_ke = .false.
         !! Enable the KE-conserving rescale of the remapped layer
         !! velocities.  After the per-face column remap (which already
         !! conserves `u·h`, i.e. momentum), rescale the BAROCLINIC
         !! velocity anomaly per column so column KE `Σ ½ h·u²` is
         !! preserved (Adcroft & Hallberg 2006 layer-velocity remap),
         !! capped at a 1.25× rescale factor.  The barotropic/depth-mean
         !! component is never touched (mode-split consistency).  Default
         !! `.false.` ⇒ velocities unchanged ⇒ bit-identical.
      logical :: check_vanished_content = .false.
         !! `&vcoord_nml check_vanished_content` — the I1′ tripwire.  Carried
         !! on this slot (rather than on `ocean_dyn_t`) because the vertical
         !! coordinate is what MAKES vanished layers, so the knob that
         !! polices them belongs beside `zstar_h_min` and the filler
         !! contract.  A plain scalar on the type: it rides the existing
         !! `copyin(this)` and adds no device array.  Read by
         !! `check_vanished_invariant_or_die` in `rdb_ocean_dyn`.  Default
         !! `.false.` ⇒ no scan, no cost.
      logical :: zfixed_closed_faces = .false.
         !! `&vcoord_nml zfixed_closed_faces` — partial-step z-level face
         !! closure.  Meaningful on the three GEOMETRIC families that
         !! vanish bed-side layers, `VCOORD_Z_FIXED`, `VCOORD_ZSTAR` (MOM6
         !! z*) and `VCOORD_ZSTAR_FULL`, where a layer whose reference range lies
         !! inside the bed (or, under `z_fixed`, the ice draft) is an
         !! inert FILLER; a velocity face at which that layer is a filler
         !! on EITHER side is a z-level WALL, not a thin passage (Adcroft,
         !! Hill & Marshall 1997; Losch 2008).
         !!
         !! The per-layer 0/1 face mask itself lives on `ocean_metrics_t`
         !! (`open_u`/`open_v`, built once at configure by
         !! `ocean_vcoord_closed_face_masks` from THIS module's target at
         !! `eta = 0`, `ocean_vcoord_eta0_target`).  The flag is carried here so the ALE
         !! remap driver — which never sees `ocean_metrics_t` — can build
         !! its FACE columns as `min(h_L, h_R)` and drop the closed
         !! layers, instead of pouring momentum into water that is not
         !! there.  Scalar on the type, reaches the device through the
         !! existing `copyin(this)`.  Default `.false.` => bit-identical.

      ! ---- Cached extents (for kernel loops + sanity checks) ----
      integer :: nx_total = 0
         !! Total i-extent of `target_h` (incl. halos).
      integer :: ny_total = 0
         !! Total j-extent of `target_h` (incl. halos).
      integer :: nz_ml = 0
         !! Number of active layers.
   contains
      procedure, non_overridable :: init => ocean_vcoord_init
      procedure, non_overridable :: destroy => ocean_vcoord_destroy
      procedure, non_overridable :: enter_data => ocean_vcoord_enter_data
      procedure, non_overridable :: exit_data => ocean_vcoord_exit_data
      procedure, non_overridable :: compute_target_h => ocean_vcoord_compute_target_h
      procedure, non_overridable :: compute_target_h_rho => ocean_vcoord_compute_target_h_rho
      procedure, non_overridable :: build_zref_full => ocean_vcoord_build_zref_full
      procedure, non_overridable :: bytes => ocean_vcoord_bytes
   end type ocean_vcoord_t