ocean_horizontal_viscosity_t Derived Type

   type :: ocean_horizontal_viscosity_t
      logical :: is_init = .false.
         !! True between `init` and `destroy`.  Tracks GPU device
         !! attachment too — prefer this to `allocated(...)`.
      real(wp) :: nu_h = 0.0_wp
         !! Constant Laplacian horizontal viscosity (m^2/s).  Phase
         !! Tier-1 default is zero (kernel becomes a no-op); set
         !! positive to enable damping.  Numerical-stability cap:
         !! `nu_h * dt * (1/dx^2 + 1/dy^2) <= 0.5`.
      logical :: stress_tensor = .false.
         !! When true, the compute step uses the MOM6-faithful
         !! thickness-weighted stress-divergence operator
         !! `diffu = (1/(h_u+h_neglect))·∇·(h·A·∇u)` with a per-cell
         !! CFL viscosity limiter and `wet_*` coast-masking instead of
         !! the velocity Laplacian `A·∇²u`.  Default false ⇒
         !! bit-identical to the historical kernel.  See the module
         !! header for the operator form.
      logical :: no_slip = .false.
         !! Coastal lateral BC selector (shared with the lateral-mix /
         !! Coriolis kernels).  `.false.` (free-slip) masks the corner
         !! shear stress by `wet_q`; `.true.` (no-slip) by `2 - wet_q` —
         !! on the `stress_tensor` path AND on the velocity-Laplacian
         !! paths (scalar `nu_h` and the flow-aware face closures), whose
         !! corner (shear) fluxes carry the same factor.  The biharmonic
         !! add-on is always free-slip (`wet_q`; MOM6 refuses NOSLIP with
         !! BIHARMONIC).  All-wet ⇒ factor ≡ 1 ⇒ bit-identical.
      real(wp) :: bound_coef = 0.8_wp
         !! CFL safety coefficient for the per-cell viscosity limiter
         !! (MOM6 `HORVISC_BOUND_COEF`).  Consulted on the
         !! `stress_tensor` path and, when `bound_kh` is set, on the
         !! velocity-Laplacian paths too.
      logical :: bound_kh = .false.
         !! MOM6 `BOUND_KH` analogue for the velocity-Laplacian paths
         !! (scalar `nu_h` and the flow-aware per-face closure).  When
         !! `.true.`, the per-face harmonic viscosity is clamped to
         !! `bound_coef·0.125/(dt·(idx²+idy²))` — ~0.25× the forward-
         !! Euler stability limit, matching MOM6's `Kh_Max_xx` margin.
         !! Load-bearing beyond simple FE stability: the barotropic
         !! mode receives the depth-mean viscous force FROZEN over the
         !! outer step (via `F_bt`), and for grid-scale gravity modes
         !! with `ω·dt ≳ 1` an unbounded `λ·dt = ν·k²·dt ≳ 0.6` frozen
         !! force is applied with reversed phase — ANTI-damping — which
         !! exponentially pumps rim-trapped barotropic modes through
         !! the Δu corrector (the 600² Lagrangian double-gyre h-guard
         !! blow-up; e-fold ~10 outer steps).  The clamp keeps
         !! `λ·dt ≤ 0.5·bound_coef` at every face so the corrector
         !! loop stays damped at any resolution.  Default `.false.`
         !! ⇒ bit-identical.

      ! ---- Anisotropic viscosity (Smith & McWilliams 2003) ----
      real(wp) :: kh_aniso = 0.0_wp
         !! Anisotropic Laplacian viscosity magnitude (m²/s).  When
         !! positive, a two-coefficient direction tensor splits the
         !! harmonic viscosity into along- and cross-direction parts.
         !! Only consulted on the `stress_tensor` path.  Default 0 ⇒
         !! isotropic, bit-identical.  MOM6 `KH_ANISO` analogue
         !! (Smith & McWilliams 2003, "Anisotropic horizontal
         !! viscosity for ocean models", Ocean Modelling 5(2), §2).
         !! Composes with the biharmonic add-on (`nu_4` / `smag_ah` /
         !! `leith_biharm`) — the two are independent linear operators
         !! that sum, matching MOM6's composition (`kh_aniso` inside
         !! the harmonic block, biharmonic added to the same tensor).
      real(wp) :: aniso_n1n2 = 0.0_wp
         !! Precomputed direction-tensor factor `2·n1·n2/(n1²+n2²)`
         !! (MOM6 `n1n2`).  Set by `ocean_hvisc_set_aniso_direction`
         !! from the `(n1,n2)` direction vector.  Default `(1,0)` =
         !! grid-i ⇒ `n1n2 = 0` (cross terms vanish; the operator just
         !! adds `kh_aniso` to the tension coefficient).
      real(wp) :: aniso_n1n1_m_n2n2 = 1.0_wp
         !! Precomputed direction-tensor factor `(n1²−n2²)/(n1²+n2²)`
         !! (MOM6 `n1n1_m_n2n2`).  Default `(1,0)` ⇒ 1.
      real(wp) :: nu_4 = 0.0_wp
         !! Constant biharmonic horizontal viscosity (m^4/s).  Scale-
         !! selective damping for stratified closed-basin runs: damps
         !! proportional to ν₄·k⁴, so it kills grid-scale baroclinic
         !! noise without bleeding into resolved scales the way a
         !! large Laplacian ν_h would.  Required to keep stratified
         !! Tasman-class runs bounded past day ~10 (Laplacian-only
         !! configurations grow exponentially via parametric
         !! amplification of roundoff seeds).  Numerical-stability
         !! cap: `ν₄ · dt · (1/dx² + 1/dy²)² <= 1/16`.  Typical 2 km
         !! resolution value ~1e10 m^4/s.  Composes with **all** of
         !! the harmonic dispatch arms, including `stress_tensor` — it
         !! is no longer silently disabled by it.

      type(scratch_3d_buffer_t) :: du_visc
         !! Per-step viscous tendency at east faces, shape
         !! (nx+1, ny, nz).  Filled by the compute step, consumed
         !! by the apply step.
      type(scratch_3d_buffer_t) :: dv_visc
         !! Per-step viscous tendency at north faces, shape
         !! (nx, ny+1, nz).
      type(scratch_3d_buffer_t) :: lap_u
         !! First-pass Laplacian buffer used by the biharmonic path —
         !! holds `∇²u_face` so the second Laplacian pass (`∇²(∇²u)`)
         !! can read it.  Same shape as `du_visc`.  Allocated
         !! unconditionally; sits inert when `nu_4 = 0`.
      type(scratch_3d_buffer_t) :: lap_v
         !! v-face counterpart of `lap_u`.

      ! ---- MOM6 stress-divergence scratch (stress_tensor path only) ----
      type(scratch_3d_buffer_t) :: str_xx
         !! Thickness-weighted tension stress `A_T·(du/dx−dv/dy)·h_T`
         !! at T-cell centres, shape `(nx, ny, nz)`.
      type(scratch_3d_buffer_t) :: str_xy
         !! Thickness-weighted shear stress
         !! `A_q·(dv/dx+du/dy)·h_q·slip` at Bu corners, shape
         !! `(nx+1, ny+1, nz)`.
      type(scratch_3d_buffer_t) :: ah_t
         !! Harmonic viscosity averaged onto T-cell centres (m²/s),
         !! shape `(nx, ny, nz)`.  CFL-clamped per cell.
      type(scratch_3d_buffer_t) :: ah_q
         !! Harmonic viscosity averaged onto Bu corners (m²/s),
         !! shape `(nx+1, ny+1, nz)`.  CFL-clamped per corner.

      ! ---- MEKE frictional-source coupling (capability [5] seam) ----
      logical :: compute_ke_diss = .false.
         !! When set (by configure when MEKE's frictional source is on),
         !! the apply step fills `ke_diss` with the lateral-viscosity KE
         !! dissipation rate.  Default off ⇒ no extra work, bit-identical.
      real(wp), allocatable :: ke_diss(:, :)
         !! Depth-integrated KE dissipation rate by the lateral viscosity,
         !! `Σ_k ρ_k h_k (u·du_visc + v·dv_visc)` (kg/s³, ≤0), at T-cell
         !! centres `(nx, ny)`.  MEKE consumes it as `-frcoeff·i_mass·ke_diss`
         !! (the mean→eddy frictional source).  Holds the most recent stage's
         !! rate (the quantity MEKE needs is a rate, so no accumulation).
   contains
      procedure, non_overridable :: init => ocean_hvisc_init
      procedure, non_overridable :: destroy => ocean_hvisc_destroy
      procedure, non_overridable :: enter_data => ocean_hvisc_enter_data
      procedure, non_overridable :: exit_data => ocean_hvisc_exit_data
      procedure, non_overridable :: bytes => ocean_horizontal_viscosity_bytes
   end type ocean_horizontal_viscosity_t