kappa_shear_column_kernel Subroutine

   pure subroutine kappa_shear_column_kernel(grid, this, ms, hT, hS, dt)
      !! Per-column JHL08 solve.  One `do concurrent (j, i)` over owned
      !! cells; every column is solved serially in surface-down order.
      !! Per-thread work is fixed-size `local()` arrays (L1 layout —
      !! all live in GPU registers, no shared memory); the column solve
      !! bodies are the same-module pure `!$acc routine seq` helpers below.
      !!
      !! Pipeline:
      !!   gather (FLIP global bottom-up -> local surface-down): h
      !!          (floored), u,v face-averaged to centre, T,S = hTr/h.
      !!   D4 (massless_merge on + column has a vanished layer): build the
      !!          merge maps, fold vanished layers onto the massive
      !!          sub-grid (nzc<=nz), solve there, interp kappa/TKE back.
      !!          Identity columns (none vanished) bypass the merge.
      !!   precompute: grids, h_Int, boundary length scale, background
      !!          kappa_0 pre-step tridiagonal, frozen EOS buoyancy
      !!          derivatives, initial N^2/S^2.
      !!   outer:  adaptive substepping with predictor-corrector and the
      !!          Picard inner (kappa,Q) solve.
      !!   scatter (FLIP back): kappa_avg -> kd_int, 0 at bed/surface.
      type(hgrid_t), intent(in) :: grid
      type(ocean_kappa_shear_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.
      real(wp), intent(in) :: dt

      integer :: i, j, k, kg, nx, ny, nz, nzc
      real(wp) :: f2_val, hk, inv_h
      logical :: merge_on, do_merge
      ! gathered surface-down column inputs
      real(wp) :: h_sd(NZL), u_sd(NZL), v_sd(NZL), t_sd(NZL), s_sd(NZL)
      ! precomputed interface/layer grids
      real(wp) :: idz_s(NZL), idz_int_s(NZLI), hint_s(NZLI), il2_s(NZLI)
      real(wp) :: kappa_avg_sd(NZLI), tke_avg_sd(NZLI)
      ! D4 massless-merge per-column scratch (+9 arrays: 8 real, 1 int).
      ! kc/kf are the merge maps; hc/uc/vc/tc/sc the merged column;
      ! kappa_c/tke_c the merged-grid interface outputs before interp.
      real(wp) :: hc(NZL), uc(NZL), vc(NZL), tc(NZL), sc(NZL)
      real(wp) :: kappa_c(NZLI), tke_c(NZLI), kf(NZL + 1)
      integer :: kc(NZL + 1)

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      merge_on = this%massless_merge

      do concurrent(j=1:ny, i=1:nx) &
         local(k, kg, nzc, f2_val, hk, inv_h, do_merge, &
               h_sd, u_sd, v_sd, t_sd, s_sd, &
               idz_s, idz_int_s, hint_s, il2_s, &
               kappa_avg_sd, tke_avg_sd, &
               hc, uc, vc, tc, sc, kappa_c, tke_c, kc, kf)

         ! ---- Default outputs (overwritten for wet columns) ----
         do k = 1, nz + 1
            kappa_avg_sd(k) = 0.0_wp
            tke_avg_sd(k) = 0.0_wp
         end do

         if (ms%wet_mask(i, j) > 0.0_wp) then
            ! ---- Gather with the index flip (local k=1 = surface) ----
            ! global layer (kg = nz+1-k) -> local layer k.  Gather RAW
            ! thickness; the I1 precheck below decides floor vs merge.
            do k = 1, nz
               kg = nz + 1 - k
               h_sd(k) = ms%h_layer(i, j, kg)
               u_sd(k) = 0.5_wp*(ms%u_face_x_layer(i, j, kg) + &
                                 ms%u_face_x_layer(i + 1, j, kg))
               v_sd(k) = 0.5_wp*(ms%v_face_y_layer(i, j, kg) + &
                                 ms%v_face_y_layer(i, j + 1, kg))
            end do

            ! ---- I1 precheck (cheap, per column): merge only when the
            ! knob is on AND the column actually carries a vanished layer.
            ! Identity columns bypass the merge machinery entirely, so
            ! knob-on is bit-identical to knob-off on healthy envelopes.
            do_merge = .false.
            if (merge_on) then
               do k = 1, nz
                  if (h_sd(k) < H_VANISHED) then
                     do_merge = .true.
                     exit
                  end if
               end do
            end if

            f2_val = this%f_centre(i, j)*this%f_centre(i, j)

            if (do_merge) then
               ! ---- D4 merge path: raw thickness, T/S back-out with the
               ! div-eps armour (hT = T*h, so hT/h recovers the layer mean
               ! even for a vanished layer).
               do k = 1, nz
                  kg = nz + 1 - k
                  inv_h = 1.0_wp/max(h_sd(k), H_DIV_EPS)
                  t_sd(k) = hT(i, j, kg)*inv_h
                  s_sd(k) = hS(i, j, kg)*inv_h
               end do

               ! Build maps -> merge fields -> solve on nzc -> interp back.
               call massless_build_maps(h_sd, nz, H_VANISHED, nzc, hc, kc, kf)
               call massless_merge_fields(h_sd, kc, nz, nzc, &
                                          u_sd, v_sd, t_sd, s_sd, &
                                          uc, vc, tc, sc)
               call ks_precompute(nzc, this%lz_rescale, hc, &
                                  idz_s, idz_int_s, hint_s, il2_s)
               call ks_solve_column(nzc, dt, f2_val, this%rho0, &
                                    this%ri_crit, this%shearmix_rate, &
                                    this%fri_curvature, this%c_n, this%c_s, &
                                    this%lambda, this%kappa_0, this%kappa_seed, &
                                    this%kappa_trunc, this%tke_bg, this%tol_err, &
                                    this%max_inner_it, this%max_substep_it, &
                                    this%src_max_chg, this%vel_underflow, &
                                    this%eos, &
                                    hc, uc, vc, tc, sc, &
                                    idz_s, idz_int_s, hint_s, il2_s, &
                                    kappa_c, tke_c)
               call massless_interp_back(kappa_c, kc, kf, nz, kappa_avg_sd)
               call massless_interp_back(tke_c, kc, kf, nz, tke_avg_sd)
            else
               ! ---- Existing path (knob off, or knob on + healthy column):
               ! blunt gather floor at H_VANISHED, solve on nz.  Bitwise
               ! identical to pre-D4 behaviour.
               do k = 1, nz
                  kg = nz + 1 - k
                  hk = max(h_sd(k), H_VANISHED)
                  inv_h = 1.0_wp/hk
                  h_sd(k) = hk
                  t_sd(k) = hT(i, j, kg)*inv_h
                  s_sd(k) = hS(i, j, kg)*inv_h
               end do

               ! ---- Precompute the thickness grids (h_Int, 1/h, L_bdry) ----
               call ks_precompute(nz, this%lz_rescale, h_sd, &
                                  idz_s, idz_int_s, hint_s, il2_s)

               ! ---- Adaptive outer solve (background pre-step + EOS
               !      buoyancy derivatives are done inside) ----
               call ks_solve_column(nz, dt, f2_val, this%rho0, &
                                    this%ri_crit, this%shearmix_rate, &
                                    this%fri_curvature, this%c_n, this%c_s, &
                                    this%lambda, this%kappa_0, this%kappa_seed, &
                                    this%kappa_trunc, this%tke_bg, this%tol_err, &
                                    this%max_inner_it, this%max_substep_it, &
                                    this%src_max_chg, this%vel_underflow, &
                                    this%eos, &
                                    h_sd, u_sd, v_sd, t_sd, s_sd, &
                                    idz_s, idz_int_s, hint_s, il2_s, &
                                    kappa_avg_sd, tke_avg_sd)
            end if
         end if

         ! ---- Scatter with the flip; force exact 0 at bed + surface ----
         ! local interface K -> global interface (nz+2-K).
         do k = 1, nz + 1
            kg = nz + 2 - k
            this%kd_int(i, j, kg) = kappa_avg_sd(k)
            this%tke_int(i, j, kg) = tke_avg_sd(k)
         end do
         this%kd_int(i, j, 1) = 0.0_wp
         this%kd_int(i, j, nz + 1) = 0.0_wp
         this%tke_int(i, j, 1) = 0.0_wp
         this%tke_int(i, j, nz + 1) = 0.0_wp
      end do
   end subroutine kappa_shear_column_kernel