ocean_vdiff_t Derived Type

   type :: ocean_vdiff_t
      logical :: is_init = .false.
         !! True between `init` and `destroy`.
      real(wp) :: K_v_momentum = 0.0_wp
         !! Constant momentum vertical viscosity (m^2/s).  Zero is a
         !! no-op.
      real(wp) :: K_v_tracer = 0.0_wp
         !! Constant tracer vertical diffusivity (m^2/s).  Zero is a
         !! no-op.
      logical :: implicit_stress = .false.
         !! Fold the surface wind stress into the backward-Euler vertical-
         !! friction tridiagonal as a Neumann top-BC right-hand-side term
         !! (`&ocean_vdiff_nml implicit_stress`).  Roundabout is bottom-up, so
         !! the SURFACE is `k = nz`: the kinematic stress `τ/ρ₀` adds to the
         !! `rhs(:, :, nz)` row only (the diagonal is unchanged), filtered
         !! through the implicit operator together with the interior shear.
         !! When `.true.` the explicit surface-stress pre-solve apply in the
         !! dyn run-stage is suppressed (the driver gates it) so the forcing
         !! is not double-counted.  Default `.false.` ⇒ matrix + RHS built
         !! exactly as before ⇒ bit-identical to the explicit path.
      logical :: implicit_drag = .false.
         !! Fold the quadratic / linear bottom drag into the backward-Euler
         !! vertical-friction tridiagonal as a stress bottom-BC diagonal
         !! coupling (`&ocean_vdiff_nml implicit_drag`).  Roundabout bed is
         !! `k = 1`: the drag rate `λ_bot` adds to the `b_diag(:, :, 1)`
         !! diagonal only (a positive add ⇒ unconditionally stable on thin
         !! bottom layers where the explicit `u·(1 − dt·c_d|U|/h)` reverses
         !! sign once `dt·c_d|U|/h > 2`).  `λ_bot` is the SAME Rayleigh rate
         !! (`c_d·|U_bbl|/h_1` quadratic, `r` linear, `|U|` frozen at uⁿ)
         !! the bottom-drag slot already forms; the row is already
         !! normalized by `h_1` so the diagonal increment is `dt·λ_bot`
         !! directly.  Mutually exclusive with `&ocean_bdrag_nml implicit`
         !! (configure fails loud); the explicit drag apply is gated off
         !! when on.  Default `.false.` ⇒ bit-identical.
      logical :: implicit_top_drag = .false.
         !! Fold the ICE-SHELF TOP drag into the `k = nz` DIAGONAL
         !! (`&ocean_vdiff_nml implicit_top_drag`) instead of the explicit
         !! pre-solve add, and MASK the wind-stress RHS off on the faces
         !! the ice covers.  The mirror of `implicit_drag` at the other
         !! end of the column: a drag is a diagonal term and the wind is
         !! an RHS term, so the two share the surface row without either
         !! being approximated, and the cover mask is what makes the
         !! sharing physical (there is no atmosphere under a shelf).
         !! `λ_top` is the SAME Rayleigh rate (`C_d·|U|/h_nz` quadratic,
         !! `r` linear, `|U|` frozen at uⁿ) the top-drag slot forms, and
         !! it arrives already cover- and wet-masked.  Requires
         !! `&ocean_tdrag_nml enable`; mutually exclusive with
         !! `&ocean_tdrag_nml implicit` and with `htbl > 0`; all fail loud
         !! at configure.  The explicit top-drag apply is gated off when
         !! on.  Default `.false.` ⇒ bit-identical.
      logical :: hvel_mom6 = .false.
         !! MOM6 `HARMONIC_VISC = True` parity for the MOMENTUM face thickness
         !! (`hvel`).  Roundabout's `h_u` is the
         !! ARITHMETIC mean unconditionally while MOM6's is the HARMONIC mean
         !! blended back toward arithmetic near the bed when the flow runs
         !! thick->thin -- and MOM6's `h_shear` is then the ARITHMETIC mean of
         !! those hvels.  The two means are effectively SWAPPED relative to
         !! MOM6.  For a grounded sliver `harm(1e-10, 20) ~ 2e-10` against
         !! `arith ~ 10`: eleven orders, and that gap IS the ~1e19 suppression
         !! that renders the spurious sliver PGF harmless in MOM6.
         !!
         !! **Both halves must move together.** LAGRANGIAN_PGF_BUG.md 6.2
         !! records harmonic-`h_u` alone (and harmonic-`h_u` + arithmetic-`dz`)
         !! each making things WORSE -- because neither carried `botfn`.  Plain
         !! harmonic over-suppresses asymmetrically.
      real(wp) :: hbbl_visc = 10.0_wp
         !! Bottom-boundary-layer scale for the `botfn` blend (MOM6 `HBBL`,
         !! 10.0 m in the double_gyre reference) when the BBL glue is off.
         !! Under the per-face glue the blend is normalised by each face's
         !! `bbl_thick` instead (MOM6 `I_Hbbl = 1/bbl_thick`), and this is
         !! only the HBBL fallback for a drag configured with
         !! `&ocean_bdrag_nml hbbl = 0` (MOM6 has ONE `HBBL`).
      logical :: use_harmonic = .false.
         !! MOM6 HARMONIC_VISC analogue.  When `.true.` the implicit
         !! vdiff solver uses the harmonic mean of adjacent layer
         !! thicknesses in the face-thickness denominator instead of
         !! the arithmetic mean.  Better-conditioned at thin /
         !! vanishing layers — `harm = 2·h_a·h_b / (h_a + h_b)` stays
         !! small when one neighbour is small, whereas
         !! `arith = 0.5·(h_a + h_b)` is dominated by the thicker
         !! neighbour and produces stiff coefficients that the Thomas
         !! solve can't handle gracefully.  Default `.false.`
         !! preserves the arithmetic-mean behaviour bit-identically.
      logical :: bbl_glue = .false.
         !! MOM6 `bottomdraglaw` parity for the interface coupling
         !! (`&ocean_vdiff_nml bbl_glue`, PGF_BUG.md §9).  MOM6 holds a
         !! layered basin at rest NOT by computing a clean PGF — its PFu
         !! over grounded/sliver layers is identical to ours — but by
         !! absorbing it viscously every step (`du_dt_visc = −PFu` to
         !! machine precision, measured): `find_coupling_coef` raises the
         !! interface coupling to `a_cpl = kv_bbl/h_shear` within the
         !! bottom boundary layer (botfn weight, `z_i` accumulated from
         !! HARMONIC thicknesses so grounded stacks sit at z≈0), and the
         !! BED coupling is the piston `kv_bbl/(h₁/2)` — unbounded as the
         !! bottom layer thins (measured a_cpl 3e7–1.3e8 m/s on the
         !! rest-state reproducer).  This knob
         !! ports both halves into the momentum tridiagonal:
         !!   * interface kv → `kv + (kv_bbl − KV)·botfn` (`KV` =
         !!     `kv_bbl_bg`), with `h_shear` capped toward `bbl_thick` under
         !!     the same botfn weight;
         !!   * bed row → `dt·kv_bbl/(hf₁·(min(hvel₁/2, bbl_thick)))` — THE
         !!     bed sink whenever the glue is on: it replaces the
         !!     `dt·λ_bot` Rayleigh fold, and the driver skips the explicit
         !!     drag apply (the explicit `du_drag` still feeds the
         !!     barotropic `F_slow`, the split `implicit_drag` takes).
         !! `kv_bbl` and `bbl_thick` are PER FACE: `vdiff_set_viscous_bbl`
         !! (MOM6 `set_viscous_BBL`, quadratic or linear drag law, KW99
         !! rotation/stratification-limited thickness) fills them once per
         !! outer step when `bbl_per_face` (configure: a bottom drag is
         !! configured).  A hand-built slot without `bbl_per_face` gets the
         !! historical constants `bbl_piston·hbbl_visc` / `hbbl_visc`.
         !! Requires `hvel_mom6` (supplies the height-above-bed
         !! bookkeeping), enforced at configure.  With no bottom drag
         !! configured the configure latch turns it off (nothing to glue).
      real(wp) :: bbl_piston = 3.0e-4_wp
         !! BBL drag piston velocity u* (m/s) for `bbl_glue` —
         !! MOM6 `CDRAG·DRAG_BG_VEL` (0.003·0.1 with reference defaults).
         !! `kv_bbl = bbl_piston·hbbl_visc`.
      logical :: hvel_upwind = .true.
         !! Near-bed upwind (arithmetic-donor) blend inside the
         !! `hvel_mom6` face-thickness build.  `.true.` = MOM6 parity
         !! (the historical hvel_mom6 behaviour, bit-identical).
         !! `.false.` = pure harmonic hvel: the blend's u-sign test
         !! flip-flops on roundoff velocities at (near-)rest and
         !! collapses the BBL glue at flipped faces (PGF_BUG.md §9.8) —
         !! the rest-state envelope runs with it off.
      logical :: hvel_harmonic = .false.
         !! Which MOM6 face-thickness branch `hvel_mom6` builds.  `.false.`
         !! (default) is MOM6's DEFAULT, `HARMONIC_VISC = False`
         !! (`vertvisc_coef`): the face thickness is the ARITHMETIC mean,
         !! blended toward the harmonic mean near the bed when the flow runs
         !! from the thin side to the thick one, and the height above the
         !! bed is `max(zh, z_clear)` — `z_clear` the height of the higher of
         !! the two cells' interfaces above the SHALLOWER of the two beds.
         !! So every face layer below the shallower bed of a step sits at
         !! `z_i ≈ 0`, inside the bottom boundary layer, and the BBL glue
         !! couples it.  `.true.` is MOM6 `HARMONIC_VISC = True`: harmonic
         !! face thickness, blended toward arithmetic when the flow runs
         !! thick -> thin, height above bed from harmonic thicknesses alone
         !! (the historical roundabout `hvel_mom6`, gated by `hvel_upwind`).
      logical :: bbl_per_face = .false.
         !! The MOM6 bottom boundary layer (`set_viscous_BBL`) is computed
         !! per face once per outer step by `vdiff_set_viscous_bbl` into
         !! `kv_bbl_u/v` and `bbl_thick_u/v`, which the glue then reads.
         !! Set at configure when `bbl_glue` is on AND a bottom drag is
         !! configured with `&ocean_bdrag_nml hbbl > 0` (MOM6 requires
         !! HBBL).  `.false.` ⇒ the glue uses the historical constants
         !! `bbl_piston·hbbl_visc` / `hbbl_visc`, refilled each call (the
         !! hand-built test path).
      integer :: bbl_form = BBL_FORM_QUADRATIC
         !! Drag law of the BBL (`BBL_FORM_*`).
      real(wp) :: bbl_cd = 0.0_wp
         !! MOM6 `CDRAG`.  Quadratic: `&ocean_bdrag_nml cd`.  Linear: the
         !! equivalent `r·hbbl/bg_vel`, so that `CDRAG·DRAG_BG_VEL = r·hbbl`
         !! — the stress of the HBBL-distributed linear drag, `r·hbbl·u`.
      real(wp) :: bbl_hbbl = 0.0_wp
         !! MOM6 `HBBL` (`&ocean_bdrag_nml hbbl`): the thickness over which
         !! the near-bottom velocity is averaged for the drag law.
      real(wp) :: bbl_bg_vel = 0.0_wp
         !! MOM6 `DRAG_BG_VEL` (`&ocean_bdrag_nml bg_vel`).
      real(wp) :: bbl_thick_min = 0.0_wp
         !! MOM6 `BBL_THICK_MIN` (`&ocean_bdrag_nml bbl_thick_min`; MOM6
         !! default 0).
      logical :: bbl_rino_cap = .false.
         !! MOM6 `RiNo_mix` (= kappa-shear on): the BBL is capped at
         !! `0.5·HBBL`.
      real(wp) :: bbl_rho0 = 1035.0_wp
         !! Boussinesq reference density of the BBL stratification measure.
      real(wp) :: kv_bbl_bg = 1.0e-4_wp
         !! MOM6 `KV`: the background viscosity the BBL viscosity REPLACES
         !! within the boundary layer,
         !! `Kv_tot = Kv_tot + (kv_bbl − KV)·botfn` (`find_coupling_coef`) —
         !! the shear / boundary-layer contributions are kept.  Set from
         !! `&ocean_vmix_nml kv_bg`.
      real(wp), allocatable :: kv_bbl_u(:, :)
         !! BBL viscosity `kv_bbl` (m²/s) at u-faces, `(nx+1, ny)`.
      real(wp), allocatable :: kv_bbl_v(:, :)
         !! BBL viscosity at v-faces, `(nx, ny+1)`.
      real(wp), allocatable :: bbl_thick_u(:, :)
         !! BBL thickness (m) at u-faces, `(nx+1, ny)`.
      real(wp), allocatable :: bbl_thick_v(:, :)
         !! BBL thickness (m) at v-faces, `(nx, ny+1)`.
      real(wp), allocatable :: bbl_conc_t(:, :, :)
         !! Cell temperature CONCENTRATION `(nx, ny, nz)` for the BBL
         !! stratification measure, read through the I1′ column rule
         !! (`rdb_vl_column_conc`: a filler reads its donor).  Workspace of
         !! `vdiff_set_viscous_bbl`; `(1,1,1)` unless `bbl_per_face`.
      real(wp), allocatable :: bbl_conc_s(:, :, :)
         !! Cell salinity concentration, as `bbl_conc_t`.
      logical :: zlevel_faces = .false.
         !! `&vcoord_nml zfixed_closed_faces` — z-level partial steps.
         !! Latched by `configure_ocean_closed_faces`, forwarded as a
         !! plain scalar to `diffuse_velocity_columns_impl`, where it
         !! switches the face column to `min(h_L, h_R)` and CUTS the
         !! tridiagonal coupling across any interface touching a closed
         !! (inert-filler-on-one-side) layer.  See that routine's
         !! `zlevel_faces` docstring for the full argument.  Default
         !! `.false.` ⇒ bit-identical.

      ! ---- Tridiagonal scratch (cell-centred for tracers) ----
      ! Reused across all column solves in one tracer call.  The
      ! Thomas algorithm modifies `c_diag` and `rhs` in place during
      ! the forward sweep, so a separate `c_prime` / `d_prime` pair
      ! isn't needed.
      type(scratch_3d_buffer_t) :: a_diag_t
         !! Sub-diagonal for the tracer solve.  Shape (nx, ny, nz).
      type(scratch_3d_buffer_t) :: b_diag_t
         !! Main diagonal.
      type(scratch_3d_buffer_t) :: c_diag_t
         !! Super-diagonal (overwritten by `c_prime`).
      type(scratch_3d_buffer_t) :: rhs_t
         !! Right-hand side / solution.

      ! ---- Momentum tridiagonal scratch (face-located) ----
      ! u-face shape is (nx+1, ny, nz); v-face shape is (nx, ny+1, nz).
      ! Separate buffers keep the descriptor extents exact so the
      ! `do concurrent` indexing doesn't drift past the allocated
      ! footprint.
      type(scratch_3d_buffer_t) :: a_diag_u, b_diag_u, c_diag_u, rhs_u
         !! East-face (u) tridiagonal scratch.
      type(scratch_3d_buffer_t) :: a_diag_v, b_diag_v, c_diag_v, rhs_v
         !! North-face (v) tridiagonal scratch.

      ! ---- Cell-centred diffusivity workspace ----
      ! Shape (nx, ny, nz+1).  Filled with `K_v_*` (broadcast from
      ! the scalar) when the apply routines are called without a
      ! `kv_source` argument; left alone when the caller passes a
      ! closure-produced 3D field directly.  k=1 is the bed
      ! interface (forced to zero — closed bed BC), k=nz+1 is the
      ! free surface (also zero — closed top).
      type(scratch_3d_buffer_t) :: kv_scalar_buf
         !! Diffusivity workspace used when no 3D source is provided.
   contains
      procedure, non_overridable :: init => ocean_vdiff_init
      procedure, non_overridable :: destroy => ocean_vdiff_destroy
      procedure, non_overridable :: enter_data => ocean_vdiff_enter_data
      procedure, non_overridable :: exit_data => ocean_vdiff_exit_data
      procedure, non_overridable :: bytes => ocean_vdiff_bytes
   end type ocean_vdiff_t