ocean_metrics_t Derived Type

   type :: ocean_metrics_t
      logical :: is_init = .false.
         !! True between `init` and `destroy`.  Guard on this, never on
         !! `allocated(...)` (host pointer only; misses GPU mapping).

      ! ---- Lengths (m) ----
      real(wp), allocatable :: dxT(:, :), dyT(:, :)
         !! Cell-centre (T) zonal/meridional grid lengths (m), `(nx,ny)`.
      real(wp), allocatable :: dxCu(:, :), dyCu(:, :)
         !! u-face (Cu) lengths (m), `(nx+1,ny)`.
      real(wp), allocatable :: dxCv(:, :), dyCv(:, :)
         !! v-face (Cv) lengths (m), `(nx,ny+1)`.
      real(wp), allocatable :: dxBu(:, :), dyBu(:, :)
         !! Corner (Bu) lengths (m), `(nx+1,ny+1)`.

      ! ---- Topography-aware face widths (m) ----
      real(wp), allocatable :: dy_cu(:, :)
         !! Open zonal width of the u-face for transport (m), `(nx+1,ny)`.
         !! v1: filled = `dyCu` but a SEPARATE array, so the transport
         !! kernels read the right name once porous/partial cells arrive.
      real(wp), allocatable :: dx_cv(:, :)
         !! Open meridional width of the v-face for transport (m),
         !! `(nx,ny+1)`.  v1: filled = `dxCv`, separate array.

      real(wp), allocatable :: dy_cu_bt(:, :)
         !! Open zonal u-face width the BAROTROPIC substep transports on
         !! (m), `(nx+1,ny)`.  ALWAYS full size and byte-equal to `dy_cu`
         !! unless porous barriers are on, in which case the per-step
         !! refresh scales it by the COLUMN-INTEGRATED open fraction
         !! `(A(eta_top) - A(eta_bed)) / (eta_top - eta_bed)`, which is
         !! identically the THICKNESS-WEIGHTED MEAN of the per-layer
         !! fractions (both are the same integral of `w` over the column,
         !! so the identity is exact, not an approximation).  Without it
         !! the barotropic solve would be porous-blind and the layer
         !! renormalisation (which drives `sum_k flux_k = uhbt`) would
         !! hand the blocked transport straight back.
         !!
         !! NOT a claim of `BT_cont` parity.  MOM6's production barotropic
         !! face area is `sum_k (dy_Cu*por_k) * h_marginal_k * visc_rem_k`
         !! — weighted by the PPM MARGINAL thickness and by `visc_rem`,
         !! neither of which appears here; its `sum_k h_k*(dy_Cu*por_k)`
         !! form is the open-boundary-segment branch only, and its
         !! `set_local_BT_cont_types` carries no `por` at all.  What this
         !! array reproduces is the telescoping identity above, applied to
         !! the plain layer thicknesses.
      real(wp), allocatable :: dx_cv_bt(:, :)
         !! v-face twin, `(nx,ny+1)`.

      ! ---- Porous barriers (Adcroft 2013; see rdb_ocean_porous) ----
      logical :: use_porous = .false.
         !! Master switch (`&ocean_porous_nml enable`).  OFF ⇒ the
         !! `por_face_area_*` arrays stay at their `(1,1,1)` placeholder
         !! size and every transport kernel takes the un-narrowed
         !! `dy_cu` / `dx_cv` branch — byte-identical to a build without
         !! porous barriers.
      integer :: porous_eta_interp = 0
         !! Interface-at-velocity-point rule, a `POROUS_ETA_*` value
         !! (`rdb_ocean_porous`).  0 = MAX (the shallower interface).
      real(wp) :: porous_mask_depth = 0.0_wp
         !! Gate HEIGHT (m, positive up, `<= 0`): faces whose mean
         !! along-face height is at or above this stay fully open
         !! (MOM6 `PORBAR_MASKING_DEPTH`, sign-flipped to a height).
      real(wp), allocatable :: por_bed(:, :)
         !! Static snapshot of the bottom topographic HEIGHT at cell
         !! centres (m, positive up — i.e. `-barotropic%b`, which is the
         !! positive-down reference depth), `(nx,ny)` when `use_porous`,
         !! `(1,1)` otherwise.
         !! The porous curve works in ABSOLUTE heights, so the recompute
         !! needs the bed on the same datum as `por_d*`; keeping a copy
         !! here makes the kernel self-contained (no barotropic-state
         !! argument threaded through the dynamics).
      real(wp), allocatable :: por_dmin_u(:, :), por_dmax_u(:, :), por_davg_u(:, :)
         !! u-face along-face deepest / shallowest / mean topographic
         !! height (m, positive up), `(nx+1,ny)` when `use_porous`,
         !! `(1,1)` otherwise.  Static — filled once at setup.
      real(wp), allocatable :: por_dmin_v(:, :), por_dmax_v(:, :), por_davg_v(:, :)
         !! v-face twins, `(nx,ny+1)` when `use_porous`, `(1,1)` otherwise.
      real(wp), allocatable :: por_face_area_u(:, :, :)
         !! u-face layer-averaged OPEN-AREA fraction (nondim, `[0,1]`),
         !! `(nx+1,ny,nz)` when `use_porous`, `(1,1,1)` otherwise.
         !! Recomputed every RK2 stage (interface-height dependent) and
         !! MULTIPLIED into `dy_cu` by the transport kernels.
      real(wp), allocatable :: por_face_area_v(:, :, :)
         !! v-face twin, `(nx,ny+1,nz)` when `use_porous`, `(1,1,1)`
         !! otherwise.

      ! ---- Partial-step z-level face closure (VCOORD_Z_FIXED) ----
      logical :: use_closed_faces = .false.
         !! Master switch (`&vcoord_nml zfixed_closed_faces`), latched by
         !! `configure_ocean_closed_faces`.  OFF ⇒ `open_u`/`open_v` stay
         !! at their `(1,1,1)` placeholder size, no kernel branch is
         !! taken, byte-identical to a build without the feature.
      real(wp), allocatable :: open_u(:, :, :)
         !! u-face per-layer 0/1 OPEN mask, `(nx+1,ny,nz)` when
         !! `use_closed_faces`, `(1,1,1)` otherwise.  1 = the layer has
         !! water on BOTH sides of the face; 0 = it is an inert `z_fixed`
         !! filler on at least one side and the face is a z-LEVEL WALL for
         !! that layer (Adcroft, Hill & Marshall 1997; Losch 2008).
         !! STATIC — built once at configure by
         !! `ocean_vcoord_closed_face_masks` from the `z_fixed` target at
         !! `η = 0`, never refreshed (the bed and the draft are static and
         !! `η` is absorbed by the first live layer).
      real(wp), allocatable :: open_v(:, :, :)
         !! v-face twin, `(nx,ny+1,nz)` when `use_closed_faces`,
         !! `(1,1,1)` otherwise.
         !!
         !! ### THE COMPOSITION RULE (stated once, here)
         !!
         !! The three face gates are INDEPENDENT and compose by
         !! multiplication — none replaces another:
         !! ```
         !! dy_eff(I,j,k) = dy_cu(I,j) · por_face_area_u(I,j,k) · open_u(I,j,k)
         !! dx_eff(i,J,k) = dx_cv(i,J) · por_face_area_v(i,J,k) · open_v(i,J,k)
         !! ```
         !! `dy_cu`/`dx_cv` carry the 2-D LAND decision (metric zeroing in
         !! `metrics_apply_land_mask`); `por_face_area_*` narrows
         !! continuously for unresolved SUBGRID sills (Adcroft 2013); and
         !! `open_*` closes per LAYER for the resolved z-level staircase.
         !! Porous barriers and closed faces are therefore NOT mutually
         !! exclusive.
         !!
         !! Every consumer applies the two 3-D factors as SEPARATE,
         !! separately host-gated, INLINE `do concurrent` passes rather
         !! than pre-composing them into a third array.  Two reasons:
         !! a composed array would have to be recomputed whenever the
         !! porous fit is refreshed (per outer step) and so could not be
         !! static; and an inert host-gated branch that never names the
         !! array costs nothing, whereas handing a state array to an
         !! external helper pessimises every `do concurrent` in the
         !! calling routine even when the branch is not taken (CLAUDE.md,
         !! measured at +4.8 % for an inert porous pass).

      ! ---- Static ice-shelf cavity geometry (P5.1; see rdb_ocean_cavity) ----
      logical :: use_cavity = .false.
         !! Master switch (`&ocean_cavity_dyn_nml enable`), latched in
         !! `ocean_state_init_from_config` BEFORE `init` so the allocation
         !! gate below can read it.  OFF ⇒ `z_draft` / `cover_frac` /
         !! `p_ice_ref` stay at their `(1,1)` placeholder size, `bt_H_ref`
         !! latches the bed as it always did, and every path is
         !! byte-identical to a build without cavities.
      real(wp), allocatable :: z_draft(:, :)
         !! Prescribed STATIC ice-base depth (m, positive DOWN, `>= 0`),
         !! `(nx_total, ny_total)` INCLUDING ghosts when `use_cavity`,
         !! `(1,1)` otherwise.  Filled by the formula setters in
         !! `rdb_ocean_cavity` immediately after the bathymetry, then
         !! carried through the SAME periodic/fold re-wrap + halo sequence
         !! `barotropic%b` gets (ordering is load-bearing: the draft must
         !! exist before the wet mask is seeded from `b - z_draft`).
         !! `z_draft = 0` is open ocean — including beyond the calving
         !! front.
      real(wp), allocatable :: cover_frac(:, :)
         !! Ice-covered area fraction (nondimensional, `[0,1]`), same
         !! shape + gating as `z_draft`.  v1 is BINARY, `merge(1, 0,
         !! z_draft > 0)`; an area-blended calving front belongs to the
         !! melt work.  Allocated alongside `z_draft` so the thermodynamic
         !! slice does not have to re-open this lifecycle.
      real(wp), allocatable :: p_ice_ref(:, :)
         !! Boussinesq-isostatic (flotation) ice load `rho_ref*GRAVITY*
         !! z_draft` (Pa, `>= 0`), same shape + gating as `z_draft`.  Built
         !! ONCE at configure from the SAME product the FV_MOM6 surface BC
         !! forms (`rho_ref*GRAVITY`), which is what makes
         !! `pa(nz+1) = rho_ref*g*eta_geo + p_ice_ref` cancel to bit-zero
         !! at rest.  Stored rather than recomputed so `GRAVITY`/`rho_ref`
         !! cannot drift between the two users.  CONSUMED as the static
         !! half of `multilayer_state_t%p_top = p_ice_ref + sf%p_surf` —
         !! seeded in `configure_ocean_cavity` and rebuilt each outer step
         !! in `ocean_dyn_step_split` whenever the psurf seam makes
         !! `sf%p_surf` live.  It is the load's route into the PRESSURE
         !! (the FV_MOM6 `pa(nz+1)` top BC and the in-situ EOS); its route
         !! into the BAROTROPIC mode is the datum `bt_H_ref = b - z_draft`
         !! and nothing else, which is why it never joins `sf%p_surf`.

      ! ---- Static land masks (real 0/1; derived in metrics_apply_land_mask) ----
      real(wp), allocatable :: wet_T(:, :)
         !! T-cell wet (1) / land (0) mask, `(nx,ny)` — the HALO-VALID
         !! working copy of `multilayer%wet_mask` (periodic-wrapped +
         !! north-folded, R5a) that `wet_u/wet_v/wet_q` are derived from.
         !! Kept device-resident so the continuity + tracer PPM
         !! reconstruction can mirror a land neighbour's thickness to the
         !! local cell (spec §14 C2 / MOM6's reflected-coast PPM).
         !! All-wet domain ⇒ `wet_T≡1` ⇒ mirror never triggers (no-op).
      real(wp), allocatable :: wet_u(:, :)
         !! u-face (Cu) open mask, `(nx+1,ny)`.  `wet_u(i,j) =
         !! wet_T(i-1,j)*wet_T(i,j)` — a u-face is open iff BOTH adjacent
         !! T-cells are wet.  `mass_flux_x(i,j)` is the west face of cell
         !! `(i,j)` (continuity divergence reads `flux(i+1)-flux(i)`), so
         !! the `i-1`/`i` pairing matches `dy_cu`'s stagger exactly.
         !! A face of ZERO width (`dy_cu = 0`, the tripolar cap's
         !! node-aligned pole columns) is closed too, wet neighbours or
         !! not — see `metrics_apply_land_mask`.
         !! All-wet domain without such faces ⇒ `wet_u≡1` ⇒ masking is a
         !! literal no-op.
      real(wp), allocatable :: wet_v(:, :)
         !! v-face (Cv) open mask, `(nx,ny+1)`.  `wet_v(i,j) =
         !! wet_T(i,j-1)*wet_T(i,j)`, and 0 on a zero-width face
         !! (`dx_cv = 0`).
      real(wp), allocatable :: wet_q(:, :)
         !! Corner (Bu) open mask, `(nx+1,ny+1)`.  Free-slip product of
         !! the 4 surrounding T-cells: `wet_q(i,j) =
         !! wet_T(i-1,j-1)*wet_T(i,j-1)*wet_T(i-1,j)*wet_T(i,j)`.  Consumed
         !! by the relative-vorticity / strain factor (CHUNK B).

      ! ---- Areas (m^2) — load-bearing; `dx*dy` is dead (D5) ----
      real(wp), allocatable :: areaT(:, :)
         !! T-cell area (m^2), `(nx,ny)`.
      real(wp), allocatable :: areaCu(:, :)
         !! Cu-cell area (m^2), `(nx+1,ny)`.
      real(wp), allocatable :: areaCv(:, :)
         !! Cv-cell area (m^2), `(nx,ny+1)`.
      real(wp), allocatable :: areaBu(:, :)
         !! Bu-cell area (m^2), `(nx+1,ny+1)`.

      ! ---- Stored inverses (Adcroft reciprocal; filled in finalize) ----
      real(wp), allocatable :: idxT(:, :), idyT(:, :)
         !! 1/dxT, 1/dyT (1/m), `(nx,ny)`.
      real(wp), allocatable :: idxCu(:, :), idyCu(:, :)
         !! 1/dxCu, 1/dyCu (1/m), `(nx+1,ny)`.
      real(wp), allocatable :: idxCv(:, :), idyCv(:, :)
         !! 1/dxCv, 1/dyCv (1/m), `(nx,ny+1)`.
      real(wp), allocatable :: iareaT(:, :)
         !! 1/areaT (1/m^2), `(nx,ny)`.
      real(wp), allocatable :: iareaBu(:, :)
         !! 1/areaBu (1/m^2), `(nx+1,ny+1)`.
      real(wp), allocatable :: iareaCu(:, :)
         !! 1/areaCu (1/m^2), `(nx+1,ny)`.
      real(wp), allocatable :: iareaCv(:, :)
         !! 1/areaCv (1/m^2), `(nx,ny+1)`.

      ! ---- Geography (degrees) ----
      real(wp), allocatable :: geolatT(:, :), geolonT(:, :)
         !! Latitude / longitude at T points (degrees), `(nx,ny)`.
      real(wp), allocatable :: geolatBu(:, :), geolonBu(:, :)
         !! Latitude / longitude at Bu corners (degrees), `(nx+1,ny+1)`.
      real(wp), allocatable :: angle_dx(:, :)
         !! Grid ROTATION at T points (RADIANS), `(nx,ny)`: the angle of the
         !! grid's +i axis measured COUNTER-CLOCKWISE from true east — MOM6's
         !! `angle_dx` convention (the mosaic stores it in degrees at every
         !! supergrid node; the T value is node `(2i,2j)`).  It rotates a
         !! geographic (east, north) vector onto the grid axes:
         !!
         !!     u_grid =  cos(angle_dx)*u_east + sin(angle_dx)*v_north
         !!     v_grid = -sin(angle_dx)*u_east + cos(angle_dx)*v_north
         !!
         !! and back with the transpose.  This is how lat-lon vector forcing
         !! (e.g. wind stress on a reanalysis grid) is put on a curvilinear
         !! grid — MOM6 does the same with `G%cos_rot` / `G%sin_rot`.
         !! Zero on Cartesian and spherical grids (the axes ARE east/north).
         !! Supergrid: read from the file's `angle_dx`, else (and on the
         !! analytic tripolar) derived from the node geography by
         !! `supergrid_angle_dx_from_geography`.  Ghosts: extrapolated,
         !! then wrapped / folded like every other metric — across the fold
         !! the conjugate cell's +i axis points the OTHER way, so the folded
         !! ghost rows carry `angle + pi`.  Static; no kernel reads it yet
         !! (the forcing regridder will).

      ! ---- hvisc ratio bundle (dimensionless / m; filled in finalize) ----
      real(wp), allocatable :: dy_dxT(:, :)
         !! dyT/dxT at T (dimensionless), `(nx,ny)`.  =1 on Cartesian.
      real(wp), allocatable :: dx_dyT(:, :)
         !! dxT/dyT at T (dimensionless), `(nx,ny)`.
      real(wp), allocatable :: dy_dxBu(:, :)
         !! dyBu/dxBu at Bu (dimensionless), `(nx+1,ny+1)`.
      real(wp), allocatable :: dx_dyBu(:, :)
         !! dxBu/dyBu at Bu (dimensionless), `(nx+1,ny+1)`.
      real(wp), allocatable :: dx2h(:, :)
         !! dxT^2 at T (m^2), `(nx,ny)`.
      real(wp), allocatable :: dy2h(:, :)
         !! dyT^2 at T (m^2), `(nx,ny)`.
      real(wp), allocatable :: dx2q(:, :)
         !! dxBu^2 at Bu (m^2), `(nx+1,ny+1)`.
      real(wp), allocatable :: dy2q(:, :)
         !! dyBu^2 at Bu (m^2), `(nx+1,ny+1)`.
   contains
      procedure, non_overridable :: init => ocean_metrics_init
      procedure, non_overridable :: destroy => ocean_metrics_destroy
      procedure, non_overridable :: enter_data => ocean_metrics_enter_data
      procedure, non_overridable :: exit_data => ocean_metrics_exit_data
      procedure, non_overridable :: bytes => ocean_metrics_bytes
   end type ocean_metrics_t