rdb_constants.F90 Source File

!! Physical constants and working precision for the tidal solver
module rdb_constants
   !! Working precision kind parameter and physical constants.
   use pic_types, only: dp, sp
   implicit none
   private

   public :: wp
   public :: GRAVITY, RHO_WATER, OMEGA_EARTH
   public :: CP_WATER
   public :: LATENT_HEAT_VAPORIZATION
   public :: DRY_TOLERANCE
   public :: NH_MIN_DEPTH
   public :: THIN_LAYER_THRESHOLD
   public :: H_VANISHED, H_DIV_EPS
   public :: LAND_DEPTH_THRESHOLD
   public :: FROUDE_CAP
   public :: PI, TWO_PI, DEG2RAD, RAD2DEG
   public :: NZ_STACK_MAX
   public :: nz_stack_required, nz_stack_is_sufficient
   public :: VCOORD_LAGRANGIAN, VCOORD_EULERIAN_Z, &
             VCOORD_SIGMA, VCOORD_ZSIGMA, VCOORD_ZSTAR, &
             VCOORD_ZSTAR_FULL, VCOORD_ZSTAR_SIGMA, VCOORD_Z_FIXED, &
             VCOORD_RHO, VCOORD_HYCOM
   public :: REMAP_PCM, REMAP_PLM, REMAP_PPM, REMAP_PPM_H4, REMAP_PQM
   public :: KEPS_STAB_CONSTANT, KEPS_STAB_GALPERIN, KEPS_STAB_CANUTO
   public :: FRHAT_ARITHMETIC, FRHAT_HYBRID, parse_frhat_scheme

#ifdef RDB_DOUBLE_PRECISION
   integer, parameter :: wp = dp
      !! Working precision kind parameter (double)
#else
   integer, parameter :: wp = sp
      !! Working precision kind parameter (single)
#endif

   real(wp), parameter :: GRAVITY = 9.80665_wp
      !! Gravitational acceleration (m/s^2)
   real(wp), parameter :: RHO_WATER = 1035.0_wp
      !! Seawater density (kg/m^3)
   real(wp), parameter :: OMEGA_EARTH = 7.2921e-5_wp
      !! Earth angular velocity (rad/s)
   real(wp), parameter :: CP_WATER = 4000.0_wp
      !! Specific heat capacity of seawater (J/(kg·K)).  Converts surface
      !! heat flux Q_net [W/m²] to kinematic heat flux Q_net/(rho_0·cp) [m·K/s].
   real(wp), parameter :: LATENT_HEAT_VAPORIZATION = 2.5e6_wp
      !! Latent heat of vaporization of water (J/kg).  A representative
      !! near-surface value (2.501e6 at 0 °C, weakly T-dependent).  Used by
      !! the evaporative latent-heat surface cooling: a net evaporation of
      !! E [m/s] withdraws a heat flux ρ_water·E·L_v [W/m²] from the top of
      !! the water column as the evaporated mass changes phase.
   real(wp), parameter :: DRY_TOLERANCE = 1.0e-6_wp
      !! Minimum depth to consider a cell wet (m)
   real(wp), parameter :: NH_MIN_DEPTH = 1.0e-1_wp
      !! Minimum depth for non-hydrostatic pressure operations (m).
      !! Cells shallower than this are treated as hydrostatic to avoid
      !! 1/h singularities in the NH Poisson operator.
   real(wp), parameter :: THIN_LAYER_THRESHOLD = 1.0e-3_wp
      !! Below this depth, momentum is zeroed to prevent
      !! huge velocities at wet/dry fronts (m)

   ! ---- Thin-layer constants of record (ocean-path semantic split, D4 taxonomy) ----
   real(wp), parameter :: H_VANISHED = 1.5e-4_wp
      !! Dynamic-vanish threshold (m).  A layer thinner than this is
      !! dynamically VANISHED — skip / merge into a neighbour, do NOT clamp.
      !! When dividing by `h_layer` for a PHYSICAL result, gate on this.
      !! Tests of it are STRICT `>`: a layer sitting exactly ON the marker
      !! reads as vanished, which is what the geometric vcoord families
      !! rely on (`rdb_vcoord :: vcoord_h_min_role`).
      !!
      !! **Do not hand-roll the test.** Tracers are stored as CONTENT
      !! (`hTr = h·c`) and the tree holds one invariant about vanished
      !! layers —
      !!
      !!     I1′:  `h <= H_VANISHED  ⇒  hTr = h·c_live`, every registered
      !!           tracer, `c_live` the concentration of the filler's donor
      !!           (nearest live layer above; topmost live layer for the
      !!           fillers above it; `hTr = 0` with no live layer)
      !!
      !! — with ONE definition of the predicate, of "the concentration of
      !! layer k", and of the column merge that restores it:
      !! `src/shared_module_utilities/rdb_vanished_layer.inc`
      !! (`rdb_vl_is_live`, `rdb_vl_conc`, `rdb_vl_column_conc`,
      !! `rdb_vl_holds_live_conc`, `rdb_vl_merge_content`), ONE
      !! enforcement point (`multilayer_state_t%enforce_vanished_content`,
      !! at the tail of the outer step) and a fail-loud tripwire
      !! (`&vcoord_nml check_vanished_content`).  The contract, and the
      !! table of consumers that deliberately substitute something ELSE
      !! (the EOS's reference T/S, the sponge's nearest massive layer, the
      !! diagnostics' NaN, the melt sampler's skip), is written down once
      !! in `src/core/ocean/README.md`, "The vanished-layer content rule".
      !! A new raw `hTr/h` divide or a new comparison against this constant
      !! outside the sanctioned modules is flagged by the `vanished-layer`
      !! pre-commit hook.
   real(wp), parameter :: H_DIV_EPS = 1.0e-20_wp
      !! Pure division-safety epsilon (m), far below any physical thickness:
      !! adding it (`1/(h+H_DIV_EPS)`) only prevents a 1/0, never changes a
      !! well-posed result.  Use where a divisor is ALREADY guaranteed
      !! positive (MOM6 `H_subroundoff` role).
   real(wp), parameter :: LAND_DEPTH_THRESHOLD = 2.0_wp
      !! Cell bathymetry depth (m, positive-down) below which an ocean cell
      !! is treated as LAND for surface-forcing masking.  Consumed by
      !! `ocean_state_seed_from_cfg` to populate `multilayer%wet_mask`.
   real(wp), parameter :: PI = 4.0_wp*atan(1.0_wp)
      !! Pi
   real(wp), parameter :: TWO_PI = 8.0_wp*atan(1.0_wp)
      !! 2*Pi
   real(wp), parameter :: DEG2RAD = PI/180.0_wp
      !! Degrees → radians conversion factor.
   real(wp), parameter :: RAD2DEG = 180.0_wp/PI
      !! Radians → degrees conversion factor.
   real(wp), parameter :: FROUDE_CAP = 10.0_wp
      !! Maximum Froude number allowed at wet/dry fronts.
      !! If |u|/sqrt(g*h) exceeds this, momentum is rescaled.
      !! Physical flows rarely exceed Fr=3; Fr=10 is very permissive.

#ifndef RDB_NZ_STACK_MAX
#define RDB_NZ_STACK_MAX 128
#endif
   integer, parameter :: NZ_STACK_MAX = RDB_NZ_STACK_MAX
      !! Maximum vertical layers for per-column stack workspace inside
      !! `do concurrent`.  Stack arrays in the BPG, remap, kappa-shear,
      !! Redi and diag-remap kernels are dimensioned by this constant so
      !! the compiler can emit fixed-size thread-local storage on GPU.
      !!
      !! **The requirement is `NZ_STACK_MAX >= nz + 1`** — see
      !! `nz_stack_required` below for the derivation and the survey of
      !! every consumer.  Enforced fail-loud by `validate_config`.
      !!
      !! History: 100 -> 256 (2026-06-19), on the belief that the Redi
      !! neutral-diffusion sweep needed `2*nz+2` *of this constant*.  It
      !! does not: Redi's `nsurf` locals are declared `2*NZ_STACK_MAX+2`
      !! (`rdb_ocean_redi.F90` `PoLc`/`PoRc`/`KoLc`/`KoRc`/`hEc`), so they
      !! scale with the constant and only need `NZ_STACK_MAX >= nz`.
      !! Lowered 256 -> 128 (2026-08-05) after that measurement: 256 was
      !! ~2x larger than any consumer needed and cost ~4 GB of CUDA
      !! per-thread local memory, which the driver reserves at each
      !! kernel's first launch and never returns.  128 covers the largest
      !! shipped case (nz=90, `acc_channel_kitchensink_xl.nml`) and the
      !! largest test (nz=100) with headroom.
      !!
      !! Heap-allocated workspace arrays (PPM snapshots, remap snapshots)
      !! have no such limit — they are allocated at runtime.

   ! ---- Vertical coordinate type constants ----

   integer, parameter :: VCOORD_LAGRANGIAN = -1
      !! Pure Lagrangian / isopycnal — `h_layer` evolves freely and is
      !! never remapped.  `compute_target_h` returns the current `h`
      !! unchanged and the ALE remap step is a no-op.  Distinct from
      !! `VCOORD_EULERIAN_Z` (which targets the static `H · dsig`
      !! decomposition): LAGRANGIAN tracks material surfaces, so for
      !! adiabatic flow it reproduces MOM6's `COORD_CONFIG="gprime"`
      !! NK=2 isopycnal layering.  Parsed from "lagrangian" / "isopycnal".
   integer, parameter :: VCOORD_EULERIAN_Z = 0
      !! Ignore SSH when forming `target_h`; layer thicknesses follow
      !! the static `H · dsig(k)` decomposition.  Used by the ocean
      !! path as the "do nothing, leave the IC layers alone" default
      !! before the user picks a remap-aware coord_type.  Was originally
      !! declared in `rdb_ocean_vcoord` as an orphan; moved here so the
      !! VCOORD enum lives in one place.
   integer, parameter :: VCOORD_SIGMA = 1
      !! Pure terrain-following sigma coordinate (default)
   integer, parameter :: VCOORD_ZSIGMA = 2
      !! z-sigma hybrid: z-levels in deep/steep, sigma in shallow
   integer, parameter :: VCOORD_ZSTAR = 4
      !! z-star "lite" (MOM6-like quasi-horizontal): single global
      !! `z_ref` pattern stretched proportionally per column by H/z_ref(nz)
   integer, parameter :: VCOORD_ZSTAR_FULL = 5
      !! z-star "full" MOM6: per-column `z_ref` anchored to local
      !! bathymetry with vanishing layers below the bed.  Surface
      !! layers stay at a fixed physical thickness independent of H
   integer, parameter :: VCOORD_ZSTAR_SIGMA = 6
      !! z*/sigma hybrid: smoothstep blend from pure sigma in shallow
      !! water (H ≤ depth_transition) to z*-lite in deep water
      !! (H ≥ depth_transition + blend_width).  Coastal wet/dry
      !! friendliness in shallow zones, SSH-tracking in deep zones
   integer, parameter :: VCOORD_Z_FIXED = 7
      !! Fixed-z interfaces with vanishing layers in shallow water.
      !! ALE remap pulls `h_layer(k) → h_target(k)` where the targets
      !! come from absolute-depth interfaces `z_target = (0, h_nominal,
      !! 2·h_nominal, ..., h_ref)` with `h_nominal = h_ref/nz_ml` — or,
      !! under `&vcoord_nml z_fixed_profile = "list" | "tanh"`, from a
      !! stretched per-layer table (`ocean_vcoord_t%z_fixed_zi`).
      !! In deep cells the interfaces stay locked at those depths; in
      !! shallow cells the bed-side layers vanish to a minimum
      !! thickness `h_min` while the surface layer absorbs the residual.
      !! Mirrors MOM6's `COORD_CONFIG = "gprime"` layer structure and
      !! the existing `seed_h_layer_z_fixed_impl` IC algorithm — used
      !! as a per-step ALE target rather than an init-only seed so the
      !! interface stays anchored as the simulation evolves.
   integer, parameter :: VCOORD_RHO = 8
      !! Isopycnal coordinate (P2).  Layer interfaces are placed on
      !! prescribed potential-density surfaces `rho_target(0:nz)`
      !! (monotone-increasing, referenced to `rho_ref_pressure`).  The
      !! ALE regrid inverts a PPM reconstruction of the column density
      !! profile for the depth where ρ equals each interior target, then
      !! remaps T/S/tracers/velocities onto the new grid.  Bottom-up
      !! convention: `rho_target(0)` (lightest) anchors the surface
      !! interface (k=nz); `rho_target(nz)` (densest) the bed (k=1).
      !! Validation-grade alone (weakly-stratified columns collapse);
      !! HYCOM hybrid is the production follow-on (P3).  See
      !! `rdb_ocean_vcoord :: ocean_vcoord_compute_target_h_rho_impl`.
   integer, parameter :: VCOORD_HYCOM = 9
      !! Hybrid z*/isopycnal coordinate (P3, Bleck 2002 / MOM6
      !! `coord_hycom`).  Runs the exact `VCOORD_RHO` density-space
      !! inversion (same `ocean_vcoord_compute_target_h_rho` kernel,
      !! `hybrid=.true.`) but adds two HYCOM deltas around it: a
      !! bottom-up density monotonize before the inversion, and a z*
      !! nominal-floor sweep after it that pushes too-shallow interfaces
      !! down to a fixed-resolution surface band (`dsig·(H+η)`).  Net
      !! effect: a fixed-res near-surface z* band (no surface collapse)
      !! with an isopycnal interior — the production GVC coordinate.
      !! Reuses `rho_target` / `rho_ref_pressure`; no new knobs.

   ! ---- Vertical remapping method constants ----

   integer, parameter :: REMAP_PCM = 1
      !! Piecewise constant (0th order, donor cell)
   integer, parameter :: REMAP_PLM = 2
      !! Piecewise linear with minmod limiter (1st order)
   integer, parameter :: REMAP_PPM = 3
      !! Piecewise parabolic (2nd order, Colella & Woodward)
   integer, parameter :: REMAP_PPM_H4 = 4
      !! Piecewise parabolic with non-uniform 4th-order (H4) edge values —
      !! thickness-weighted edge estimate (White & Adcroft 2008) feeding the
      !! same Colella & Woodward parabola + monotonicity limiter as REMAP_PPM
   integer, parameter :: REMAP_PQM = 5
      !! Piecewise quartic (4th-order, PQM_IH4IH3 — White & Adcroft 2008):
      !! implicit-h4 edge values + implicit-h3 edge slopes + a quartic per-cell
      !! reconstruction with the White & Adcroft monotonicity limiter and a
      !! conservative quartic overlap integral.  Requires nz >= 5; falls back to
      !! REMAP_PPM for thinner columns.  Boundary cells are PCM.

   ! ---- k-epsilon stability-function scheme constants ----

   integer, parameter :: KEPS_STAB_CONSTANT = 0
      !! Constant c_mu (=0.09) / Pr_t (=1) -- standard high-Re k-epsilon
      !! (default; bit-identical to the original closure).
   integer, parameter :: KEPS_STAB_GALPERIN = 1
      !! Galperin et al. (1988) quasi-equilibrium stratification-aware
      !! stability functions (Mellor-Yamada level-2.5 algebraic form).
   integer, parameter :: KEPS_STAB_CANUTO = 2
      !! Canuto et al. (2001) Model A second-moment-closure stability functions
      !! (rational alpha_M / alpha_N form per Umlauf & Burchard 2003).

   ! ---- Barotropic face-layer-weight ("frhat") scheme constants ----

   integer, parameter :: FRHAT_ARITHMETIC = 0
      !! `&ocean_bt_nml frhat_scheme = "arithmetic"` (DEFAULT pre-port
      !! behaviour).  Every barotropic depth mean
      !! (`rdb_barotropic_coupling::derive_bt_from_layers`,
      !! `face_depth_mean_u/v`, `face_depth_mean_rem_u/v`,
      !! `apply_bt_correction`'s open/visc_rem folds,
      !! `compute_bt_rem_from_visc_rem`'s `av_rem`) weights each layer's
      !! face thickness by the plain two-abutting-cell arithmetic mean,
      !! `h_face(k) = 0.5*(h_L(k) + h_R(k))`.  Bit-identical to the
      !! pre-frhat-port tree.
   integer, parameter :: FRHAT_HYBRID = 1
      !! `&ocean_bt_nml frhat_scheme = "hybrid"` — MOM6's
      !! `HVEL_SCHEME = HYBRID` face-thickness closure (`btcalc`'s
      !! default scheme there).
      !! Above the shallower of the two abutting columns' bed depths the
      !! face thickness is the arithmetic mean (as `FRHAT_ARITHMETIC`);
      !! below it, the harmonic mean `h_L*h_R/(h_L+h_R)` — which
      !! vanishes with the thinner side instead of being dragged up by
      !! the thicker one — with a linear blend across the transition
      !! layer.  Suppresses a thin partial-bed layer's weight next to a
      !! thick abutting layer, the mechanism `docs/
      !! visc_rem_bt_rem_plan.md` Section 7 identifies as still missing
      !! from every roundabout barotropic depth mean.  See
      !! `rdb_barotropic_coupling::frhat_h_face_step`
      !! (`src/shared_module_utilities/rdb_frhat_face.inc`) for the
      !! bottom-up port of MOM6's `e_u`/`D_shallow_u` recursion.

contains

   pure integer function parse_frhat_scheme(str, default_code) result(scheme)
      !! Convert an `&ocean_bt_nml frhat_scheme` namelist string to a
      !! `FRHAT_*` constant.  Unrecognised ⇒ `default_code` if given,
      !! else `FRHAT_ARITHMETIC`.
      character(len=*), intent(in) :: str
      integer, intent(in), optional :: default_code
      integer :: fallback
      fallback = FRHAT_ARITHMETIC
      if (present(default_code)) fallback = default_code
      select case (trim(adjustl(str)))
      case ("arithmetic", "ARITHMETIC")
         scheme = FRHAT_ARITHMETIC
      case ("hybrid", "HYBRID")
         scheme = FRHAT_HYBRID
      case default
         scheme = fallback
      end select
   end function parse_frhat_scheme

   pure function nz_stack_required(nz) result(req)
      !! Smallest `NZ_STACK_MAX` that safely covers a run of `nz` layers.
      !!
      !! `nz + 1`.  Derived from a survey of every `NZ_STACK_MAX`-dimensioned
      !! array in the tree; the binding consumers are the ones that index a
      !! plainly-`NZ_STACK_MAX`-declared array to `nz+1` rather than the
      !! `NZ_STACK_MAX+1` an interface array would use:
      !!
      !! Interface-indexed column workspaces are the ones to watch: a
      !! tridiagonal over layer INTERFACES runs to `nz + 1`, so a
      !! plainly-`NZ_STACK_MAX`-declared scratch array must still admit
      !! that extra slot.  The remap kernel's own `z_old(0:NZ_STACK_MAX)` /
      !! `q_L`/`q_R`/`q6` are sized on the same rule.
      !!   * `rdb_ocean_diag_fills` `z_iface` — needs `n_bin + 2`, which the
      !!     `n_bin <= NZ_STACK_MAX - 1` output-level guard bounds by `nz + 1`.
      !!
      !! NOT `2*nz + 2`: the Redi neutral-surface locals that motivated the
      !! 256 bump are declared `2*NZ_STACK_MAX + 2` and so scale with the
      !! constant — they need only `NZ_STACK_MAX >= nz`.
      integer, intent(in) :: nz
      integer :: req

      req = nz + 1
   end function nz_stack_required

   pure function nz_stack_is_sufficient(nz) result(ok)
      !! `.true.` when the COMPILED `NZ_STACK_MAX` covers `nz` layers.
      !! Callers that get `.false.` must refuse the run — the overflow is
      !! a silent thread-local-storage overrun, not a crash.
      integer, intent(in) :: nz
      logical :: ok

      ok = nz_stack_required(nz) <= NZ_STACK_MAX
   end function nz_stack_is_sufficient

end module rdb_constants