vmix_kpp_overlay_impl Subroutine

   pure subroutine vmix_kpp_overlay_impl(grid, this, ms, ss, sf, &
                                         nx_arg, ny_arg, nz_arg, sw_src, &
                                         temp_h, salt_h, sw_active, &
                                         sw_pen_frac, sw_R, sw_zeta1, sw_zeta2, &
                                         sw_method)
      !! KPP boundary-layer overlay on top of the interior closure
      !! already in `this%kv` / `this%kt`.  Phase 1 was shear-driven
      !! only; Phase 2 added the convective velocity scale `w_*` and
      !! γ_T/γ_S non-local transport; Phase 3 (this revision) adds
      !! the V_t² unresolved-turbulence term in the bulk-Ri
      !! denominator (LMD94 eq 23).  Surface BL is now feature-
      !! complete except for Langmuir / Stokes enhancement.
      !!
      !! Two passes per column:
      !!
      !!   1. Bulk-Ri sweep top-down (from the top layer centre) to
      !!      find the boundary-layer depth h_b.  Reference values
      !!      are u_ref, v_ref, B_ref at the top layer centre.  At
      !!      every subsequent layer centre k:
      !!         Ri_bulk(k) = (B_ref - B(k)) * (z_k - z_ref) /
      !!                      max(|V(k) - V_ref|² + V_t²(d_k),
      !!                          shear²_floor)
      !!      with B(k) = -g·ρ(k)/ρ_0.  When Ri_bulk first crosses
      !!      `ri_crit` (default 0.3), linearly interpolate between
      !!      the previous and current layer centres to recover the
      !!      crossing depth.  If we walk to the bed without
      !!      crossing, h_b = total column depth.
      !!
      !!      V_t²(d) = (c_vt2 / ri_crit) · d · N(d) · w_s_col
      !!      adds an unresolved-turbulence contribution that
      !!      sharpens the BL depth diagnosis when grid-scale shear
      !!      is weak (LMD94 eq 23).  N(d) = √(max(0, ΔB/Δd)) is the
      !!      local buoyancy frequency.  w_s_col is computed once
      !!      per column using `this%bl_depth(i, j)` from the
      !!      previous step (lagged h_b) → self-bootstrapping on
      !!      the first call.  Setting `c_vt2 = 0` disables V_t²
      !!      bit-identically.
      !!
      !!   2. KPP overlay: at every interface k with depth from the
      !!      surface d(k) < h_b, compute
      !!         σ = d(k) / h_b
      !!         G(σ) = σ · (1 - σ)²
      !!         w_s = √(u_*² + w_*²)
      !!         kv_kpp = h_b · w_s · G(σ)
      !!      and `max` it into `this%kv(:, :, k)` and
      !!      `this%kt(:, :, k)`.  Outside the BL the interior
      !!      values (PP81 etc.) carry through unchanged.
      !!
      !!   u_* = √(|τ|/ρ_0) is the shear-driven scale.
      !!   w_* = max(0, -B_0·h_b)^(1/3) is the convective scale —
      !!         non-zero only when the surface buoyancy flux is
      !!         destabilizing (B_0 < 0).
      !!   B_0 = (g/ρ_0)·(α_T·F_T - β_S·F_S) where F_T = Q_heat/
      !!         (ρ_0·cp) and F_S = Q_salt/ρ_0 are the kinematic
      !!         surface heat / salt fluxes.  B_0 > 0 = stabilizing
      !!         (heating / freshening), B_0 < 0 = destabilizing
      !!         (cooling / salting).  The `1/ρ_0` is load-bearing:
      !!         `α_T` / `β_S` here are the DIMENSIONAL linear-EOS
      !!         sensitivities (kg m^-3 per degC / psu), not the
      !!         fractional ones — see `kpp_surface_buoyancy_flux`,
      !!         which owns the expression and the convention.  B_0 is
      !!         persisted per column into `this%b0` on the second pass
      !!         (diagnostic; the same quantity as `epbl%b0`).
      !!
      !! `sf` is REQUIRED (A7): B_0 is computed PER COLUMN from the 2D
      !! `sf%Q_heat / Q_salt` fields inside the kernel, so file-driven
      !! spatially-varying forcing (data-override, A3) feeds KPP with no
      !! further change.  Zero-flux behaviour = zero-filled fields.
      !!
      !! Penetrating shortwave (PR-21): when `sw_active`, `B_0` is charged
      !! only for the SW ABSORBED INSIDE the boundary layer.  The heat the
      !! BL actually feels is `Q_bl = Q_heat - I0·T(depth)` (Large,
      !! McWilliams & Doney 1994, App. B), where `I0 = sw_pen_frac·sw_src`
      !! and `T` is the shared two-band `sw_transmission` (Paulson &
      !! Simpson 1977).  `sw_method` selects the reference depth:
      !! `all_sw` ⇒ no correction (legacy), `mxl_sw` ⇒ `T(h_b)`,
      !! `lv1_sw` ⇒ `T(h_layer(nz))`.  `Q_heat`/`Q_salt` stay on `sf`;
      !! the only new array dummy is the explicit-shape `sw_src`.
      !! NOTE (Roundabout KPP fidelity, roadmap trap #4 / PR-11): the overlay
      !! has no Monin-Obukhov length, so `B_0 > 0` ⇒ `w_* ≡ 0` and the SW
      !! method is a numerical no-op under net heating; it bites only in
      !! the sunny-but-net-cooling regime (`B_0 <= 0`, `q_sw > 0`).
      type(hgrid_t), intent(in) :: grid
      type(ocean_vmix_t), intent(inout) :: this
      type(multilayer_state_t), intent(in) :: ms
      type(ocean_surface_stress_t), intent(in) :: ss
      type(ocean_surface_flux_t), intent(in) :: sf
      integer, intent(in) :: nx_arg, ny_arg, nz_arg
         !! Grid extents — declared before `sw_src` (decl-order hook).
      real(wp), intent(in) :: sw_src(nx_arg, ny_arg)
         !! Caller-selected irradiance source (`sf%Q_heat` or `sf%q_sw`),
         !! explicit-shape (per-RK2-stage kernel: no assumed-shape waiver).
      real(wp), intent(in) :: temp_h(nx_arg, ny_arg, nz_arg)
         !! `hTr` of the temperature tracer (degC·m), flattened off the
         !! registry by the shim.  Read ONLY under
         !! `buoyancy_coeffs == BUOY_COEFFS_EOS`, and only at `k = nz`.
      real(wp), intent(in) :: salt_h(nx_arg, ny_arg, nz_arg)
         !! `hTr` of the salinity tracer (PSU·m).  Same contract.
      logical, intent(in) :: sw_active
         !! Host-side `sf%has_sw` gate — false ⇒ `B_0` uses the unmodified
         !! legacy source line (bit-identity).
      real(wp), intent(in) :: sw_pen_frac, sw_R, sw_zeta1, sw_zeta2
         !! Two-band SW parameters (from `sf`), by value.
      integer, intent(in) :: sw_method
         !! `KPP_SW_ALL` | `KPP_SW_MXL` | `KPP_SW_LV1`.

      integer :: i, j, k, nx, ny, nz
      real(wp) :: tau_mag, u_star
      real(wp) :: B_0, wstar3, w_star, w_s
      real(wp) :: u_ref, v_ref, b_ref, d_centre_ref
      real(wp) :: u_k, v_k, b_k, d_centre_k, d_running
      real(wp) :: shear2, ri_bulk, ri_prev, d_prev
      real(wp) :: delta_d, frac, denom, h_b
      real(wp) :: sigma, g_shape, kv_kpp, d_face_k
      real(wp) :: q_T_kin, q_S_kin, gamma_factor
      real(wp) :: h_b_lagged, wstar3_lagged, w_s_col
      real(wp) :: delta_b, n_brunt, vt2
      real(wp) :: q_bl
      real(wp) :: a_buoy, b_buoy, t_sfc, s_sfc, h_sfc, p_buoy, p_buoy_ref
      logical :: crossing_found
      logical :: destabilizing

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

      ! E4 — the pressure the `buoyancy_coeffs = "eos"` α/β are evaluated
      ! at.  `B_0` is a SURFACE buoyancy flux, so the natural pressure is
      ! the one at the TOP of the column:
      !
      !   * `&ocean_psurf_nml in_eos` on  -> `ms%p_top(i,j)`, the ice /
      !     atmospheric load in Pa (the E3 seam), and nothing else.
      !   * off                           -> `eos%p_ref`, the pressure the
      !     model's own `ms%rho_layer` is referenced to, so α stays
      !     consistent with the density field the rest of the closure
      !     differences (a σ₂ run gets its α at 2000 dbar, not at 0).
      !
      ! The gate is `in_eos`, NOT `p_top /= 0`: a cavity run fills `p_top`
      ! with the ice load whether or not `in_eos` is set, and `in_eos` is
      ! the single switch for the whole seam (same note EPBL carries at
      ! its `epbl%in_eos` assignment).  The select below is on a
      ! domain-uniform logical, so it is warp-uniform — free.
      !
      ! This is an α, never a density: it is not differenced along a layer
      ! and so does not violate the "nothing horizontally varying may
      ! enter rho_layer" contract (src/core/ocean/README.md, `p_top` seam).
      p_buoy_ref = this%eos%p_ref

      ! Surface buoyancy flux B_0 is computed PER COLUMN inside the
      ! BL-depth loop from the 2D flux fields (A7): identical arithmetic
      ! per column to the old host scalar when the fields are constant
      ! (bitwise), spatially varying when data-override fills them.

      ! Zero γ before re-populating.
      do concurrent(k=1:nz + 1, j=1:ny, i=1:nx)
         this%gamma_t(i, j, k) = 0.0_wp
         this%gamma_s(i, j, k) = 0.0_wp
      end do

      ! ---- Pass 1: per-column BL depth ----
      do concurrent(j=1:ny, i=1:nx) &
         local(k, u_ref, v_ref, b_ref, d_centre_ref, &
               u_k, v_k, b_k, d_centre_k, &
               shear2, ri_bulk, ri_prev, d_prev, &
               delta_d, frac, denom, h_b, crossing_found, &
               tau_mag, u_star, &
               h_b_lagged, wstar3_lagged, w_s_col, &
               delta_b, n_brunt, vt2, B_0, destabilizing, q_T_kin, q_S_kin, q_bl, &
               a_buoy, b_buoy, t_sfc, s_sfc, h_sfc, p_buoy)
         u_ref = 0.5_wp*(ms%u_face_x_layer(i, j, nz) + ms%u_face_x_layer(i + 1, j, nz))
         v_ref = 0.5_wp*(ms%v_face_y_layer(i, j, nz) + ms%v_face_y_layer(i, j + 1, nz))
         b_ref = -GRAVITY*ms%rho_layer(i, j, nz)/this%rho0
         d_centre_ref = 0.5_wp*ms%h_layer(i, j, nz)

         ! PR-12 dedup: |tau| at cell centres is now a shared field
         ! (ocean_surface_stress_set_derived, same 3-line FP op order as
         ! the inline computation this replaces — bit-identical, §7.5).
         ! Phase 4b: under an ice shelf the wind is masked out of `tau`
         ! and the boundary layer is driven by the ICE-OCEAN stress
         ! instead, which is NOT in `tau` — `stress_shelf` carries it.
         ! The two supports are disjoint (cover mask vs `cover_frac`
         ! weight), so the sum is the total upper-boundary momentum flux.
         ! Always allocated; the zero array without a cavity, and
         ! `x + 0.0` is `x` bit-for-bit.  See the `stress_mag` /
         ! `stress_shelf` contract in `rdb_ocean_surface_stress`.
         tau_mag = ss%stress_mag(i, j) + ss%stress_shelf(i, j)
         u_star = sqrt(tau_mag/this%rho0)
         h_b_lagged = this%bl_depth(i, j)
         ! B_0 charges only the SW absorbed inside the (lagged) BL depth:
         !   Q_bl = Q_heat - I0·T(h_b)  (MXL_SW).  Gated on sw_active so
         !   sw_pen_frac=0 keeps the unmodified legacy source line.
         if (sw_active) then
            q_bl = sf%Q_heat(i, j)
            if (sw_method == KPP_SW_MXL) then
               q_bl = q_bl - sw_pen_frac*sw_src(i, j)* &
                      sw_transmission(h_b_lagged, sw_R, sw_zeta1, sw_zeta2)
            else if (sw_method == KPP_SW_LV1) then
               q_bl = q_bl - sw_pen_frac*sw_src(i, j)* &
                      sw_transmission(ms%h_layer(i, j, nz), sw_R, sw_zeta1, sw_zeta2)
            end if
            q_T_kin = q_bl/(this%rho0*sf%cp)
         else
            q_T_kin = sf%Q_heat(i, j)/(this%rho0*sf%cp)
         end if
         q_S_kin = sf%Q_salt(i, j)/this%rho0
         ! E4 — α/β for B_0.  CONSTANT reproduces the pre-knob line
         ! byte-for-byte (the multiply is on the same two handle members,
         ! via a local, which is an FP no-op); EOS evaluates the ACTIVE
         ! equation of state at this column's surface (T, S) and at the
         ! top-of-column pressure.  A vanished surface layer has no
         ! meaningful (T, S), so it keeps the constants rather than
         ! dividing by ~0 — `H_VANISHED` (dynamic-vanish), not
         ! `H_DIV_EPS`, because that is a real skip, not 1/0 armour.
         a_buoy = this%eos%alpha_T
         b_buoy = this%eos%beta_S
         if (this%buoyancy_coeffs == BUOY_COEFFS_EOS) then
            h_sfc = ms%h_layer(i, j, nz)
            ! vanished-ok: falls back to the CONSTANT `eos%alpha_T`/`beta_S`, not to a zero
            ! concentration: a vanished surface layer must not hand the
            ! buoyancy flux fresh / 0 degC coefficients.
            if (h_sfc > H_VANISHED) then
               t_sfc = temp_h(i, j, nz)/h_sfc
               s_sfc = salt_h(i, j, nz)/h_sfc
               p_buoy = p_buoy_ref
               if (this%p_top_in_eos) p_buoy = ms%p_top(i, j)
               call eos_buoyancy_coeffs(this%eos, t_sfc, s_sfc, p_buoy, &
                                        a_buoy, b_buoy)
            end if
         end if
         B_0 = kpp_surface_buoyancy_flux(a_buoy, b_buoy, this%rho0, &
                                         q_T_kin, q_S_kin)
         destabilizing = (B_0 < 0.0_wp)
         wstar3_lagged = max(0.0_wp, -B_0)*h_b_lagged
         w_s_col = sqrt(u_star*u_star + wstar3_lagged**(2.0_wp/3.0_wp))

         ri_prev = 0.0_wp
         d_prev = d_centre_ref
         h_b = 0.0_wp
         crossing_found = .false.

         d_centre_k = d_centre_ref
         do k = nz - 1, 1, -1
            if (.not. crossing_found) then
               d_centre_k = d_centre_k + 0.5_wp*(ms%h_layer(i, j, k + 1) + ms%h_layer(i, j, k))
               u_k = 0.5_wp*(ms%u_face_x_layer(i, j, k) + ms%u_face_x_layer(i + 1, j, k))
               v_k = 0.5_wp*(ms%v_face_y_layer(i, j, k) + ms%v_face_y_layer(i, j + 1, k))
               b_k = -GRAVITY*ms%rho_layer(i, j, k)/this%rho0

               delta_d = d_centre_k - d_centre_ref
               delta_b = b_ref - b_k

               n_brunt = sqrt(max(0.0_wp, delta_b/delta_d))
               vt2 = (this%c_vt2/this%ri_crit)*d_centre_k*n_brunt*w_s_col

               shear2 = max((u_k - u_ref)*(u_k - u_ref) + (v_k - v_ref)*(v_k - v_ref) + vt2, &
                            this%shear2_floor)
               ri_bulk = delta_b*delta_d/shear2

               if (ri_bulk >= this%ri_crit) then
                  denom = ri_bulk - ri_prev
                  if (abs(denom) > 1.0e-12_wp) then
                     frac = (this%ri_crit - ri_prev)/denom
                  else
                     frac = 0.5_wp
                  end if
                  h_b = d_prev + frac*(d_centre_k - d_prev)
                  crossing_found = .true.
               end if

               ri_prev = ri_bulk
               d_prev = d_centre_k
            end if
         end do

         if (.not. crossing_found) then
            h_b = 0.0_wp
            do k = 1, nz
               h_b = h_b + ms%h_layer(i, j, k)
            end do
         end if

         this%bl_depth(i, j) = h_b
      end do

      ! ---- Pass 2: overlay kv_kpp inside the BL ----
      do concurrent(j=1:ny, i=1:nx) &
         local(k, d_running, d_face_k, sigma, g_shape, kv_kpp, h_b, &
               tau_mag, u_star, &
               wstar3, w_star, w_s, gamma_factor, B_0, destabilizing, &
               q_T_kin, q_S_kin, q_bl, &
               a_buoy, b_buoy, t_sfc, s_sfc, h_sfc, p_buoy)
         ! PR-12 dedup + the ice-shelf stress — see Pass 1's comment.
         tau_mag = ss%stress_mag(i, j) + ss%stress_shelf(i, j)
         u_star = sqrt(tau_mag/this%rho0)
         h_b = this%bl_depth(i, j)
         ! B_0 charges only the SW absorbed inside the (current) BL depth
         !   Q_bl = Q_heat - I0·T(h_b)  (MXL_SW).  See Pass 1.  Gated on
         !   sw_active so sw_pen_frac=0 keeps the legacy source line.
         if (sw_active) then
            q_bl = sf%Q_heat(i, j)
            if (sw_method == KPP_SW_MXL) then
               q_bl = q_bl - sw_pen_frac*sw_src(i, j)* &
                      sw_transmission(h_b, sw_R, sw_zeta1, sw_zeta2)
            else if (sw_method == KPP_SW_LV1) then
               q_bl = q_bl - sw_pen_frac*sw_src(i, j)* &
                      sw_transmission(ms%h_layer(i, j, nz), sw_R, sw_zeta1, sw_zeta2)
            end if
            q_T_kin = q_bl/(this%rho0*sf%cp)
         else
            q_T_kin = sf%Q_heat(i, j)/(this%rho0*sf%cp)
         end if
         q_S_kin = sf%Q_salt(i, j)/this%rho0
         ! E4 — same α/β selection as pass 1; see its comment block.
         a_buoy = this%eos%alpha_T
         b_buoy = this%eos%beta_S
         if (this%buoyancy_coeffs == BUOY_COEFFS_EOS) then
            h_sfc = ms%h_layer(i, j, nz)
            ! vanished-ok: falls back to the CONSTANT `eos%alpha_T`/`beta_S`, not to a zero
            ! concentration: a vanished surface layer must not hand the
            ! buoyancy flux fresh / 0 degC coefficients.
            if (h_sfc > H_VANISHED) then
               t_sfc = temp_h(i, j, nz)/h_sfc
               s_sfc = salt_h(i, j, nz)/h_sfc
               p_buoy = p_buoy_ref
               if (this%p_top_in_eos) p_buoy = ms%p_top(i, j)
               call eos_buoyancy_coeffs(this%eos, t_sfc, s_sfc, p_buoy, &
                                        a_buoy, b_buoy)
            end if
         end if
         B_0 = kpp_surface_buoyancy_flux(a_buoy, b_buoy, this%rho0, &
                                         q_T_kin, q_S_kin)
         this%b0(i, j) = B_0
         destabilizing = (B_0 < 0.0_wp)
         wstar3 = max(0.0_wp, -B_0)*h_b
         w_star = wstar3**(1.0_wp/3.0_wp)
         w_s = sqrt(u_star*u_star + w_star*w_star)
         if (h_b > 0.0_wp .and. w_s > 0.0_wp) then
            d_running = 0.0_wp
            do k = nz, 1, -1
               d_running = d_running + ms%h_layer(i, j, k)
               d_face_k = d_running
               if (k > 1 .and. d_face_k < h_b) then
                  sigma = d_face_k/h_b
                  g_shape = sigma*(1.0_wp - sigma)*(1.0_wp - sigma)
                  kv_kpp = h_b*w_s*g_shape
                  this%kv(i, j, k) = max(this%kv(i, j, k), kv_kpp)
                  this%kt(i, j, k) = max(this%kt(i, j, k), kv_kpp)
                  if (destabilizing) then
                     gamma_factor = this%cs_nonlocal*g_shape
                     this%gamma_t(i, j, k) = gamma_factor*q_T_kin
                     this%gamma_s(i, j, k) = gamma_factor*q_S_kin
                  end if
               end if
            end do
         end if
      end do
   end subroutine vmix_kpp_overlay_impl