ocean_epbl_t Derived Type

   type :: ocean_epbl_t
      logical :: is_init = .false.
         !! True between `init` and `destroy`.

      ! ---- Scheme selection + master switch ----
      logical :: enable = .false.
         !! Master switch.  Default off — all existing namelists and
         !! tests stay bit-identical.  Requires `vmix%use_closure`;
         !! disables the KPP overlay (configure logs the override).
      integer :: mstar_scheme = EPBL_MSTAR_OM4
      integer :: vstar_scheme = EPBL_VSTAR_CUBE_ROOT
      integer :: combine_mode = EPBL_COMBINE_ADD

      ! ---- mstar knobs ----
      real(wp) :: mstar_const = 1.2_wp
         !! Constant-scheme mstar (MOM6 MSTAR).
      real(wp) :: mstar_cap = -1.0_wp
         !! Cap on mstar for the OM4/RH18 schemes; off when < 0.
      real(wp) :: mstar_coef1 = 0.3_wp
         !! OM4 stabilizing-balance coefficient (MOM6 MSTAR2_COEF1).
      real(wp) :: c_ek = 0.085_wp
         !! OM4 Ekman-limit coefficient (MOM6 MSTAR2_COEF2).
      real(wp) :: mstar_conv_adj = 0.0_wp
         !! Reduce mechanical mstar when convection dominates, in
         !! [0,1]; 0 = off (MOM6 MSTAR_CONV_ADJ).
      real(wp) :: rh18_cn1 = 0.275_wp
      real(wp) :: rh18_cn2 = 8.0_wp
      real(wp) :: rh18_cn3 = -5.0_wp
      real(wp) :: rh18_cs1 = 0.2_wp
      real(wp) :: rh18_cs2 = 0.4_wp
         !! RH18 mstar fit coefficients (paper values).

      ! ---- Energetics knobs ----
      real(wp) :: nstar = 0.2_wp
         !! Fraction of convectively released PE that becomes
         !! entrainment-driving TKE (MOM6 NSTAR).
      real(wp) :: tke_decay = 2.5_wp
         !! Ratio of the natural Ekman depth to the mechanical-TKE
         !! decay scale (MOM6 TKE_DECAY).
      real(wp) :: wstar_ustar_coef = 1.0_wp
         !! Weight of the convective reservoir in the velocity scale.
      real(wp) :: vstar_scale_fac = 1.0_wp
         !! Overall vstar multiplier (∝ diffusivity).
      real(wp) :: vstar_surf_fac = 1.2_wp
         !! RH18 vstar scheme: mechanical surface vstar ∝ u* factor.
      real(wp) :: von_karman = 0.41_wp
         !! kappa in Kd = vstar * kappa * mixing_length.
      real(wp) :: ekman_scale_coef = 1.0_wp
         !! Rotational inhibition of the mixing length.
      real(wp) :: min_mix_len = 0.0_wp
         !! Mixing-length floor (m).
      real(wp) :: mixlen_exponent = 2.0_wp
         !! Shape-function exponent (2 = KPP-like; fast path).
      real(wp) :: translay_scale = 0.1_wp
         !! Transition-layer floor of the shape function; must be in
         !! [0,1) when `mld_iteration` is on.

      ! ---- MLD iteration knobs ----
      logical :: mld_iteration = .true.
         !! Iterate the boundary-layer depth to self-consistency with
         !! the mixing-length shape function (MOM6 USE_MLD_ITERATION).
      real(wp) :: mld_tol = 1.0_wp
         !! Convergence tolerance on the MLD root-find (m).
      integer :: mld_max_its = 20
      logical :: mld_bisection = .false.
         !! Plain bisection instead of false-position.
      logical :: mld_use_prev_guess = .false.
         !! Seed the iteration from the previous step's MLD (faster,
         !! ~1-3 iterations; off = always start at half depth).

      ! ---- Environment ----
      real(wp) :: rho0 = 1035.0_wp
         !! Boussinesq reference density (kg/m^3).
      real(wp) :: omega = 7.2921e-5_wp
         !! Earth rotation rate (1/s), for the omega_frac blend.
      real(wp) :: omega_frac = 0.0_wp
         !! Blend |f| with 2*Omega: absf = sqrt((1-of) f^2 + of 4 Om^2).
      real(wp) :: ustar_min = 1.0e-8_wp
         !! Floor on u* (deliberate rdb floor — pure denominator
         !! guard in the TKE decay scale; MOM6 derives ~1e-10).
      real(wp) :: prandtl = 1.0_wp
         !! Kv = prandtl * Kd into the momentum solve.
      logical :: in_eos = .false.
         !! `&ocean_psurf_nml in_eos` (E3): start this scheme's column
         !! pressure stack at `multilayer_state_t%p_top` (the ice-shelf
         !! load + surface pressure, Pa) instead of at 0 Pa.
         !!
         !! The stack has TWO consumers here and the knob moves BOTH,
         !! deliberately: the IN-SITU pressure argument of
         !! `eos_specvol_derivs`, and the PE weight
         !! `dpe_* = dmass*p_mid*dsv_*`.  They are the same pressure --
         !! the weight is the hydrostatic load the layer's centre of mass
         !! has to lift, and under a floating shelf the ice is part of
         !! that load.  Splitting them would put two pressure conventions
         !! in one column, which is the failure the `p_top` seam contract
         !! exists to prevent.
         !!
         !! Host scalar, assigned by `configure_ocean_epbl` BEFORE
         !! `enter_data`, and passed BY VALUE into the column kernel --
         !! never read through the device-mapped handle.
         !!
         !! `.false.` (default) ⇒ the stack starts at 0 Pa, character for
         !! character the pre-E3 arithmetic.  It is NOT enough that
         !! `p_top` be the zero array: under a cavity `p_top` is the ice
         !! load whether or not `in_eos` is set, so this gate is what
         !! keeps an existing cavity + EPBL run bit-identical.
      logical :: tke_diags = .false.
         !! Compute + store the per-column TKE budget terms (W/m^2).
         !! The column ledger closes to round-off — see the design
         !! doc §7; `test_ocean_epbl` asserts it.

      ! ---- Langmuir turbulence (LF17 wind-only path) ----
      ! One scalar La per column per MLD iteration, computed from u*
      ! and the current MLD guess (Li & Fox-Kemper 2017 statistical
      ! waves; no wave model, no Stokes arrays).  Enhancement applies
      ! to mstar only (Reichl & Li 2019).  NOTE: the LF17 branch uses
      ! NEITHER MOM6's MIN_LANGMUIR nor LA_DEPTH_MIN floors — those
      ! belong to the profile-averaging wave paths only.
      logical :: use_lt = .false.
         !! Master Langmuir switch (default off — bit-identity).
      integer :: lt_scheme = EPBL_LT_RESCALE
         !! Enhancement form: rescale (multiplicative) or additive.
      real(wp) :: lt_enhance_coef = 0.447_wp
         !! Enhancement coefficient (MOM6 LT_ENHANCE_COEF).
      real(wp) :: lt_enhance_exp = -1.33_wp
         !! La exponent (MOM6 LT_ENHANCE_EXP).
      real(wp) :: lt_max_enhance = 5.0_wp
         !! Cap on the multiplicative enhancement factor.
      real(wp) :: la_frac_hbl = 0.04_wp
         !! Stokes surface-layer-average depth as a fraction of the
         !! BLD (MOM6 LA_DEPTH_RATIO).
      real(wp) :: lt_lac1 = -0.87_wp
      real(wp) :: lt_lac2 = 0.0_wp
      real(wp) :: lt_lac3 = 0.0_wp
      real(wp) :: lt_lac4 = 0.95_wp
      real(wp) :: lt_lac5 = 0.95_wp
         !! Stability-modified Langmuir number coefficients (MOM6
         !! LT_MOD_LAC1..5: MLD/Ekman, MLD/Obukhov stable, unstable,
         !! Ekman/Obukhov stable, unstable).  All-zero => La_mod = La.

      ! ---- EOS hookup (shared handle from the eos slot) ----
      ! The energy weights use the SAME EOS the dyn-core runs.  We
      ! carry a value copy of the flat-POD `eos_t` (no
      ! allocatable) set once at configure from `ocean_state%eos`;
      ! it maps onto the device with the parent `this` for free and
      ! the per-column `eos_specvol_derivs(this%eos, ...)` call reads
      ! it from registers.  One source of truth — no private scalar
      ! copies that can drift from the dyn-core EOS.
      type(eos_t) :: eos

      ! ---- Persistent fields ----
      real(wp), allocatable :: f_centre(:, :)
         !! |f| at cell centres (1/s); filled by `set_f_centre`.
      real(wp), allocatable :: mld(:, :)
         !! Converged active-mixing-layer depth (m) from the last
         !! call — the next step's first guess (when
         !! `mld_use_prev_guess`) and the primary diagnostic.
      real(wp), allocatable :: b0(:, :)
         !! Surface buoyancy flux (m^2/s^3, > 0 stabilizing) from the last
         !! solve — the same `B0 = g·rho0·(dSV/dT·q_T + dSV/dS·q_S)` the
         !! mstar scaling uses, persisted so the Bodner (2023) MLE
         !! restratification can form its convective velocity scale
         !! `w*^3 = max(0,-b0)·mld`.  Zero until the first EPBL step.
      real(wp), allocatable :: kd_int(:, :, :)
         !! EPBL diapycnal diffusivity at interfaces (m^2/s),
         !! (nx, ny, nz+1); zero at bed (k=1) and surface (k=nz+1).
      real(wp), allocatable :: la(:, :)
         !! Turbulent Langmuir number from the last call (diagnostic;
         !! 0 where `use_lt` is off or the column is dry).

      ! ---- Per-step TKE budget diagnostics (W/m^2) ----
      ! Overwritten each `epbl_compute` (final MLD iteration).  The
      ! ledger: wind + conv + forcing - mixing - mech_decay -
      ! conv_decay = 0 exactly (forcing <= 0, all others >= 0).
      real(wp), allocatable :: tke_wind(:, :)
      real(wp), allocatable :: tke_conv(:, :)
      real(wp), allocatable :: tke_forcing(:, :)
      real(wp), allocatable :: tke_mixing(:, :)
      real(wp), allocatable :: tke_mech_decay(:, :)
      real(wp), allocatable :: tke_conv_decay(:, :)

      ! ---- Iteration-invariant column workspaces ----
      ! Filled once per call by the column prep sweep; read by every
      ! MLD iteration.  All (nx, ny, nz).  Everything else the sweep
      ! needs is carried as scalars (see the design doc D6) — no
      ! per-column work arrays, no NZ_STACK_MAX locals.
      type(scratch_3d_buffer_t) :: t0
         !! Pre-mixing layer temperature (degC).
      type(scratch_3d_buffer_t) :: s0
         !! Pre-mixing layer salinity (PSU).
      type(scratch_3d_buffer_t) :: dpe_t
         !! d(column PE)/d(T_k): rho0*h_k * p_mid * dSV_dT (J/m^2/degC).
      type(scratch_3d_buffer_t) :: dpe_s
         !! d(column PE)/d(S_k) (J/m^2/PSU).
      type(scratch_3d_buffer_t) :: dcolht_t
         !! d(column height)/d(T_k): rho0*h_k * dSV_dT (m/degC) —
         !! steric sensitivity for the gravity-wave radiation term.
      type(scratch_3d_buffer_t) :: dcolht_s
         !! d(column height)/d(S_k) (m/PSU).
      type(scratch_3d_buffer_t) :: ctke_sw
         !! (PR-21) Per-layer TKE cost (J/m^2, <= 0 for solar heating) of
         !! homogenising the penetrating shortwave absorbed IN that layer.
         !! Filled once per call by the prep sweep from the shared two-band
         !! transmission + the in-layer PE-cost shape `Phi(h/zeta)`; the
         !! surface layer's share folds into `ctke_sfc`, the sub-surface
         !! layers drain the sweep's TKE reservoirs.  Zero (never filled)
         !! unless `epbl_sw_ctke .and. sf%has_sw`.  Same shape / lifetime /
         !! device contract as `dcolht_s`.
      logical :: epbl_sw_ctke = .true.
         !! (PR-21) Charge the EPBL TKE ledger for penetrating shortwave
         !! (`&ocean_thermo_nml epbl_sw_ctke`).  Default `.true.`; inert at
         !! `sw_pen_frac = 0` (⇒ `sf%has_sw = .false.`) ⇒ bit-identical.
   contains
      procedure, non_overridable :: init => ocean_epbl_init
      procedure, non_overridable :: destroy => ocean_epbl_destroy
      procedure, non_overridable :: enter_data => ocean_epbl_enter_data
      procedure, non_overridable :: exit_data => ocean_epbl_exit_data
      procedure, non_overridable :: set_f_centre => ocean_epbl_set_f_centre
      procedure, non_overridable :: bytes => ocean_epbl_bytes
   end type ocean_epbl_t