ocean_pressure_force_t Derived Type

   type :: ocean_pressure_force_t
      logical :: is_init = .false.
         !! True between `init` and `destroy`.  Prefer this to
         !! `allocated(...)` — tracks GPU device attachment too.
      integer :: variant = OPGF_VARIANT_MONT
         !! Active pressure-force scheme variant.  MONT is the default
         !! here AND the `&ocean_pgf_nml form` default, so a state built
         !! straight from this type and one built through the config
         !! agree.  MONT is general-purpose: valid over variable
         !! bathymetry and every vcoord.
      logical :: use_eos_along_path = .false.
         !! Convenience flag — re-evaluate density at each integration
         !! sample rather than at layer centres.  Phase 5d.
      logical :: scratch_gated = .false.
         !! ALLOCATION GATE for the variant-specific scratch buffers.
         !!
         !! `.false.` (default): `init` allocates all 16 buffers, whatever
         !! `variant` is — the historical behaviour every direct
         !! `pgf%init(grid, nz_ml=...)` call site (tests, benchmarks) relies
         !! on, since those set `variant` AFTER `init`.
         !!
         !! `.true.`: `init` allocates ONLY the buffers the active
         !! `variant` / `reconstruct_for_pressure` actually touch (see the
         !! per-buffer gate comments in `ocean_pressure_force_init`), saving
         !! up to 11 of the 16 (~3.5 GB at 1000x800x50).  Latched together
         !! with `variant` + `reconstruct_for_pressure` BEFORE `init(grid)`
         !! by `ocean_state_init_from_config` — the same conditional-
         !! allocation contract as the default-off closures and
         !! `continuity_t%windowed_advection`.
         !!
         !! Safety contract: every gated-off buffer is unreachable on its
         !! gated-off path because `ocean_pressure_force_compute` dispatches
         !! on `variant` at the TOP and `return`s out of each branch, and the
         !! one external `pgf%e_face` consumer (`compute_pbce`) `error stop`s
         !! unless `variant == OPGF_VARIANT_FV_MOM6`.  `enter_data` /
         !! `exit_data` / `bytes` all key off `allocated(...)`, so a gated-off
         !! buffer is neither mapped nor counted.
      real(wp) :: rho0 = 1035.0_wp
         !! BOUSSINESQ reference density (kg/m³) — the `ρ₀` that divides the
         !! pressure gradient into an acceleration, `du/dt = −(1/ρ₀)·∂p/∂x`.
         !! Read by EVERY variant (it is `inv_rho0` in the face passes and
         !! `g_over_rho0` in the Montgomery recursion), and by
         !! `compute_pbce` in the barotropic coupling.
         !!
         !! ASSIGNED FROM CONFIG by `configure_ocean_pgf`
         !! (`&ocean_ic_nml rho_0` → `eos%rho0`, the single ρ₀ of record —
         !! the EOS, the `eta_ib` surface-pressure seam, EPBL, kappa-shear,
         !! tidal mixing, MEKE/GM and the isopycnal slopes all take the same
         !! scalar).  The literal here is only the pre-configure default for
         !! direct `pgf%init(...)` call sites (tests, benchmarks) that never
         !! run the configure pass; it matches the `&ocean_ic_nml rho_0`
         !! default so a state built either way agrees.
         !!
         !! Host scalar: every read is host-side (into a local `inv_rho0`,
         !! or passed by value into a `*_impl`), so `configure_ocean_pgf`
         !! needs no `!$acc update device` — nothing reads it through the
         !! device-mapped `pgf` handle.

      ! ---- gprime / reduced-gravity knobs (OPGF_VARIANT_GPRIME) ----
      real(wp) :: gprime_gfs = 9.81_wp
         !! Free-surface reduced gravity (m/s²). Full physical g reduces to
         !! a standard free-surface model; reducing it slows the BT mode
         !! (wave speed `sqrt(gfs · H)`) for a slower BT CFL.
      real(wp) :: gprime_gint = 0.0098_wp
         !! Internal reduced gravity (m/s²) at the single layer-1/layer-2
         !! interface. Only used by `OPGF_VARIANT_GPRIME`.

      ! ---- FV_MOM6 reference state (OPGF_VARIANT_FV_MOM6) ----
      real(wp) :: rho_ref = 1035.0_wp
         !! ANOMALY reference density (kg/m³) subtracted from layer densities
         !! when building the `pa` pressure-anomaly stack (FV_MOM6), and the
         !! surface-layer `g·ρ_ref/ρ₀` in `compute_pbce`.
         !! `rho_ref = rho0 = ρ_surf` makes the surface layer's anomaly
         !! vanish; only denser bed layers contribute.
         !!
         !! A DISTINCT ROLE from `rho0` — `rho0` scales the gradient into an
         !! acceleration, `rho_ref` only shifts the baseline the anomaly is
         !! measured from — kept as a separate member so the two are never
         !! silently interchanged (MOM6 carries the same pair, `GV%Rho0` vs
         !! the `PressureForce_FV` `rho_ref`).  Both are nonetheless SOURCED
         !! FROM THE SAME CONFIGURED ρ₀ by `configure_ocean_pgf`
         !! (`&ocean_ic_nml rho_0` → `eos%rho0`): roundabout has no separate
         !! anomaly-reference knob, and `rho_ref ≠ rho0` would put a constant
         !! `g·(ρ_ref−ρ₀)/ρ₀` offset in `pa(top)` that the Boussinesq
         !! divisor no longer cancels.  Same host-scalar note as `rho0` —
         !! no device update needed.
      real(wp) :: h_neglect = 1.0e-10_wp
         !! Face-thickness floor in the FV_MOM6 denominator
         !! `(h_L + h_R + h_neglect)`. Prevents division by zero when both
         !! adjacent cells have vanishing bed layers.
      logical :: skip_nonoverlap = .false.
         !! Grounded-layer PGF gate (`&ocean_isopycnal_nml pgf_skip_nonoverlap`,
         !! set by the driver ONLY under VCOORD_LAGRANGIAN).  When `.true.` the
         !! face PGF is zeroed wherever the layer's z-extents in the two
         !! abutting columns do NOT overlap — a layer that has wedged out
         !! against the bed on one side, where the two-point Jacobian
         !! `Δp_centre + g·ρ_layer·Δz_centre` has no common depth to difference
         !! across and leaves `g·(ρ_layer − ρ̄_ambient)·∂z/∂x` of acceleration on
         !! a RESTING ocean.  `.false.` ⇒ the gate branch is never taken ⇒
         !! bit-identical.
         !!
         !! Honoured by `mont`, `fv_lite`, `fv_wright` and `fv_mom6`.  The gate
         !! reads `z_centre`; `ocean_pressure_force_init`'s allocation gate
         !! provides that buffer unconditionally for MONT and the two
         !! FV_LITE-family variants (all three also use it in the face passes)
         !! and, for FV_MOM6, exactly when this flag is set — so it MUST be
         !! latched before `init` (the ocean path does that in
         !! `ocean_state_init_from_config`, and `configure_ocean_pgf`
         !! re-checks that the latch did not drift).  `gprime` differences
         !! interface positions directly and has no such Jacobian, so it is
         !! the one variant left N/A (the driver warns).
         !!
         !! `mont` needs the gate for the SAME geometric reason the FV forms
         !! do, even though its face expression is not a two-point Jacobian:
         !! a grounded layer sits at the bed on the shallow side and at its
         !! flat-isopycnal height on the deep side, so the two `e_edge`
         !! values entering the `M` recursion are hundreds of metres apart
         !! and `M` stops being horizontally uniform at rest.
      real(wp) :: nonoverlap_vanish_tol = 2.0_wp*H_VANISHED
         !! Thickness (m) at or below which a layer counts as GROUNDED for the
         !! `skip_nonoverlap` gate: a face is zeroed only where the layer's
         !! z-extents do not overlap AND the layer is this thin on at least
         !! one side.  Non-overlap alone is not grounding — a layer that is
         !! MASSIVE on both sides but sits at different depths (a
         !! sigma-seeded stack over a bathymetric step, where a 10 m layer at
         !! 90-100 m faces a 25 m layer at 225-250 m) is a steep coordinate
         !! surface the FV forms are built for, exactly as under
         !! `vcoord_type="sigma"`.  Zeroing its PGF while continuity keeps
         !! moving its mass across the face breaks the PGF-work / PE
         !! exchange and grows energy without bound (the compat-matrix
         !! staircase, MaxCFL panic at step 178).  Set by the driver from
         !! `nonoverlap_vanish_tol_for(angstrom_h)`; the default matches
         !! `angstrom_h = 0`.
      logical :: mass_weight = .false.
         !! FV_MOM6 shelf-break `hWght` mass-weighting toggle. When `.true.`
         !! Pass-2's horizontal pressure integral biases the face density
         !! toward the thinner column at unequal-depth faces, cancelling the
         !! spurious bottom-layer shelf-break PGF. Equal-depth columns ⇒
         !! reduces exactly to the midpoint average ⇒ bit-identical.
      logical :: reconstruct_for_pressure = .false.
         !! FV_MOM6 in-layer T/S reconstruction toggle. `.false.` (default):
         !! layer-mean (PCM) density ⇒ bit-identical. `.true.`: per-layer
         !! `dpa` / `intz_dpa` from a 5-point Boole quadrature of a monotone
         !! PLM/PPM sub-layer T/S profile (Adcroft, Hallberg & Harrison 2008;
         !! White, Adcroft & Hallberg 2009), removing the spurious-PGF error
         !! on thick/sloped layers. Only consulted by FV_MOM6.
      integer :: recon_scheme = PGF_RECON_PLM
         !! In-layer reconstruction scheme: 1 = PLM, 2 = PPM. Only consulted
         !! when `reconstruct_for_pressure = .true.`.
      logical :: insitu_density = .true.
         !! FV_MOM6 constant-by-layer (PCM) density at its IN-SITU pressure
         !! (`&ocean_pgf_nml insitu_density`, MOM6 parity).  `.true.`
         !! (default): each layer's `dpa`/`intz_dpa` and the cross-face
         !! `intx_dpa`/`inty_dpa` integrate `EOS(T, S, p = -g·rho0·z)` with
         !! the layer-mean T/S (`compute_fv_mom6_insitu_pcm_impl`).  Under
         !! Wright the vertical integral is ANALYTIC (MOM6
         !! `int_density_dz_wright`: one polynomial evaluation per layer and
         !! per lateral sub-column); under Roquet the 5-point Boole rule of
         !! MOM6 `int_density_dz_generic_pcm` with the (T, S) part of the
         !! EOS hoisted out of the pressure points (`roquet_pcm_dpa_intz`: 1
         !! T/S polynomial + 5 pressure Horners per layer, 3 + 15 per face).
         !! `.false.`: the legacy PCM integral of `ms%rho_layer`, a
         !! POTENTIAL density at the single `&ocean_eos_nml p_ref`, which
         !! drops the pressure dependence of the horizontal density
         !! gradient below the reference level.  Consulted only by FV_MOM6
         !! with `reconstruct_for_pressure = .false.`, an EOS handle and
         !! T/S, and only for a PRESSURE-DEPENDENT EOS (Wright, Roquet):
         !! for the linear EOS in-situ and potential density coincide, so
         !! the legacy path runs and answers are bit-identical.
      logical :: p_top_in_bc = .false.
         !! FV_MOM6 top-of-column pressure in the surface boundary
         !! condition (`&ocean_pgf_nml p_top_in_bc`). `.false.` (default):
         !! `pa(nz+1) = rho_ref·g·eta_geo`, bit-identical. `.true.`: the
         !! load `multilayer_state_t%p_top` (Pa) is ADDED there, so the
         !! pressure stack measures down from the loaded surface —
         !! `pa(nz+1) = rho_ref·g·eta_geo + p_top`. Only consulted by
         !! FV_MOM6 (both the PCM and the `reconstruct_for_pressure`
         !! branch); fail-loud at configure for any other variant, which
         !! carries no injectable `pa` stack.
         !!
         !! A DEPTH-UNIFORM `p_top` perturbs every layer's `PFu` by the
         !! SAME `−(1/ρ₀)∇p_top` (Theorem 1 in the Pass-1 docstring
         !! below), and the split solver replaces the depth mean of the
         !! layer PGF with the barotropic solution, so the baroclinic
         !! operator does not see it and this does NOT double-count the
         !! `eta_forcing` seam. See the `p_top` seam contract in
         !! `src/core/ocean/README.md`.
      real(wp) :: gfs_scale = 1.0_wp
         !! Free-surface gravity scaling (= GFS / G_EARTH). Default 1.0 ⇒
         !! pure FV_MOM6, bit-identical. When < 1, Pass 5 applies the
         !! Montgomery `dM` correction subtracting
         !! `(1 - gfs_scale)·(g/ρ₀)·ρ_surf·∇η` from every layer's PGF
         !! (depth-independent, so the BT mass-flux invariant survives); the
         !! driver also drops `bt_work%g_bt` to `gfs_scale·GRAVITY`.
         !! Combined ⇒ wave speed `sqrt(gfs_scale · g · H)`.

      ! ---- Bathymetry (copy of state%barotropic%b) ----
      ! Used by the gprime PGF to recover ∇η from ∇(sum h_layer) - ∇b;
      ! without it variable-bathymetry runs see a spurious PGF dominated
      ! by ∇H_bathy (~10⁴× larger than ∇η at shelf-breaks). Default zero =
      ! flat bed. Set by the driver via `set_bathymetry()`.
      real(wp), allocatable :: b(:, :)

      ! ---- Workspace ----
      type(scratch_3d_buffer_t) :: p_edge
         !! Layer-edge pressure stack at cell centres.  Shape
         !! (nx, ny, nz_ml+1).  k=1 bed, k=nz_ml+1 surface.
         !! Reused across the two horizontal-gradient passes.
      type(scratch_3d_buffer_t) :: z_centre
         !! Layer-centre z at cell centres, from the free surface downward
         !! (z=0 surface, negative below). Shape (nx, ny, nz_ml). Filled for
         !! MONT and the FV_LITE-family branches. Surface-relative (NOT
         !! bed-relative) so the σ-coord Jacobian cancels the spurious
         !! cross-bathymetry pressure gradient — and so that neither form
         !! double-counts the barotropic `-g·∇η` the BT substep already
         !! carries. MONT reads it as the interface height
         !! `e_edge(k) = z_centre(k) + 0.5·h_layer(k)`.
      type(scratch_3d_buffer_t) :: mont_M
         !! Boussinesq Montgomery potential `M` at layer centres (m² s⁻²).
         !! Shape (nx, ny, nz_ml). Written by the MONT column recursion and
         !! read by its two face passes; no other variant touches it.
      type(scratch_3d_buffer_t) :: rho_insitu
         !! In-situ layer-centre density from the Wright
         !! column-sweep Picard step.  Shape (nx, ny, nz_ml).
         !! Populated only by the FV_WRIGHT branch; the other
         !! variants source ρ from `ms%rho_layer`.
      type(scratch_3d_buffer_t) :: dpdx_face
         !! East-face PGF acceleration: -(1/rho0) * dp/dx at
         !! u-face.  Shape (nx+1, ny, nz_ml).  Apply adds dt*dpdx
         !! to u_face_x_layer.
      type(scratch_3d_buffer_t) :: dpdy_face
         !! North-face PGF acceleration.  Shape (nx, ny+1, nz_ml).

      ! ---- FV_MOM6 scratch (OPGF_VARIANT_FV_MOM6) ----
      !! Interface heights `e` (positive-up).  `e(:,:,1)` is the bed
      !! (= -b), `e(:,:,nz+1)` is the free surface (= η).  Same
      !! sign convention as MOM6, but bottom-up indexing to match the
      !! Roundabout bed-up layer convention.
      type(scratch_3d_buffer_t) :: e_face
         !! Interface heights at cell centres, shape (nx, ny, nz+1).
      !! Pressure anomaly stack at interfaces, units Pa.  Relative to
      !! the `rho_ref · g · z` baseline so the surface-pressure
      !! contribution is just `pa(top) = rho_ref · g · η`.  Built by
      !! marching down from the surface; `pa(k) − pa(k+1) = dpa(k)`
      !! where `dpa(k) = (Rlay(k) − rho_ref) · g · h(k)`.
      type(scratch_3d_buffer_t) :: pa
         !! Pressure anomaly at interfaces, shape (nx, ny, nz+1).
      !! Per-layer vertical integral of `dpa` from layer top inward.
      !! For Boussinesq Rlay path: `intz_dpa(k) = 0.5 · (Rlay(k) −
      !! rho_ref) · g · h(k)²` (mid-point rule).
      type(scratch_3d_buffer_t) :: intz_dpa
         !! Per-layer ∫ dpa dz, shape (nx, ny, nz), units Pa·m.
      ! Horizontal integrals at faces — average of the two adjacent
      ! cells.  `intx_pa(K) = 0.5·(pa_L + pa_R)` at interface K;
      ! `intx_dpa(k) = 0.5·(Rlay(k) − rho_ref) · g · (h_L + h_R)`
      ! within layer k.  Surface BC fixes intx_pa at top; deeper
      ! values come from `intx_pa(K+1) = intx_pa(K) + intx_dpa(k)`.
      type(scratch_3d_buffer_t) :: intx_pa
         !! u-face ∫ pa across x, shape (nx+1, ny, nz+1), units Pa·m.
      type(scratch_3d_buffer_t) :: inty_pa
         !! v-face ∫ pa across y, shape (nx, ny+1, nz+1), units Pa·m.
      type(scratch_3d_buffer_t) :: intx_dpa
         !! u-face ∫ dpa, shape (nx+1, ny, nz), units Pa·m.
      type(scratch_3d_buffer_t) :: inty_dpa
         !! v-face ∫ dpa, shape (nx, ny+1, nz), units Pa·m.
      type(scratch_3d_buffer_t) :: conc_T
         !! Layer-mean temperature the in-situ PCM and reconstruction
         !! kernels read, shape (nx, ny, nz): `hT/h` on a live layer, the
         !! I1′ donor's on a vanished one (Pass C of both kernels).  Filled
         !! once per call by a per-column pass, so no later pass divides by
         !! a thickness or walks a column.  Unused by the `rho_layer` path.
      type(scratch_3d_buffer_t) :: conc_S
         !! Layer-mean salinity, as `conc_T`.

      ! ---- In-layer reconstruction scratch (reconstruct_for_pressure) ----
      !! Per-column PLM/PPM top (shallower) and bottom (deeper) edge
      !! values of the layer-mean T and S, shape (nx, ny, nz).  Filled
      !! by `compute_fv_mom6_reconstruct_impl`'s edge-build pass; consumed
      !! by the 5-point Boole quadrature.  Allocated only when
      !! `reconstruct_for_pressure` (or `scratch_gated = .false.`).
      type(scratch_3d_buffer_t) :: recon_T_t
      type(scratch_3d_buffer_t) :: recon_T_b
      type(scratch_3d_buffer_t) :: recon_S_t
      type(scratch_3d_buffer_t) :: recon_S_b
   contains
      procedure, non_overridable :: init => ocean_pressure_force_init
      procedure, non_overridable :: destroy => ocean_pressure_force_destroy
      procedure, non_overridable :: enter_data => ocean_pressure_force_enter_data
      procedure, non_overridable :: exit_data => ocean_pressure_force_exit_data
      procedure, non_overridable :: set_bathymetry => ocean_pressure_force_set_bathymetry
      procedure, non_overridable :: bytes => ocean_pressure_force_bytes
   end type ocean_pressure_force_t