barotropic_substep_linear Subroutine

   pure subroutine barotropic_substep_linear(grid, metrics, bt_work, force_u, force_v, n_steps, dt_inner, &
                                             eta_forcing)
      !! Public only for the unit-test suite (no production module imports it);
      !! ignore when developing production code in other modules.
      !! Linearized barotropic-substep kernel.  Forward-Euler steps
      !! the barotropic state (eta, ubt, vbt) at dt_inner for n_steps
      !! against the constant slow forcing `force_u`, `force_v` (each
      !! at the same C-grid location as ubt, vbt).  Accumulates the
      !! per-step running sums into `ubt_sum`, `vbt_sum`, `eta_sum`;
      !! at the end divides by `n_steps` and stores the time-mean
      !! back into `bt_eta`, `bt_ubt`, `bt_vbt` for the caller.
      !!
      !! Dynamics (closed walls, no Coriolis — linearized gravity-
      !! wave only for unit tests):
      !!
      !!   ∂η/∂t   = -∂(H_ref · u_bt)/∂x - ∂(H_ref · v_bt)/∂y
      !!   ∂u_bt/∂t = -g · ∂η/∂x + force_u
      !!   ∂v_bt/∂t = -g · ∂η/∂y + force_v
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
      type(barotropic_workstate_t), intent(inout) :: bt_work
      ! Explicit-shape (not assumed-shape): a `(:, :)` dummy carries an array
      ! descriptor that must be device-resident inside the `do concurrent`
      ! kernel.  stdpar manages that, but the OpenMP-target variant reads a
      ! non-mapped descriptor → CUDA illegal/misaligned access in the fast
      ! loop.  Explicit shape passes base + dims (no descriptor) and also
      ! avoids the per-launch descriptor-walk memcpys.
      real(wp), intent(in) :: force_u(grid%nx_total + 1, grid%ny_total)
      real(wp), intent(in) :: force_v(grid%nx_total, grid%ny_total + 1)
      integer, intent(in) :: n_steps
      real(wp), intent(in) :: dt_inner
      real(wp), intent(in), optional :: eta_forcing(grid%nx_total, grid%ny_total)
         !! Optional equilibrium-tide elevation (m); when present the PGF
         !! drives grad(eta - eta_forcing).  Absent => bit-identical.

      integer :: i, j, n, nx, ny
      logical :: tide_on
      real(wp) :: inv_n
      real(wp) :: d_eta
      real(wp) :: h_face_E, h_face_W, h_face_N, h_face_S
      real(wp) :: flux_x_R, flux_x_L, flux_y_N, flux_y_S, div_h_u
      real(wp) :: G   !! η-gradient PGF coefficient, sourced from
                      !! `bt_work%g_bt`.  Defaults to 9.81 (full gravity)
                      !! but the gprime path can lower it to g_FS.

      nx = grid%nx_total
      ny = grid%ny_total
      G = bt_work%g_bt
      tide_on = present(eta_forcing)

      ! Reset running sums.
      do concurrent(j=1:ny, i=1:nx)
         bt_work%eta_sum(i, j) = 0.0_wp
      end do
      do concurrent(j=1:ny, i=1:nx + 1)
         bt_work%ubt_sum(i, j) = 0.0_wp
      end do
      do concurrent(j=1:ny + 1, i=1:nx)
         bt_work%vbt_sum(i, j) = 0.0_wp
      end do

      do n = 1, n_steps
         ! ---- Pass 1: eta update at every cell ----
         do concurrent(j=1:ny, i=1:nx) &
            local(h_face_E, h_face_W, h_face_N, h_face_S, &
                  flux_x_R, flux_x_L, flux_y_N, flux_y_S, div_h_u)
            if (i < nx) then
               h_face_E = 0.5_wp*(bt_work%bt_H_ref(i, j) + bt_work%bt_H_ref(i + 1, j))
            else
               h_face_E = bt_work%bt_H_ref(i, j)
            end if
            if (i > 1) then
               h_face_W = 0.5_wp*(bt_work%bt_H_ref(i - 1, j) + bt_work%bt_H_ref(i, j))
            else
               h_face_W = bt_work%bt_H_ref(i, j)
            end if
            if (j < ny) then
               h_face_N = 0.5_wp*(bt_work%bt_H_ref(i, j) + bt_work%bt_H_ref(i, j + 1))
            else
               h_face_N = bt_work%bt_H_ref(i, j)
            end if
            if (j > 1) then
               h_face_S = 0.5_wp*(bt_work%bt_H_ref(i, j - 1) + bt_work%bt_H_ref(i, j))
            else
               h_face_S = bt_work%bt_H_ref(i, j)
            end if
            if (bt_work%use_bt_cont_type) then
               ! MOM6 piecewise-cubic flux closure: transport is capped
               ! by the upstream column's per-layer h sum so the BT mode
               ! can't pump mass through a face where the upstream
               ! column lacks the height to supply it.  BTCL_u/v were
               ! built from the slow ML snapshot in
               ! `set_local_BT_cont_types`.
               flux_x_R = find_uhbt(bt_work%bt_ubt(i + 1, j), bt_work%BTCL_u(i + 1, j))*metrics%dy_cu_bt(i + 1, j)
               flux_x_L = find_uhbt(bt_work%bt_ubt(i, j), bt_work%BTCL_u(i, j))*metrics%dy_cu_bt(i, j)
               flux_y_N = find_vhbt(bt_work%bt_vbt(i, j + 1), bt_work%BTCL_v(i, j + 1))*metrics%dx_cv_bt(i, j + 1)
               flux_y_S = find_vhbt(bt_work%bt_vbt(i, j), bt_work%BTCL_v(i, j))*metrics%dx_cv_bt(i, j)
            else if (bt_work%use_upstream_h_face) then
               ! Upstream-PPM h_face from the slow ML snapshot.  Same
               ! convention slow continuity uses, so the BT mode's
               ! mass flux and the per-layer mass flux carry the
               ! same face thickness — eliminates the centred-vs-
               ! upstream mismatch that leaves phantom bed-layer
               ! velocity at slopes.  Built once per outer step in
               ! `compute_h_face_upstream`.  bt_H_ref unused on
               ! this branch — the upstream column sum carries the
               ! total thickness (incl. η at top of stage).
               flux_x_R = bt_work%h_face_up_x(i + 1, j)*bt_work%bt_ubt(i + 1, j)*metrics%dy_cu_bt(i + 1, j)
               flux_x_L = bt_work%h_face_up_x(i, j)*bt_work%bt_ubt(i, j)*metrics%dy_cu_bt(i, j)
               flux_y_N = bt_work%h_face_up_y(i, j + 1)*bt_work%bt_vbt(i, j + 1)*metrics%dx_cv_bt(i, j + 1)
               flux_y_S = bt_work%h_face_up_y(i, j)*bt_work%bt_vbt(i, j)*metrics%dx_cv_bt(i, j)
            else
               flux_x_R = h_face_E*bt_work%bt_ubt(i + 1, j)*metrics%dy_cu_bt(i + 1, j)
               flux_x_L = h_face_W*bt_work%bt_ubt(i, j)*metrics%dy_cu_bt(i, j)
               flux_y_N = h_face_N*bt_work%bt_vbt(i, j + 1)*metrics%dx_cv_bt(i, j + 1)
               flux_y_S = h_face_S*bt_work%bt_vbt(i, j)*metrics%dx_cv_bt(i, j)
            end if
            ! Conservative transport divergence · iareaT (= inv_dx/inv_dy on uniform).
            div_h_u = ((flux_x_R - flux_x_L) + (flux_y_N - flux_y_S))*metrics%iareaT(i, j)
            bt_work%bt_eta(i, j) = bt_work%bt_eta(i, j) - dt_inner*div_h_u
         end do

         ! ---- Pass 2: ubt update at interior east faces ----
         ! Single-rank reference kernel; physical-edge gating lives in
         ! the nonlinear production variant (barotropic_substep_nonlinear).
         do concurrent(j=1:ny, i=2:nx) local(d_eta)
            d_eta = bt_work%bt_eta(i, j) - bt_work%bt_eta(i - 1, j)
            if (tide_on) d_eta = d_eta - (eta_forcing(i, j) - eta_forcing(i - 1, j))
            bt_work%bt_ubt(i, j) = bt_work%bt_ubt(i, j) + dt_inner*( &
                                   -G*d_eta*metrics%idxCu(i, j) + &
                                   force_u(i, j))
         end do
         do concurrent(j=1:ny)
            bt_work%bt_ubt(1, j) = 0.0_wp
            bt_work%bt_ubt(nx + 1, j) = 0.0_wp
            ! Physical-wall closure — slow continuity blocks flow at
            ! these faces, so the substep loop must too.  See header.
            bt_work%bt_ubt(grid%nghost + 1, j) = 0.0_wp
            bt_work%bt_ubt(grid%nghost + grid%nx_phys + 1, j) = 0.0_wp
         end do

         ! ---- Pass 3: vbt update at interior north faces ----
         do concurrent(j=2:ny, i=1:nx) local(d_eta)
            d_eta = bt_work%bt_eta(i, j) - bt_work%bt_eta(i, j - 1)
            if (tide_on) d_eta = d_eta - (eta_forcing(i, j) - eta_forcing(i, j - 1))
            bt_work%bt_vbt(i, j) = bt_work%bt_vbt(i, j) + dt_inner*( &
                                   -G*d_eta*metrics%idyCv(i, j) + &
                                   force_v(i, j))
         end do
         do concurrent(i=1:nx)
            bt_work%bt_vbt(i, 1) = 0.0_wp
            bt_work%bt_vbt(i, ny + 1) = 0.0_wp
            bt_work%bt_vbt(i, grid%nghost + 1) = 0.0_wp
            bt_work%bt_vbt(i, grid%nghost + grid%ny_phys + 1) = 0.0_wp
         end do

         ! ---- Accumulate running sums ----
         do concurrent(j=1:ny, i=1:nx)
            bt_work%eta_sum(i, j) = bt_work%eta_sum(i, j) + bt_work%bt_eta(i, j)
         end do
         do concurrent(j=1:ny, i=1:nx + 1)
            bt_work%ubt_sum(i, j) = bt_work%ubt_sum(i, j) + bt_work%bt_ubt(i, j)
         end do
         do concurrent(j=1:ny + 1, i=1:nx)
            bt_work%vbt_sum(i, j) = bt_work%vbt_sum(i, j) + bt_work%bt_vbt(i, j)
         end do
      end do

      ! Snapshot end-of-loop η/u/v BEFORE the time-mean overwrite.
      ! `apply_bt_correction` uses all three so the recombined
      ! per-layer momentum and the h_layer rescale both live at
      ! t + dt_outer (Hallberg 2009).  Mixing end-step velocity with
      ! time-mean SSH puts the two legs half a substep out of phase
      ! and corrupts the gravity-wave dispersion.
      do concurrent(j=1:ny, i=1:nx)
         bt_work%bt_eta_end(i, j) = bt_work%bt_eta(i, j)
      end do
      do concurrent(j=1:ny, i=1:nx + 1)
         bt_work%bt_ubt_end(i, j) = bt_work%bt_ubt(i, j)
      end do
      do concurrent(j=1:ny + 1, i=1:nx)
         bt_work%bt_vbt_end(i, j) = bt_work%bt_vbt(i, j)
      end do

      ! Time-mean: divide running sum by n_steps, store back into bt_*.
      inv_n = 1.0_wp/real(n_steps, wp)
      do concurrent(j=1:ny, i=1:nx)
         bt_work%bt_eta(i, j) = bt_work%eta_sum(i, j)*inv_n
      end do
      do concurrent(j=1:ny, i=1:nx + 1)
         bt_work%bt_ubt(i, j) = bt_work%ubt_sum(i, j)*inv_n
      end do
      do concurrent(j=1:ny + 1, i=1:nx)
         bt_work%bt_vbt(i, j) = bt_work%vbt_sum(i, j)*inv_n
      end do
   end subroutine barotropic_substep_linear