build_factorize_tracer_matrix Subroutine

   pure subroutine build_factorize_tracer_matrix(nx, ny, nz, use_harmonic, dt, kv_centre, &
                                                 h_layer, a_diag, b_diag, c_diag)
      !! Build the backward-Euler tridiagonal per cell column and run the
      !! Thomas forward factorization — the tracer-INDEPENDENT half of the
      !! vertical-diffusion solve (depends only on kv/h/dt).  Run once per
      !! stage; every registered tracer then reuses the factored
      !! coefficients via `apply_factored_tracer`.
      !!
      !! `kv_centre(:, :, k)` is the diffusivity at the BOTTOM interface
      !! of layer `k`; `kv_centre(:, :, nz+1)` is the surface interface.
      !! On exit:
      !!   a_diag = sub-diagonal (unchanged), for the per-tracer RHS sweep
      !!   b_diag = the Thomas pivots (b(1) = raw diagonal, b(k>1) = denom)
      !!   c_diag = super-diagonal already divided by its pivot
      integer, intent(in) :: nx, ny, nz
      logical, intent(in) :: use_harmonic
      real(wp), intent(in) :: dt
      real(wp), intent(in) :: kv_centre(nx, ny, nz + 1)
      real(wp), intent(in) :: h_layer(nx, ny, nz)
      real(wp), intent(inout) :: a_diag(nx, ny, nz)
      real(wp), intent(inout) :: b_diag(nx, ny, nz)
      real(wp), intent(inout) :: c_diag(nx, ny, nz)

      integer :: i, j, k
      real(wp) :: dz_face, alpha, beta, denom
      real(wp) :: hc, hm, hp

      ! Vanishing-layer handling (load-bearing for conservation under
      ! Fox-Kemper × windowed tracer advect): a thin z* surface layer can be
      ! driven to h ≤ 0 on an intermediate RK2 stage by the combined FK +
      ! resolved transport (Lagrangian, before the ALE remap).  The raw
      ! 1/h in α/β would then go negative/Inf and corrupt the WHOLE column's
      ! solve, dropping tracer mass.  We floor the per-layer thickness to
      ! H_VANISHED in the denominators (the SAME h̃ apply_factored_tracer
      ! uses) and ZERO the diffusive flux at any interface touching a
      ! vanishing layer (h ≤ H_VANISHED) — so a collapsed layer is decoupled
      ! (identity row) and its frozen tracer mass is preserved exactly,
      ! while its thick neighbours conserve among themselves.  For
      ! h ≫ H_VANISHED (sigma / double_gyre) h̃ = h and every interface is
      ! active ⇒ bit-identical.
      do concurrent(j=1:ny, i=1:nx) local(dz_face, alpha, beta, denom, k, hc, hm, hp)
         ! k = 1: bed BC (no flux below).  α uses the interface above
         ! layer 1 (= kv_centre(:, :, 2)).
         hc = max(h_layer(i, j, 1), H_VANISHED)
         ! SINGLE-LAYER COLUMN (nz = 1): the bed row IS the surface row and
         ! there is no interior interface, so alpha is identically zero.
         ! Without the gate this reads `h_layer(i, j, 2)` — past the end of
         ! a (nx, ny, 1) array — and the k = nz block below then reads
         ! `h_layer(i, j, 0)` and overwrites this row.  Same defect (and
         ! same shape of fix) as the velocity tridiagonal below.
         ! Loop-invariant gate INSIDE the single DC ⇒ nz >= 2 bit-identical.
         alpha = 0.0_wp
         if (nz > 1) then
            hp = max(h_layer(i, j, 2), H_VANISHED)
            dz_face = face_thick(hc, hp, use_harmonic)
            alpha = dt*kv_centre(i, j, 2)/(hc*dz_face)
            if (h_layer(i, j, 1) <= H_VANISHED .or. h_layer(i, j, 2) <= H_VANISHED) alpha = 0.0_wp
         end if
         a_diag(i, j, 1) = 0.0_wp
         c_diag(i, j, 1) = -alpha
         b_diag(i, j, 1) = 1.0_wp + alpha

         ! k = 2..nz-1: interior.  β = kv_centre(k), α = kv_centre(k+1).
         do k = 2, nz - 1
            hm = max(h_layer(i, j, k - 1), H_VANISHED)
            hc = max(h_layer(i, j, k), H_VANISHED)
            hp = max(h_layer(i, j, k + 1), H_VANISHED)
            beta = dt*kv_centre(i, j, k)/(hc*face_thick(hm, hc, use_harmonic))
            alpha = dt*kv_centre(i, j, k + 1)/(hc*face_thick(hc, hp, use_harmonic))
            if (h_layer(i, j, k) <= H_VANISHED .or. h_layer(i, j, k - 1) <= H_VANISHED) beta = 0.0_wp
            if (h_layer(i, j, k) <= H_VANISHED .or. h_layer(i, j, k + 1) <= H_VANISHED) alpha = 0.0_wp
            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 = nz: surface BC (no flux above).  β uses kv_centre(nz).
         ! SINGLE-LAYER COLUMN (nz = 1): already built as the bed row above;
         ! skip (see the k = 1 gate).  Otherwise this reads h_layer(:, :, 0).
         if (nz > 1) then
            hm = max(h_layer(i, j, nz - 1), H_VANISHED)
            hc = max(h_layer(i, j, nz), H_VANISHED)
            dz_face = face_thick(hm, hc, use_harmonic)
            beta = dt*kv_centre(i, j, nz)/(hc*dz_face)
            if (h_layer(i, j, nz) <= H_VANISHED .or. h_layer(i, j, nz - 1) <= H_VANISHED) beta = 0.0_wp
            a_diag(i, j, nz) = -beta
            c_diag(i, j, nz) = 0.0_wp
            b_diag(i, j, nz) = 1.0_wp + beta
         end if

         ! ---- Thomas forward factorization (matrix only) ----
         ! Store the pivot `denom` back into b_diag and the eliminated
         ! super-diagonal `c/denom` into c_diag.  b_diag(1) keeps the raw
         ! diagonal (its own pivot); a_diag stays the raw sub-diagonal.
         c_diag(i, j, 1) = c_diag(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
            b_diag(i, j, k) = denom
         end do
      end do
   end subroutine build_factorize_tracer_matrix