diffuse_velocity_columns_impl Subroutine

   pure subroutine diffuse_velocity_columns_impl(nu, nv, nz, dt, &
                                                 u_face, h_layer, kv_centre, &
                                                 wet_cell, &
                                                 x_face, nx_cells, ny_cells, &
                                                 use_harmonic, &
                                                 a_diag, b_diag, c_diag, rhs, &
                                                 do_stress, do_drag, rho0, &
                                                 tau_face, lambda_bot, &
                                                 do_top, lambda_top, cover_face, &
                                                 solve_momentum, &
                                                 do_remnant, visc_rem_out, &
                                                 hvel_mom6, hbbl_visc, &
                                                 bbl_glue, hvel_upwind, hvel_harmonic, &
                                                 kv_bbl_face, bbl_thick_face, kv_bbl_bg, &
                                                 do_corner, kv_prandtl, kv_corner, &
                                                 zlevel_faces, k_top_face, k_bot_face)
      !! Build + solve the tridiagonal system per face column.
      !! `u_face` is either `u_face_x_layer` (`x_face = .true.`,
      !! shape (nx+1, ny)) or `v_face_y_layer` (`x_face = .false.`,
      !! shape (nx, ny+1)).  `h_layer` is cell-centred (nx, ny, nz)
      !! and we average across the face direction to get the face
      !! thickness.  `kv_centre` is the same cell-centred diffusivity
      !! field used by the tracer kernel — we average it across the
      !! face to get a face-located value at each interface.
      !!
      !! For wall faces (i = 1 / i = nx_cells+1 for x_face;
      !! j = 1 / j = ny_cells+1 for y-face), there's no neighbour
      !! cell to average with — we fall back to the single available
      !! cell's thickness + diffusivity.  Mass through walls is
      !! forced to zero by the continuity kernel anyway, so the
      !! wall-face viscosity is only there to keep the system
      !! non-singular.
      integer, intent(in) :: nu, nv, nz
      integer, intent(in) :: nx_cells, ny_cells
      real(wp), intent(in) :: dt
      real(wp), intent(inout) :: u_face(nu, nv, nz)
      real(wp), intent(in) :: h_layer(nx_cells, ny_cells, nz)
      real(wp), intent(in) :: kv_centre(nx_cells, ny_cells, nz + 1)
      real(wp), intent(in) :: wet_cell(nx_cells, ny_cells)
         !! Cell-centred wet/dry land mask (1 wet, 0 land).  The wind-stress
         !! fold multiplies `tau_face` by the face mask `min` of the two
         !! bounding cells — matching the explicit surface-stress kernel —
         !! so no stress is injected at a no-normal-flow land face.  All-wet
         !! (mask ≡ 1) ⇒ bit-identical.  Only read when `do_stress`.
      logical, intent(in) :: x_face
      logical, intent(in) :: use_harmonic
      real(wp), intent(inout) :: a_diag(nu, nv, nz)
      real(wp), intent(inout) :: b_diag(nu, nv, nz)
      real(wp), intent(inout) :: c_diag(nu, nv, nz)
      real(wp), intent(inout) :: rhs(nu, nv, nz)
      logical, intent(in) :: do_stress
         !! Fold the surface wind stress into the `k = nz` RHS row.
      logical, intent(in) :: do_drag
         !! Fold the bottom drag into the `k = k_bot_face` diagonal.
      real(wp), intent(in) :: rho0
         !! Boussinesq reference density for the `τ/ρ₀` stress conversion.
      real(wp), intent(in), optional :: tau_face(nu, nv)
         !! Wind stress (N/m²) on this face.  Present iff `do_stress`.
      real(wp), intent(in), optional :: lambda_bot(nu, nv)
         !! Bottom-drag Rayleigh rate λ (1/s) on this face.  Present iff
         !! `do_drag`.
      logical, intent(in) :: do_top
         !! Fold the ice-shelf top drag into the `k = nz` diagonal AND
         !! mask the wind-stress RHS by `(1 − cover_face)`.
      real(wp), intent(in), optional :: lambda_top(nu, nv)
         !! Top-drag Rayleigh rate λ (1/s) on this face.  Present iff
         !! `do_top`.
      real(wp), intent(in), optional :: cover_face(nu, nv)
         !! Face ice-cover mask (0 open, 1 under ice), the OR of the two
         !! abutting cells (`rdb_ocean_top_drag`).  Present iff `do_top`.
      logical, intent(in) :: hvel_mom6
         !! MOM6 HARMONIC_VISC parity for `h_u` + `h_shear` (see the slot-type
         !! docstring).  `.false.` => the historical arithmetic-`h_u` /
         !! `face_thick`-`dz` pair, bit-identical.
      real(wp), intent(in) :: hbbl_visc
         !! Bottom-layer scale for the `botfn` blend (MOM6 HBBL).
      logical, intent(in) :: bbl_glue
         !! MOM6 `bottomdraglaw` coupling parity (see the slot-type
         !! docstring).  Configure guarantees `hvel_mom6` when this is on;
         !! the bed piston then applies whether or not `do_drag`.
         !! `.false.` ⇒ bit-identical.
      logical, intent(in) :: hvel_upwind
         !! `.false.` = skip the near-bed upwind blend (pure harmonic
         !! hvel).  See the slot-type docstring.  `hvel_harmonic` branch only.
      logical, intent(in) :: hvel_harmonic
         !! MOM6 `HARMONIC_VISC`: `.false.` = MOM6's default `z_clear`
         !! face-thickness branch, `.true.` = the harmonic branch (slot-type
         !! docstring).  Only read when `hvel_mom6`.
      real(wp), intent(in) :: kv_bbl_face(nu, nv)
         !! BBL viscosity `kv_bbl` (m²/s) of this face (MOM6
         !! `visc%Kv_bbl_u`).  Only read when `bbl_glue`.
      real(wp), intent(in) :: bbl_thick_face(nu, nv)
         !! BBL thickness (m) of this face (MOM6 `visc%bbl_thick_u`); also
         !! the height-above-bed normaliser `1/I_Hbbl` under the glue, as in
         !! MOM6's `bottomdraglaw` branch.  Only read when `bbl_glue`.
      real(wp), intent(in) :: kv_bbl_bg
         !! MOM6 `KV`: the background viscosity the BBL viscosity replaces
         !! within botfn reach of the bed, `Kv_tot + (kv_bbl − KV)·botfn`.
      logical, intent(in) :: solve_momentum
         !! `.false.` = remnant-only mode: build the matrix and compute
         !! `visc_rem_out` WITHOUT touching `u_face` (the pre-substep
         !! visc_rem refresh, PGF_BUG.md §9 — MOM6 computes vertvisc_coef
         !! before btstep every stage, so its BT weights never lag; the
         !! stage-end producer alone leaves visc_rem ≡ 1 for the whole
         !! first stage, and the Δu corrector then deposits the spurious
         !! column-mean into wet layers).  `.true.` = the normal solve.
      logical, intent(in) :: do_remnant
         !! Fill `visc_rem_out` with the viscous remnant γ_k — the
         !! sensitivity of the post-friction layer velocity to a uniform
         !! barotropic acceleration, `γ_k ≡ (1/Δt)·∂u_k^{n+1}/∂Ā`.  γ
         !! solves the SAME tridiagonal system as the momentum solve
         !! (same `a_diag`/`b_diag`, and the already-factorized `c_diag`
         !! left behind by the momentum forward sweep) with RHS ≡ 1 —
         !! Roundabout's rows are pre-normalized by `h_k`, so the remnant RHS
         !! is the unit vector, NOT `h_u(k)` (MOM6's un-normalized
         !! convention).  Reusing the surviving factorization means γ
         !! cannot drift from the operator actually solved for momentum.
      real(wp), intent(inout), optional :: visc_rem_out(nu, nv, nz)
         !! Output γ_k, bottom-up (`k = 1` bed → γ smallest, `k = nz`
         !! surface → γ → 1).  Clamped to `min(γ, 1.0)` at production
         !! (cheap FP insurance on a quantity the maximum principle
         !! already bounds to (0, 1]).  Present iff `do_remnant`.
      integer, intent(in) :: k_top_face(nu, nv)
         !! `ms%k_top_u` / `k_top_v` — the index of the first layer that
         !! is LIVE on BOTH sides of this face, `nz` wherever nothing
         !! vanishes against the top.  The surface row of the column
         !! solve, the ice-shelf top-drag diagonal fold and the
         !! wind-stress RHS fold all sit on THIS row instead of on `nz`:
         !! under a quasi-geopotential coordinate beneath a shelf the
         !! rows above it are inert fillers, decoupled to the identity by
         !! the `zlevel_faces` gates, so a `dt*lambda_top` added at `nz`
         !! would damp a row that is not coupled to the ocean at all —
         !! and the drag rate `lambda_top` is itself captured at
         !! `k_top_face` by `rdb_ocean_top_drag`, so the two MUST agree.
         !! With `k_top_face ≡ nz` every expression below is the one it
         !! replaced, bit for bit.
      integer, intent(in) :: k_bot_face(nu, nv)
         !! `ms%k_bot_u` / `k_bot_v` — the first layer LIVE on BOTH sides
         !! of this face counting UP from the bed (`max` of its two
         !! columns), `1` wherever nothing vanishes against the bed.  The
         !! bed row of the column solve (no-flux-below BC, the implicit
         !! bottom-drag diagonal fold and the `bbl_glue` piston) sits on
         !! THIS row instead of on `k = 1`, and the rows below it — inert
         !! `z_fixed` bed fillers — are the identity.  The drag rate
         !! `lambda_bot` is itself captured at `k_bot_face` by
         !! `ocean_bottom_drag_compute_tendencies`, so the two MUST agree.
         !! The height-above-bed stack `zint` (BBL glue, upwind blend)
         !! also starts accumulating here.  With `k_bot_face ≡ 1` every
         !! expression below is the one it replaced, bit for bit.
      logical, intent(in) :: zlevel_faces
         !! `&vcoord_nml zfixed_closed_faces` — z-level partial steps.
         !!
         !! Changes TWO things about the face column, and only when on
         !! (`.false.` ⇒ every expression below is textually and
         !! bit-identically what it has always been):
         !!
         !! 1. the face thickness becomes `min(h_L, h_R)` instead of the
         !!    arithmetic (or MOM6-harmonic) mean — the standard
         !!    partial-cell choice, and the one that makes "this layer is
         !!    an inert FILLER on at least one side" visible locally as
         !!    `hvel(k) <= H_VANISHED`;
         !! 2. the tridiagonal coupling across an interface touching such
         !!    a layer is CUT (`alpha = beta = 0`), and the layer's own
         !!    row is forced to the identity.
         !!
         !! Without (2) an OPEN layer is frictionally coupled, every
         !! single step, to the CLOSED layer above or below it — whose
         !! velocity `mask_layer_velocities` has just zeroed — so the
         !! wall acts as a spurious side drag instead of as free-slip.
         !! The floor `max(hvel, H_VANISHED)` does not save it: it makes
         !! the coupling huge-but-finite (`dt·nu/1.5e-4`), i.e. a rigid
         !! glue to zero, which is the worst of the three options.
         !!
         !! This is the DYNAMIC twin of `metrics%open_u/open_v`: "filler
         !! on either side" is `min(h_L,h_R) <= H_VANISHED` here and
         !! `target_h <= H_VANISHED` on either side there, and under
         !! `z_fixed` the two agree (η is absorbed by the first LIVE
         !! layer).  It is derived locally rather than threaded as a 3-D
         !! mask because this kernel is reached through four flat-impl
         !! call sites and five dispatcher call sites, none of which sees
         !! `ocean_metrics_t`, and a scalar costs nothing.
         !! The equivalent tracer decoupling already exists, unconditionally,
         !! in `build_factorize_tracer_matrix`.
      logical, intent(in) :: do_corner
         !! Add the corner-staggered viscosity to every face interface.
         !! Gated INSIDE the single DC (no split loop — NVHPC penalty);
         !! `kv_corner` is guaranteed present when `do_corner`.
      real(wp), intent(in) :: kv_prandtl
         !! Scale on `kv_corner` (Kv = Pr·Kd).  Only read when
         !! `do_corner`.
      real(wp), intent(in), optional :: kv_corner(nx_cells + 1, ny_cells + 1, nz + 1)
         !! Corner-staggered interface viscosity (SW-corner convention:
         !! corner (i,j) is the SW corner of cell (i,j)).  A face reads
         !! the 2-point average of its two END corners: u-face (i,j) →
         !! corners (i,j)/(i,j+1); v-face (i,j) → corners (i,j)/(i+1,j)
         !! — the direct corner→face route, never via a tracer point.
         !! Present iff `do_corner`.

      integer :: i, j, k, ktop, kb, i_left, i_right, j_below, j_above, i_c2, j_c2
      real(wp) :: hf_km1, hf_k, hf_kp1, dz_bot, dz_top
      real(wp) :: hvel(NZ_STACK_MAX)
      real(wp) :: zint(NZ_STACK_MAX)
         !! Normalized height (units of `hbbl_visc`) of the TOP interface
         !! of layer k above the bed, accumulated from HARMONIC face
         !! thicknesses — the mirror of MOM6's `z_i`.
         !! Grounded sliver stacks therefore sit at zint ≈ 0 no matter how
         !! thick their arithmetic-mean face layers are.  Only filled when
         !! `hvel_mom6` (configure guarantees that whenever `bbl_glue`).
      real(wp) :: zacc, z2, botfn, h_harm, h_arith, h_delta, hl_c, hr_c, i_hbbl
      real(wp) :: zcol_l, zcol_r, d_min, z_clear
      real(wp) :: nu_face_k, nu_face_kp1, alpha, beta, denom
      real(wp) :: botfn_int, kv_bbl, bbl_thick
      real(wp) :: inv_rho0, wind_open
      real(wp), parameter :: EPS_HVEL = 1.0e-30_wp
         !! MOM6 `h_neglect` analogue in the harmonic mean / HBBL inverse.
      real(wp), parameter :: STRESS_H_MIN = 1.0e-3_wp
         !! Floor on the surface-layer face thickness in the `τ/(ρ₀·h)`
         !! stress conversion — matches `ocean_surface_stress_t%h_min` so
         !! the implicit fold reduces to the explicit kernel in the
         !! well-resolved limit.

      inv_rho0 = 1.0_wp/rho0

      do concurrent(j=1:nv, i=1:nu) &
         local(i_left, i_right, j_below, j_above, i_c2, j_c2, &
               hf_km1, hf_k, hf_kp1, dz_bot, dz_top, &
               nu_face_k, nu_face_kp1, alpha, beta, denom, k, ktop, kb, &
               hvel, zint, zacc, z2, botfn, botfn_int, kv_bbl, bbl_thick, &
               h_harm, h_arith, h_delta, hl_c, hr_c, wind_open, &
               zcol_l, zcol_r, d_min, z_clear)
         ! Neighbour cell indices for averaging h.  (i_c2, j_c2) is the
         ! SECOND end corner of this face for the corner-viscosity
         ! add-on (the first is (i, j) on both staggerings): a u-face
         ! runs south→north (corners (i,j)/(i,j+1)), a v-face west→east
         ! (corners (i,j)/(i+1,j)).
         if (x_face) then
            i_left = max(1, i - 1)
            i_right = min(nx_cells, i)
            j_below = j
            j_above = j
            i_c2 = i
            j_c2 = j + 1
         else
            i_left = i
            i_right = i
            j_below = max(1, j - 1)
            j_above = min(ny_cells, j)
            i_c2 = i + 1
            j_c2 = j
         end if

         ! First live layer of this face counting up from the bed (`1`
         ! everywhere but a `z_fixed` bed filler band).
         kb = k_bot_face(i, j)

         ! Seed RHS with u^n.
         do k = 1, nz
            rhs(i, j, k) = u_face(i, j, k)
         end do

         ! ---- This face's bottom boundary layer (MOM6 visc%Kv_bbl_u,
         ! visc%bbl_thick_u).  Under the glue the height above the bed is
         ! normalised by the BBL thickness, as MOM6's bottomdraglaw branch
         ! does (`I_Hbbl = 1/bbl_thick`, vertvisc_coef); otherwise by
         ! `hbbl_visc` (MOM6 HBBL).  The legacy glue fills `bbl_thick_face`
         ! with `hbbl_visc`, so its arithmetic is unchanged.
         kv_bbl = 0.0_wp
         bbl_thick = hbbl_visc
         if (bbl_glue) then
            kv_bbl = kv_bbl_face(i, j)
            bbl_thick = bbl_thick_face(i, j)
         end if
         i_hbbl = 1.0_wp/(bbl_thick + EPS_HVEL)

         ! ---- Face thickness per layer (MOM6 `hvel`) ----
         ! k = 1 is the BED here (MOM6 counts from the surface), so the height
         ! above bed accumulates UPWARD.
         zacc = 0.0_wp
         ! MOM6 `HARMONIC_VISC = False` (`z_clear`) branch: the two columns'
         ! interface heights, from their beds up.  MOM6 starts them at the
         ! STATIC depth `-bathyT`; here the column sums stand in for it, which
         ! differs only by the free surface and cancels in `z_clear` wherever
         ! the two cells' eta agree (an O(d eta) difference across a face).
         zcol_l = 0.0_wp
         zcol_r = 0.0_wp
         if (hvel_mom6 .and. .not. hvel_harmonic) then
            do k = 1, nz
               zcol_l = zcol_l - h_layer(i_left, j_below, k)
               zcol_r = zcol_r - h_layer(i_right, j_above, k)
            end do
         end if
         d_min = min(-zcol_l, -zcol_r)
         do k = 1, nz
            hl_c = h_layer(i_left, j_below, k)
            hr_c = h_layer(i_right, j_above, k)
            if (hvel_mom6 .and. .not. hvel_harmonic) then
               ! MOM6 vertvisc_coef, HARMONIC_VISC = False branch.  `zcol_*`
               ! is the height of the TOP of layer k in each cell; `z_clear`
               ! is the height of the higher of the two above the SHALLOWER
               ! bed (`-d_min`), so a
               ! face layer below the shallower bed of a step sits at
               ! z_clear <= 0 and only `zh` (harmonic, ~0 against a filler)
               ! lifts it.  `harm_BL_val = 0` (MOM6 HARMONIC_BL_SCALE
               ! default) ⇒ `z2_wt = 1`.
               h_harm = 2.0_wp*hl_c*hr_c/(hl_c + hr_c + EPS_HVEL)
               h_arith = 0.5_wp*(hl_c + hr_c)
               h_delta = hr_c - hl_c
               zcol_l = zcol_l + hl_c
               zcol_r = zcol_r + hr_c
               if (k >= kb) zacc = zacc + h_harm
               z_clear = max(zcol_l, zcol_r) + d_min
               zint(k) = max(zacc, z_clear)*i_hbbl
               hvel(k) = h_arith
               ! Flow from the THIN side to the THICK side: near the bed the
               ! face takes the harmonic (thin-side) thickness, so a
               ! near-massless layer cannot be drained through a thick face.
               if (u_face(i, j, k)*h_delta > 0.0_wp) then
                  z2 = zint(k)
                  botfn = 1.0_wp/(1.0_wp + 0.09_wp*z2*z2*z2*z2*z2*z2)
                  hvel(k) = (1.0_wp - botfn)*h_arith + botfn*h_harm
               end if
            else if (hvel_mom6) then
               h_harm = 2.0_wp*hl_c*hr_c/(hl_c + hr_c + EPS_HVEL)
               h_arith = 0.5_wp*(hl_c + hr_c)
               h_delta = hr_c - hl_c
               hvel(k) = h_harm
               ! Upwind bias: only when the flow runs from the THICK side to
               ! the THIN side does the near-bed face take the arithmetic
               ! (donor) thickness.  Without this the harmonic mean
               ! over-suppresses asymmetrically -- the reason bare-harmonic
               ! attempts regressed (LAGRANGIAN_PGF_BUG.md 6.2).
               !
               ! `hvel_upwind = .false.` disables the blend (pure harmonic
               ! hvel).  The sign test keys on u itself, so at (near-)rest it
               ! flip-flops faces between harmonic and arithmetic on
               ! ROUNDOFF-sign velocities, collapsing the BBL glue at
               ! whichever faces flip — measured ×200-800/stage residual
               ! amplification on the rest-state reproducer (PGF_BUG.md
               ! §9.8).  MOM6 carries the same test and the same hazard; it
               ! never bites there only because its state stays at 1e-14.
               if (hvel_upwind .and. u_face(i, j, k)*h_delta < 0.0_wp) then
                  z2 = zacc
                  botfn = 1.0_wp/(1.0_wp + 0.09_wp*z2*z2*z2*z2*z2*z2)
                  hvel(k) = (1.0_wp - botfn)*h_harm + botfn*h_arith
               end if
               ! Height above the LIVE bed: the inert fillers below `kb`
               ! do not count (bit-identical at kb = 1).
               if (k >= kb) zacc = zacc + h_harm*i_hbbl
               zint(k) = zacc
            else
               hvel(k) = 0.5_wp*(hl_c + hr_c)
            end if
            ! z-level partial steps: the face column is the OVERLAP of the
            ! two cell columns, so its thickness is the min.  A filler on
            ! either side then reads at or below H_VANISHED and the
            ! coupling gates below cut it out of the solve.
            if (zlevel_faces) hvel(k) = min(hl_c, hr_c)
         end do

         ! ---- Rows BELOW the first live layer: the identity ----
         ! Empty whenever `kb = 1`, so bit-identical everywhere a face has
         ! no bed-side filler.  Where it is not empty those rows are inert
         ! `z_fixed` fillers inside the bed: no mass, not coupled to
         ! anything, and `rhs` holds `u^n`, so the Thomas sweep returns
         ! them unchanged (the mirror of the identity rows above `k_top`).
         do k = 1, kb - 1
            a_diag(i, j, k) = 0.0_wp
            b_diag(i, j, k) = 1.0_wp
            c_diag(i, j, k) = 0.0_wp
         end do

         ! ---- k = kb: bed BC, no flux below ----
         ! Floor the face thicknesses at H_VANISHED before they enter the
         ! α/β denominators (dt·ν/(hf_k·dz)).  For an exactly-collapsed layer
         ! (both neighbour cells h=0 ⇒ hvel=0) with kv>0, a raw hf_k=0 (and the
         ! face_thick-derived dz=0) would make α=Inf ⇒ b_diag=Inf ⇒ NaN γ.
         ! Mirrors the tracer path's max(h, H_VANISHED) floor.  dz_top/dz_bot
         ! are built FROM these floored hf_* so they inherit the floor.  In
         ! every valid config hvel ≫ H_VANISHED ⇒ the max() is a no-op ⇒
         ! bit-identical.  Floor at the READ site (not by reassigning hvel(k))
         ! — gfortran mis-optimizes a local() array element reassigned across
         ! branches.
         hf_k = max(hvel(kb), H_VANISHED)
         ! SINGLE-LAYER COLUMN (nz = 1, or kb = nz: one live layer).  The
         ! bed row IS the surface row: there is no interior interface above it, so the interior
         ! coupling is identically zero and the surface is a pure
         ! stress-Neumann BC (a RHS source, added below with the k = nz
         ! block's).  Without this gate the code reads `hvel(2)` — one
         ! past the layer loop that filled `hvel`, i.e. uninitialised
         ! `local()` stack — and the k = nz block below then OVERWRITES
         ! this row from `hvel(0)`/`zint(0)`, an out-of-bounds read that
         ! is a DETERMINISTIC `CUDA_ERROR_ILLEGAL_ADDRESS` on the NVHPC
         ! GPU build and silent stack garbage on the host.  Gated INSIDE
         ! the single `do concurrent` (a split loop costs an extra launch
         ! — see CLAUDE.md) and loop-invariant, so nz >= 2 is
         ! bit-identical.  Written as init-then-conditional-overwrite, not
         ! if/else: gfortran 15.1 miscompiles a `local()` scalar reassigned
         ! across the two arms of an if/else inside `do concurrent`.
         alpha = 0.0_wp
         if (kb < nz) then
            hf_kp1 = max(hvel(kb + 1), H_VANISHED)
            if (hvel_mom6) then
               dz_top = 0.5_wp*(hf_k + hf_kp1)   ! MOM6 h_shear: arithmetic of hvels
            else
               dz_top = face_thick(hf_k, hf_kp1, use_harmonic)
            end if
            nu_face_kp1 = 0.5_wp*(kv_centre(i_left, j_below, kb + 1) + &
                                  kv_centre(i_right, j_above, kb + 1))
            ! Corner-viscosity add-on (vertex kappa-shear Kv seam): direct
            ! 2-point end-corner average onto this face, BEFORE the BBL
            ! glue (all viscosity contributions fold ahead of the
            ! coupling-coefficient blend).
            if (do_corner) then
               nu_face_kp1 = nu_face_kp1 + &
                             kv_prandtl*(0.5_wp*(kv_corner(i, j, kb + 1) + kv_corner(i_c2, j_c2, kb + 1)))
            end if
            ! BBL glue at the interface above the bed layer (MOM6
            ! find_coupling_coef): within botfn reach of the bed the BBL
            ! viscosity REPLACES the background one, `Kv_tot +
            ! (kv_bbl − KV)·botfn` (shear / boundary-layer
            ! contributions kept; floored at 0 for a `kv_bbl` below `KV`,
            ! where MOM6 relies on Kv_tot >= KV), and the shear distance is
            ! capped toward bbl_thick — grounded stacks (zint ≈ 0) become
            ! rigidly coupled, which is the mechanism that absorbs the
            ! spurious grounded-layer PGF every step (PGF_BUG.md §9).
            if (bbl_glue) then
               botfn_int = 1.0_wp/(1.0_wp + 0.09_wp*zint(kb)*zint(kb)*zint(kb)* &
                                   zint(kb)*zint(kb)*zint(kb))
               nu_face_kp1 = max(0.0_wp, nu_face_kp1 + (kv_bbl - kv_bbl_bg)*botfn_int)
               if (dz_top > bbl_thick) then
                  dz_top = (1.0_wp - botfn_int)*dz_top + botfn_int*bbl_thick
               end if
            end if
            alpha = dt*nu_face_kp1/(hf_k*dz_top)
            ! z-level partial steps: no shear across an interface that
            ! touches a CLOSED face layer.  The tracer twin does exactly
            ! this, unconditionally, in `build_factorize_tracer_matrix`.
            if (zlevel_faces) then
               ! vanished-ok: row decoupling across a z-level wall (the closed-face rule) —
               ! a momentum-grid question, not a tracer concentration.
               if (hvel(kb) <= H_VANISHED .or. hvel(kb + 1) <= H_VANISHED) alpha = 0.0_wp
            end if
         end if
         a_diag(i, j, kb) = 0.0_wp
         c_diag(i, j, kb) = -alpha
         b_diag(i, j, kb) = 1.0_wp + alpha
         ! Bottom-drag stress BC (Roundabout bed = k=kb; MIRROR of MOM6 k=nz).
         ! λ_bot is the Rayleigh RATE (c_d·|U|/h_1 quadratic, r linear) the
         ! bottom-drag slot already forms — the row is pre-normalized by
         ! h_1, so the MOM6 `h_1 + dt·a_bot` diagonal becomes `1 + dt·λ_bot`
         ! here.  A drag is a SINK ⇒ POSITIVE diagonal add ⇒ |amplification|
         ! ≤ 1 for ANY h/dt (unconditionally stable on thin shelf bottoms).
         ! Gated INSIDE the single DC (no split loop — NVHPC penalty); the
         ! optional `lambda_bot` is guaranteed present when `do_drag`.
         if (bbl_glue) then
            ! MOM6 bed coupling (`a_cpl(nz+1)`): the drag is a viscous PISTON
            ! `kv_bbl/(min(hvel₁/2, bbl_thick))` — it DIVERGES as the bottom
            ! layer thins (a sliver bed layer is anchored rigidly to rest),
            ! where the Rayleigh `dt·λ` fold below is h-independent and lets
            ! grounded stacks reach the ballistic balance u_eq = PGF·h/r
            ! (PGF_BUG.md §9.4).  With the per-face BBL `kv_bbl =
            ! sqrt(CDRAG)·u*·bbl_thick`, so a resolved bed layer (hvel₁ >=
            ! 2·bbl_thick) carries the stress `CDRAG·u_bbl·u₁`.  It IS the
            ! bed sink whenever the glue is on — `implicit_drag` or not —
            ! and replaces (not augments) the λ fold; the driver skips the
            ! explicit drag apply.
            b_diag(i, j, kb) = b_diag(i, j, kb) + dt*kv_bbl/ &
                               (hf_k*(min(0.5_wp*hvel(kb), bbl_thick) + EPS_HVEL))
         else if (do_drag) then
            b_diag(i, j, kb) = b_diag(i, j, kb) + dt*lambda_bot(i, j)
         end if

         ! ---- k = kb+1..nz-1 ----
         do k = kb + 1, nz - 1
            ! Floor the face thicknesses (see the k=1 block) — keeps the α/β
            ! denominators non-zero for a collapsed interior layer; no-op for
            ! hvel ≫ H_VANISHED ⇒ bit-identical.
            hf_km1 = max(hvel(k - 1), H_VANISHED)
            hf_k = max(hvel(k), H_VANISHED)
            hf_kp1 = max(hvel(k + 1), H_VANISHED)
            if (hvel_mom6) then
               dz_bot = 0.5_wp*(hf_km1 + hf_k)
               dz_top = 0.5_wp*(hf_k + hf_kp1)
            else
               dz_bot = face_thick(hf_km1, hf_k, use_harmonic)
               dz_top = face_thick(hf_k, hf_kp1, use_harmonic)
            end if
            nu_face_k = 0.5_wp*(kv_centre(i_left, j_below, k) + kv_centre(i_right, j_above, k))
            nu_face_kp1 = 0.5_wp*(kv_centre(i_left, j_below, k + 1) + &
                                  kv_centre(i_right, j_above, k + 1))
            ! Corner-viscosity add-on (see the k=1 block).
            if (do_corner) then
               nu_face_k = nu_face_k + &
                           kv_prandtl*(0.5_wp*(kv_corner(i, j, k) + kv_corner(i_c2, j_c2, k)))
               nu_face_kp1 = nu_face_kp1 + &
                             kv_prandtl*(0.5_wp*(kv_corner(i, j, k + 1) + &
                                                 kv_corner(i_c2, j_c2, k + 1)))
            end if
            ! BBL glue (see the k=1 block).  The (nu, dz) transform uses the
            ! interface's own zint, so the SAME effective coupling lands in
            ! this row's alpha and the row-above's beta (symmetric matrix).
            if (bbl_glue) then
               botfn_int = 1.0_wp/(1.0_wp + 0.09_wp*zint(k - 1)*zint(k - 1)*zint(k - 1)* &
                                   zint(k - 1)*zint(k - 1)*zint(k - 1))
               nu_face_k = max(0.0_wp, nu_face_k + (kv_bbl - kv_bbl_bg)*botfn_int)
               if (dz_bot > bbl_thick) then
                  dz_bot = (1.0_wp - botfn_int)*dz_bot + botfn_int*bbl_thick
               end if
               botfn_int = 1.0_wp/(1.0_wp + 0.09_wp*zint(k)*zint(k)*zint(k)* &
                                   zint(k)*zint(k)*zint(k))
               nu_face_kp1 = max(0.0_wp, nu_face_kp1 + (kv_bbl - kv_bbl_bg)*botfn_int)
               if (dz_top > bbl_thick) then
                  dz_top = (1.0_wp - botfn_int)*dz_top + botfn_int*bbl_thick
               end if
            end if
            beta = dt*nu_face_k/(hf_k*dz_bot)
            alpha = dt*nu_face_kp1/(hf_k*dz_top)
            ! z-level partial steps — see the bed row.
            if (zlevel_faces) then
               ! vanished-ok: row decoupling across a z-level wall (the closed-face rule) —
               ! a momentum-grid question, not a tracer concentration.
               if (hvel(k) <= H_VANISHED .or. hvel(k - 1) <= H_VANISHED) beta = 0.0_wp
               ! vanished-ok: row decoupling across a z-level wall (the closed-face rule) —
               ! a momentum-grid question, not a tracer concentration.
               if (hvel(k) <= H_VANISHED .or. hvel(k + 1) <= H_VANISHED) alpha = 0.0_wp
            end if
            a_diag(i, j, k) = -beta
            c_diag(i, j, k) = -alpha
            b_diag(i, j, k) = 1.0_wp + alpha + beta
         end do

         ! ---- k = k_top: surface BC, no flux above ----
         ! `k_top` is `nz` on every column with no top-side filler, so
         ! this block is the historical `k = nz` block verbatim there.
         ! Floor the face thicknesses (see the k=1 block) — keeps the β
         ! denominator non-zero for a collapsed surface layer; no-op for
         ! hvel ≫ H_VANISHED ⇒ bit-identical.  (The stress-BC conversion
         ! below keeps its own STRESS_H_MIN floor on hf_k.)
         ktop = k_top_face(i, j)
         hf_k = max(hvel(ktop), H_VANISHED)
         ! SINGLE-LAYER COLUMN (nz = 1): this row has already been built as
         ! the bed row above (a = c = 0, b = 1 + bottom drag).  Skip the
         ! interface-below build entirely — at nz = 1 it would read
         ! `hvel(0)` / `zint(0)` (out of bounds) and overwrite the bed
         ! row's drag.  Only the stress RHS below still applies, which is
         ! exactly right: one layer carries BOTH the wind stress and the
         ! bottom drag.  The gate is `ktop > kb` (was `ktop > 1`), which
         ! also covers a face with ONE live layer (`kb = ktop`, a z_fixed
         ! column shallower than one nominal layer): the bed row built
         ! above is then the surface row too.  kb = 1 ⇒ bit-identical.
         if (ktop > kb) then
            hf_km1 = max(hvel(ktop - 1), H_VANISHED)
            if (hvel_mom6) then
               dz_bot = 0.5_wp*(hf_km1 + hf_k)
            else
               dz_bot = face_thick(hf_km1, hf_k, use_harmonic)
            end if
            nu_face_k = 0.5_wp*(kv_centre(i_left, j_below, ktop) &
                                + kv_centre(i_right, j_above, ktop))
            ! Corner-viscosity add-on (see the k=1 block).
            if (do_corner) then
               nu_face_k = nu_face_k + &
                           kv_prandtl*(0.5_wp*(kv_corner(i, j, ktop) + kv_corner(i_c2, j_c2, ktop)))
            end if
            ! BBL glue (see the k=1 block) — this row's beta is the twin of
            ! row ktop-1's alpha and must see the same transform.
            if (bbl_glue) then
               botfn_int = 1.0_wp/(1.0_wp + 0.09_wp*zint(ktop - 1)*zint(ktop - 1)*zint(ktop - 1)* &
                                   zint(ktop - 1)*zint(ktop - 1)*zint(ktop - 1))
               nu_face_k = max(0.0_wp, nu_face_k + (kv_bbl - kv_bbl_bg)*botfn_int)
               if (dz_bot > bbl_thick) then
                  dz_bot = (1.0_wp - botfn_int)*dz_bot + botfn_int*bbl_thick
               end if
            end if
            beta = dt*nu_face_k/(hf_k*dz_bot)
            ! z-level partial steps — see the bed row.
            if (zlevel_faces) then
               ! vanished-ok: row decoupling across a z-level wall (the closed-face rule) —
               ! a momentum-grid question, not a tracer concentration.
               if (hvel(ktop) <= H_VANISHED .or. hvel(ktop - 1) <= H_VANISHED) beta = 0.0_wp
            end if
            a_diag(i, j, ktop) = -beta
            c_diag(i, j, ktop) = 0.0_wp
            b_diag(i, j, ktop) = 1.0_wp + beta
         end if
         ! ---- Rows ABOVE the first live layer: the identity ----
         ! Empty whenever `k_top = nz`, so bit-identical everywhere a
         ! column has no top-side filler.  Where it is not empty those
         ! rows are inert fillers inside the ice draft: they carry no
         ! mass, they are not coupled to anything (the interior loop's
         ! `zlevel_faces` gates already zero their alpha/beta when the
         ! partial-step knob is on, and this makes the statement hold
         ! with the knob OFF too), and `rhs` still holds `u^n`, so the
         ! Thomas sweep returns them unchanged.  Written explicitly
         ! rather than left to the gates because the surface-BC block
         ! above no longer writes row `nz` when `k_top < nz`.
         do k = ktop + 1, nz
            a_diag(i, j, k) = 0.0_wp
            b_diag(i, j, k) = 1.0_wp
            c_diag(i, j, k) = 0.0_wp
         end do
         ! Ice-shelf top-drag stress BC (Roundabout surface = k=nz; the
         ! MIRROR of the bed row's `dt·λ_bot` fold above).  λ_top is the
         ! Rayleigh RATE the top-drag slot already formed and already
         ! cover- and wet-masked, so an OPEN face contributes an exact
         ! zero here.  A drag is a SINK ⇒ POSITIVE diagonal add ⇒
         ! |amplification| ≤ 1 for ANY h/dt, which is the whole point on
         ! the thin top layers a sigma coordinate leaves near a grounding
         ! line.  At nz = 1 the `if (nz > 1)` block above is skipped and
         ! this row IS the bed row: one layer then carries the bottom
         ! drag, the top drag and the wind, which is exactly right.
         ! Gated INSIDE the single DC (no split loop — NVHPC penalty);
         ! `lambda_top` is guaranteed present when `do_top`.
         if (do_top) then
            b_diag(i, j, ktop) = b_diag(i, j, ktop) + dt*lambda_top(i, j)
         end if
         ! Surface wind-stress Neumann BC (Roundabout surface = k=nz; MIRROR of
         ! MOM6 k=1).  The free-surface flux is the prescribed kinematic
         ! stress τ/ρ₀, a RHS source — the DIAGONAL is unchanged.  The row
         ! is pre-normalized by h_nz (= hf_k here), so the MOM6 surface_stress
         ! term `dt·τ/ρ₀` becomes `dt·(τ/ρ₀)/h_nz_face`.  In the zero-
         ! interior-coupling limit this reduces EXACTLY to the explicit
         ! `du = dt·τ/(ρ₀·h_top)` increment.  `STRESS_H_MIN` matches the
         ! explicit surface-stress kernel's `h_min` so the two paths agree
         ! in the well-resolved limit.  `tau_face` is masked by the face
         ! `min` of the two bounding cells' wet/dry mask (as the explicit
         ! kernel does) so a land face receives no stress.  Gated INSIDE
         ! the single DC.
         ! ICE COVER: under a shelf there is no atmosphere, so the wind
         ! RHS is scaled by `(1 − cover_face)` — a covered face takes the
         ! top DRAG (the diagonal add above) and not the wind.
         !
         ! This is now BELT AND BRACES, and deliberately kept.  P2c masks
         ! the wind at its SOURCE: `ocean_surface_stress_apply_cover`
         ! zeroes the `tau` pair on every face touching a covered cell
         ! (the same OR rule `top_drag_fill_face_cover_impl` uses — one
         ! rule, enforced face-for-face by
         ! `cavity_cover_face_rule_matches_top_drag`), so `tau_face` is
         ! already exactly zero here and this factor multiplies zero by
         ! zero.  It stays because it costs one multiply on a row that is
         ! already being assembled and it keeps the fold correct
         ! STANDALONE — a future forcing path that writes `tau` after the
         ! configure-time mask (a file reader, say) would be caught here
         ! rather than blowing an atmosphere through the ice.  The
         ! asymmetry with the explicit path is GONE: the explicit surface
         ! stress reads the same masked `tau`.
         if (do_stress) then
            wind_open = 1.0_wp
            if (do_top) wind_open = 1.0_wp - cover_face(i, j)
            rhs(i, j, ktop) = rhs(i, j, ktop) + &
                              dt*tau_face(i, j)*wind_open &
                              *min(wet_cell(i_left, j_below), wet_cell(i_right, j_above)) &
                              *inv_rho0/max(hf_k, STRESS_H_MIN)
         end if

         ! ---- Thomas forward sweep ----
         ! remnant-only mode (`solve_momentum = .false.`) runs the c'
         ! recurrence WITHOUT the rhs leg and skips the back-substitution
         ! — u_face is untouched, but c_diag ends up holding exactly the
         ! c' values the remnant solve below needs.
         c_diag(i, j, 1) = c_diag(i, j, 1)/b_diag(i, j, 1)
         if (solve_momentum) then
            rhs(i, j, 1) = rhs(i, j, 1)/b_diag(i, j, 1)
            do k = 2, nz
               denom = b_diag(i, j, k) - a_diag(i, j, k)*c_diag(i, j, k - 1)
               c_diag(i, j, k) = c_diag(i, j, k)/denom
               rhs(i, j, k) = (rhs(i, j, k) - a_diag(i, j, k)*rhs(i, j, k - 1))/denom
            end do

            ! ---- Back-substitution: write directly into u_face ----
            u_face(i, j, nz) = rhs(i, j, nz)
            do k = nz - 1, 1, -1
               u_face(i, j, k) = rhs(i, j, k) - c_diag(i, j, k)*u_face(i, j, k + 1)
            end do
         else
            do k = 2, nz
               denom = b_diag(i, j, k) - a_diag(i, j, k)*c_diag(i, j, k - 1)
               c_diag(i, j, k) = c_diag(i, j, k)/denom
            end do
         end if

         ! ---- Viscous remnant γ_k: second solve, same operator, RHS ≡ 1 ----
         ! γ is DEFINED as the momentum solve's own sensitivity to a uniform
         ! barotropic acceleration, so it is built against the SAME
         ! factorization rather than a freshly-assembled matrix: `a_diag`/
         ! `b_diag` are untouched by the momentum forward sweep above,
         ! and `c_diag` already holds c' (the sweep overwrote it in place),
         ! so `denom = b_diag(k) - a_diag(k)*c_diag(k-1)` recomputes the
         ! IDENTICAL value the momentum sweep used at this k (bit-identical
         ! denominators — no new scratch, no re-derivation from h/kv).
         ! RHS is the unit vector, NOT h_u(k) — Roundabout's rows are already
         ! normalized by h_k (MOM6's un-normalized rows are why its RHS is
         ! h_u(k); dividing MOM6's h_1 remnant RHS by h_1 gives exactly 1).
         if (do_remnant) then
            visc_rem_out(i, j, 1) = 1.0_wp/b_diag(i, j, 1)
            do k = 2, nz
               denom = b_diag(i, j, k) - a_diag(i, j, k)*c_diag(i, j, k - 1)
               visc_rem_out(i, j, k) = (1.0_wp - a_diag(i, j, k)*visc_rem_out(i, j, k - 1))/denom
            end do
            ! Clamp at production (MOM6 clamps at
            ! consumption; the maximum principle bounds γ to (0, 1] so this
            ! is cheap FP insurance, applied as each element is finalized).
            visc_rem_out(i, j, nz) = min(visc_rem_out(i, j, nz), 1.0_wp)
            do k = nz - 1, 1, -1
               visc_rem_out(i, j, k) = visc_rem_out(i, j, k) &
                                       - c_diag(i, j, k)*visc_rem_out(i, j, k + 1)
               visc_rem_out(i, j, k) = min(visc_rem_out(i, j, k), 1.0_wp)
            end do
         end if
      end do
   end subroutine diffuse_velocity_columns_impl