tidal_mixing_column_kernel Subroutine

   pure subroutine tidal_mixing_column_kernel(grid, this, ms, hT, hS)
      !! Per-column St-Laurent flux-bookkeeping sweep.  One
      !! `do concurrent (j, i)` over owned cells; each column is a serial
      !! upward sweep k=1(bed) -> nz(surface) over fixed-size `local()`
      !! arrays (register-resident).  Bottom-up global indexing
      !! throughout — no surface-down flip (the decay anchors at the bed).
      !!
      !! "The bed" is `kbed = ms%k_bot(i,j)`, the first LIVE layer counting
      !! up (`1` off `z_fixed`, ⇒ bit-identical there).  Under `z_fixed`
      !! the layers below it are inert fillers carrying the donor's T/S, so
      !! anchoring at `k = 1` read `N_bot = 0` across two fillers (killing
      !! the `e_compute` energy input on every column shallower than the
      !! nominal stack), exempted a FILLER from deposition instead of the
      !! bed layer, and handed the fillers a TKE share.  The sweep, the
      !! `N_bot` sample, the bed-layer exclusion and the bed end-cap all
      !! start at `kbed`; layers below it get no `Kd`.
      !!
      !! Pipeline per column:
      !!   gather: h (floored), per-layer T,S = hTr/h.
      !!   N^2(k): layer-centred buoyancy gradient (St-Laurent #4 —
      !!           N2_lay, not interface N^2), clamped >= 0.
      !!   E:      prescribed `e_in(i,j)` (v1) or state-dependent
      !!           TKE_coef*N_bot (v1.1, `e_compute`), masked where
      !!           H < min_zbot.
      !!   sweep:  Inv_int normalization, TKE flux bookkeeping, per-layer
      !!           Kd_add capped at kd_max, deposited 50/50 at the two
      !!           bounding interfaces.
      type(hgrid_t), intent(in) :: grid
      type(ocean_tidal_mixing_t), intent(inout) :: this
      type(multilayer_state_t), intent(in) :: ms
      ! assumed-shape-ok: tracer registry outer-shim — caller host-dereferences
      ! ms%tracers(idx)%hTr before passing; size varies per tracer slot
      ! (see CLAUDE.md "outer-shim + flat-impl" pattern); thermo cadence.
      real(wp), intent(in) :: hT(:, :, :)
         !! Temperature tracer hTr (degC*m), host-dereferenced.
      real(wp), intent(in) :: hS(:, :, :)  ! assumed-shape-ok: tracer registry outer-shim; thermo cadence
         !! Salinity tracer hTr (PSU*m), host-dereferenced.

      integer :: i, j, k, nx, ny, nz, kbed
      real(wp) :: e_col, h_tot, n_bot, tke_coef, h2c, e_cap
      real(wp) :: inv_int, hz, tke_bot, tke_rem, z_top, frac_top, tke_lay
      real(wp) :: dz_eff, denom, kd_lay
      real(wp) :: h_col(NZL), t_col(NZL), s_col(NZL)
      real(wp) :: n2_col(NZL), kd_lay_arr(NZL)
      real(wp) :: dbuoy_t, dbuoy_s, p_int
      real(wp) :: dsv_dt_k, dsv_ds_k, t_int, s_int
      real(wp) :: gamma_l, mu_l, zeta_l, kd_max_l, omega2_l, rho0_l, izeta

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml

      gamma_l = this%gamma
      mu_l = this%mu
      zeta_l = this%zeta
      kd_max_l = this%kd_max
      omega2_l = this%omega2
      rho0_l = this%rho0
      izeta = 1.0_wp/zeta_l

      do concurrent(j=1:ny, i=1:nx) &
         local(k, kbed, e_col, h_tot, n_bot, tke_coef, h2c, e_cap, &
               inv_int, hz, tke_bot, tke_rem, z_top, frac_top, tke_lay, &
               dz_eff, denom, kd_lay, h_col, t_col, s_col, n2_col, kd_lay_arr, &
               dbuoy_t, dbuoy_s, p_int, dsv_dt_k, dsv_ds_k, t_int, s_int)

         ! ---- Default output: zero contribution (overwritten if wet) ----
         do k = 1, nz + 1
            this%kd_int(i, j, k) = 0.0_wp
         end do

         if (ms%wet_mask(i, j) > 0.0_wp) then
            kbed = ms%k_bot(i, j)
            ! ---- Gather (bottom-up; floor thickness, back out T,S) ----
            ! `h_tot` is the LIVE column (the fillers below `kbed` are not
            ! water the internal tide can mix).
            h_tot = 0.0_wp
            do k = 1, nz
               dz_eff = max(ms%h_layer(i, j, k), H_VANISHED)
               h_col(k) = dz_eff
               denom = 1.0_wp/max(ms%h_layer(i, j, k), H_DIV_EPS)
               t_col(k) = hT(i, j, k)*denom
               s_col(k) = hS(i, j, k)*denom
               if (k >= kbed) h_tot = h_tot + ms%h_layer(i, j, k)
            end do

            ! ---- Layer-centred N^2 (St-Laurent #4: N2_lay) ----
            ! Buoyancy gradient between the centres of layer k and the
            ! layer above (k+1, bottom-up), over their centre spacing
            ! 0.5*(h(k)+h(k+1)).  Pressure accumulates downward from the
            ! surface for the EOS in-situ derivatives.  Surface layer
            ! (k=nz) has no overlying layer -> N2=0 there.
            p_int = 0.0_wp
            n2_col(nz) = 0.0_wp
            do k = nz - 1, 1, -1
               ! Walk pressure down from the surface to the k/k+1 interface.
               p_int = p_int + GRAVITY*rho0_l*h_col(k + 1)
               t_int = 0.5_wp*(t_col(k) + t_col(k + 1))
               s_int = 0.5_wp*(s_col(k) + s_col(k + 1))
               call eos_specvol_derivs(this%eos, t_int, s_int, p_int, &
                                       dsv_dt_k, dsv_ds_k)
               dbuoy_t = GRAVITY*rho0_l*dsv_dt_k
               dbuoy_s = GRAVITY*rho0_l*dsv_ds_k
               dz_eff = 0.5_wp*(h_col(k) + h_col(k + 1))
               n2_col(k) = (dbuoy_t*(t_col(k + 1) - t_col(k)) + &
                            dbuoy_s*(s_col(k + 1) - s_col(k)))/ &
                           max(dz_eff, H_VANISHED)
               if (n2_col(k) < 0.0_wp) n2_col(k) = 0.0_wp
            end do

            ! ---- Energy input E(x,y) ----
            ! v1: prescribed field.  v1.1 (`e_compute`): state-dependent
            ! internal-tide generation, Jayne & St Laurent (2001):
            !   E = min( 0.5*rho0*kappa_h2*kappa_itides*<h^2>*U_tide^2 * N_bot,
            !            e_max )
            ! with the roughness clamp <h^2> <= (frac_rough*H)^2.  The rho0
            ! factor sets E in W/m^2 (the same units as the prescribed e_in),
            ! so it lands correctly in the rho0*dz*(N^2+Omega^2) Kd divisor.
            ! N_bot = sqrt(N^2) at the bed-most interior interface
            ! (n2_col(kbed), the gradient across the two deepest LIVE layers).
            if (this%e_compute) then
               h2c = min(this%h2_rough, (this%frac_rough*h_tot)**2)
               tke_coef = 0.5_wp*rho0_l*this%kappa_h2*this%kappa_itides*h2c*this%utide**2
               n_bot = sqrt(max(n2_col(kbed), 0.0_wp))
               e_col = tke_coef*n_bot
               e_cap = this%e_max
               if (e_col > e_cap) e_col = e_cap
            else
               e_col = this%e_in(i, j)
            end if
            ! Mask off shallow columns (H < min_zbot).
            if (h_tot < this%min_zbot) e_col = 0.0_wp

            ! ---- Inv_int normalization (Adcroft reciprocal) ----
            ! L'Hospital degenerate-thin-column guard: H/zeta -> 0 => 1.
            hz = h_tot*izeta
            if (hz < 1.0e-14_wp) then
               inv_int = 1.0_wp
            else
               inv_int = 1.0_wp/(1.0_wp - exp(-hz))
            end if

            ! ---- Flux bookkeeping sweep (St-Laurent 2002), bed -> surf ----
            ! TKE_bot = q*mu*E is the total locally-dissipated power; the
            ! 1/rho0 Boussinesq factor lives in the TKE->Kd divisor
            ! (rho0*dz*(N^2+Omega^2)), so Sum(TKE_lay) == q*mu*E exactly
            ! (energy conservation is independent of the rho placement).
            tke_bot = gamma_l*mu_l*e_col
            tke_rem = inv_int*tke_bot
            z_top = 0.0_wp
            do k = 1, kbed - 1
               kd_lay_arr(k) = 0.0_wp   ! inert bed filler: no share of the TKE
            end do
            do k = kbed, nz
               z_top = z_top + h_col(k)
               frac_top = inv_int*exp(-z_top*izeta)
               tke_lay = tke_rem - tke_bot*frac_top
               tke_rem = tke_rem - tke_lay

               ! TKE_to_Kd = 1/(rho0*dz*(N^2+Omega^2)).  dz floored at
               ! H_VANISHED (the vdiff vanishing-layer mass-drop
               ! precedent) and the divisor armoured with H_DIV_EPS.
               dz_eff = max(h_col(k), H_VANISHED)
               denom = rho0_l*dz_eff*(n2_col(k) + omega2_l)
               denom = max(denom, H_DIV_EPS)
               kd_lay = tke_lay/denom
               ! Per-layer physical cap (St-Laurent).  DIVERGENCE D3:
               ! MOM6's max_TKE limiter re-injects the un-used power into
               ! the layer above (energy-conserving under clipping); this
               ! bare Kd_max clamp DISCARDS the clipped power.  v1 accepts
               ! the discard — it is a rarely-active safety cap, not the
               ! deposition mechanism.  The energy-conservation invariant
               ! is asserted on the PRE-clip TKE_lay accordingly.
               if (kd_max_l >= 0.0_wp .and. kd_lay > kd_max_l) kd_lay = kd_max_l
               kd_lay_arr(k) = kd_lay
            end do

            ! ---- D1: exclude the BED and SURFACE LAYERS from deposition ----
            ! MOM6 (St Laurent/Simmons; MOM_set_diffusivity find_TKE_to_Kd,
            ! MOM-inspired) sets TKE_to_Kd(i,1)=TKE_to_Kd(i,nz)=0 — BOTH
            ! endpoint layers get zero Kd_add, and the deposit loop runs the
            ! interior only (top-down do k=nz-1,2,-1).  In MOM6 top-down
            ! indexing layer 1 is the surface and layer nz the bottom; in
            ! Roundabout bottom-up k=nz is the surface and k=1 the bed, so we
            ! zero kd_lay_arr at BOTH ends.  The bed layer (`kbed`: k=1, or
            ! the first LIVE layer `k_bot` under z_fixed) is owned by the BBL
            ! drag; the surface layer (k=nz) is the mixed layer,
            ! owned downstream by EPBL/KPP, and additionally has no overlying
            ! layer so its kernel N^2 is 0 (its Omega^2-only inflated Kd would
            ! over-mix interface K=nz on shallow energetic columns).  Zeroing
            ! the deposited Kd here does NOT touch the TKE energy bookkeeping
            ! above (tke_lay/tke_rem are unchanged) — the energy invariant is
            ! asserted on the pre-deposit TKE_lay, mirroring the original
            ! surface-only exclusion.  Zero both before the 50/50 split.
            kd_lay_arr(kbed) = 0.0_wp
            kd_lay_arr(nz) = 0.0_wp

            ! ---- Deposit per-layer Kd_add 50/50 at the two bounding
            ! interfaces (St-Laurent 2002).  DIVERGENCE D2: MOM6 splits
            ! the layer power 0.5/0.5 across its bounding interfaces; we
            ! mirror that interface deposition.  Global interface K is the
            ! BOTTOM interface of layer K, so layer k contributes to
            ! interfaces k (its bed) and k+1 (its top).
            do k = 1, nz
               this%kd_int(i, j, k) = this%kd_int(i, j, k) + 0.5_wp*kd_lay_arr(k)
               this%kd_int(i, j, k + 1) = this%kd_int(i, j, k + 1) + 0.5_wp*kd_lay_arr(k)
            end do
            ! D1 (cont.): zero the bed (K=kbed) and surface (K=nz+1) interface
            ! end-caps — `vmix_assemble` / EPBL / BBL own them.  With BOTH
            ! the bed LAYER (k=kbed) and surface LAYER (k=nz) now excluded
            ! above, interface K=kbed+1 receives nothing from k=kbed (only
            ! k=kbed+1's bed half)
            ! and interface K=nz nothing from k=nz (only k=nz-1's top half) —
            ! the full MOM6 endpoint exclusion.  DIVERGENCE D4
            ! (BBL N^2 override): MOM6 replaces the near-bed N^2 with a
            ! roughness-height BBL average; v1 uses the raw per-layer N^2
            ! (a v1.1 refinement).
            this%kd_int(i, j, kbed) = 0.0_wp
            this%kd_int(i, j, nz + 1) = 0.0_wp
         end if
      end do
   end subroutine tidal_mixing_column_kernel