ocean_slopes_vert_fill_ts Subroutine

   subroutine ocean_slopes_vert_fill_ts(nx, ny, nz, h_layer, t_htr, s_htr, &
                                        kd_smooth, dt, t_fill, s_fill)
      !! Fill massless layers in T/S with sensible values via one pass of
      !! constant-`kappa·dt` vertical diffusion — a SINGLE forward-elim +
      !! back-sub Thomas sweep per column (no iteration).  Operates on the
      !! tracer-from-hTr conversion (`T = hTr/h`) and writes the scratch
      !! `t_fill`/`s_fill`; the prognostic tracers are untouched.
      !!
      !! `kap_dt_x2 = 2·kappa·dt`; the inter-layer entrainment is
      !! `ent(K) = kap_dt_x2 / ((h(k)+h(k+1)) + h_neglect)`.  Surface +
      !! bed boundary rows close the tridiagonal exactly.  Column locals
      !! are fixed-size (`NZ_STACK_MAX`) so the `local()` clause is legal
      !! on `-stdpar=gpu` (dummy-sized automatics crash).
      integer, intent(in) :: nx, ny, nz
      real(wp), intent(in) :: h_layer(nx, ny, nz)
      real(wp), intent(in) :: t_htr(nx, ny, nz)
      real(wp), intent(in) :: s_htr(nx, ny, nz)
      real(wp), intent(in) :: kd_smooth, dt
      real(wp), intent(out) :: t_fill(nx, ny, nz)
      real(wp), intent(out) :: s_fill(nx, ny, nz)

      integer :: i, j, k
      real(wp) :: kap_dt_x2, h_neglect, h0c
      real(wp) :: ent(NZ_STACK_MAX + 1)
      real(wp) :: c1(NZ_STACK_MAX)
      real(wp) :: b1, d1, h_tr, h_eff
      real(wp) :: t_in(NZ_STACK_MAX), s_in(NZ_STACK_MAX)

      kap_dt_x2 = 2.0_wp*kd_smooth*dt
      h_neglect = H_DIV_EPS
      h0c = h_neglect

      if (kap_dt_x2 <= 0.0_wp .or. nz < 2) then
         ! No smoothing — pass the raw tracer-from-hTr through.
         do concurrent(k=1:nz, j=1:ny, i=1:nx)
            h_eff = max(h_layer(i, j, k), H_VANISHED)
            t_fill(i, j, k) = t_htr(i, j, k)/h_eff
            s_fill(i, j, k) = s_htr(i, j, k)/h_eff
         end do
         return
      end if

      ! Per-column Thomas sweep.  k=1 bed ... k=nz surface (bottom-up);
      ! the tridiagonal couples layer k to k±1 identically regardless of
      ! orientation, so the bottom-up index runs the published sweep with
      ! "k=1" as the first boundary row.
      do concurrent(j=1:ny, i=1:nx) &
         local(k, ent, c1, b1, d1, h_tr, h_eff, t_in, s_in)
         ! T,S from hTr/h (floor h consistently — vanished layers feed a
         ! near-zero raw T that the diffusion then overwrites).
         do k = 1, nz
            h_eff = max(h_layer(i, j, k), H_VANISHED)
            t_in(k) = t_htr(i, j, k)/h_eff
            s_in(k) = s_htr(i, j, k)/h_eff
         end do

         ! Forward elimination — first (bed) boundary row at k=1.
         ent(2) = kap_dt_x2/((h_layer(i, j, 1) + h_layer(i, j, 2)) + h0c)
         h_tr = h_layer(i, j, 1) + h_neglect
         b1 = 1.0_wp/(h_tr + ent(2))
         d1 = b1*h_tr
         t_fill(i, j, 1) = (b1*h_tr)*t_in(1)
         s_fill(i, j, 1) = (b1*h_tr)*s_in(1)
         do k = 2, nz - 1
            ent(k + 1) = kap_dt_x2/((h_layer(i, j, k) + h_layer(i, j, k + 1)) + h0c)
            h_tr = h_layer(i, j, k) + h_neglect
            c1(k) = ent(k)*b1
            b1 = 1.0_wp/((h_tr + d1*ent(k)) + ent(k + 1))
            d1 = b1*(h_tr + d1*ent(k))
            t_fill(i, j, k) = b1*(h_tr*t_in(k) + ent(k)*t_fill(i, j, k - 1))
            s_fill(i, j, k) = b1*(h_tr*s_in(k) + ent(k)*s_fill(i, j, k - 1))
         end do
         ! Last (surface) boundary row at k=nz.
         c1(nz) = ent(nz)*b1
         h_tr = h_layer(i, j, nz) + h_neglect
         b1 = 1.0_wp/(h_tr + d1*ent(nz))
         t_fill(i, j, nz) = b1*(h_tr*t_in(nz) + ent(nz)*t_fill(i, j, nz - 1))
         s_fill(i, j, nz) = b1*(h_tr*s_in(nz) + ent(nz)*s_fill(i, j, nz - 1))
         ! Back substitution.
         do k = nz - 1, 1, -1
            t_fill(i, j, k) = t_fill(i, j, k) + c1(k + 1)*t_fill(i, j, k + 1)
            s_fill(i, j, k) = s_fill(i, j, k) + c1(k + 1)*s_fill(i, j, k + 1)
         end do
      end do
   end subroutine ocean_slopes_vert_fill_ts