ocean_apply_ale_remap_centres Subroutine

   pure subroutine ocean_apply_ale_remap_centres(grid, vcoord, ms, bt_eta, bt_H_ref, method, eos, dt)
      !! Orchestrate the centre-cell pass of the ALE remap step (h_layer +
      !! tracers). Public only for the unit-test suite.
      !! Sequence: skip if EULERIAN_Z/LAGRANGIAN; snapshot column total + h_old;
      !! build `vcoord%target_h`; remap each tracer h_old→target_h via the PPM
      !! column kernel; set h_layer = target_h; recompute
      !! bt_eta = sum_k(h_layer) - bt_H_ref. Face velocities remapped separately.
      !! Takes bt_eta/bt_H_ref directly (not ocean_dyn_t) so this module sits
      !! below the split driver in the dependency tree.
      type(hgrid_t), intent(in) :: grid
      type(ocean_vcoord_t), intent(inout) :: vcoord
      type(multilayer_state_t), intent(inout) :: ms
      ! assumed-shape-ok: outer-driver allocatable; grid%nx_total available; deferred.
      real(wp), intent(inout) :: bt_eta(:, :)
         !! Free-surface anomaly η at cell centres (m).  Updated to
         !! `sum_k(h_layer) - bt_H_ref` on exit (step 7).
      real(wp), intent(in) :: bt_H_ref(:, :)  ! assumed-shape-ok: same as bt_eta above
         !! Reference column depth H (m, positive-down).  Constant for
         !! the ocean path; passed in so the routine doesn't need a
         !! handle on the split-driver state.
      type(eos_t), intent(in), optional :: eos
         !! Device-resident EOS handle — required only for `VCOORD_RHO`
         !! (isopycnal density inversion); ignored by the geometric coords.
      real(wp), intent(in), optional :: dt
         !! Outer/thermo timestep (s) for the grid time-filter.  Absent or
         !! `regrid_time_scale = 0` (default) ⇒ filter skipped ⇒
         !! bit-identical.  See `ocean_apply_ale_remap_step`.
      integer, intent(in), optional :: method
         !! REMAP_PCM / REMAP_PLM / REMAP_PPM. Defaults to PPM (parabolic stencil
         !! cuts the spurious vertical mixing a 1st-order limiter introduces).

      integer :: t, m, nx, ny, nz, i, j, k
      real(wp) :: wtd
      logical :: do_tfilter

      ! Eulerian-z and Lagrangian/isopycnal both skip remap: the former
      ! holds h at H·dsig via vert-advection cancellation, the latter
      ! lets h evolve freely (target = current h).
      if (vcoord%coord_type == VCOORD_EULERIAN_Z .or. &
          vcoord%coord_type == VCOORD_LAGRANGIAN) return
      if (.not. vcoord%is_init) return
      if (.not. allocated(ms%h_layer)) return

      m = REMAP_PPM
      if (present(method)) m = method
      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml

      ! --- 2. Current column total (persistent device-mapped scratch on vcoord)
      do concurrent(j=1:ny, i=1:nx) local(k)
         vcoord%remap_total_h(i, j) = 0.0_wp
         do k = 1, nz
            vcoord%remap_total_h(i, j) = vcoord%remap_total_h(i, j) &
                                         + ms%h_layer(i, j, k)
         end do
      end do

      ! --- 3. Snapshot h_old on-device (host source= reads stale memory after a
      ! dynamics-side OpenACC kernel). Before target_h so VCOORD_RHO can read it.
      do concurrent(k=1:nz, j=1:ny, i=1:nx)
         vcoord%remap_h_old(i, j, k) = ms%h_layer(i, j, k)
      end do

      ! --- 4. Populate target_h (geometric vs isopycnal dispatch) ---
      do concurrent(j=1:ny, i=1:nx)
         vcoord%remap_h_ref(i, j) = vcoord%remap_total_h(i, j) - bt_eta(i, j)
      end do
      if ((vcoord%coord_type == VCOORD_RHO .or. vcoord%coord_type == VCOORD_HYCOM) &
          .and. present(eos) .and. allocated(ms%tracers) &
          .and. ms%idx_temperature > 0 .and. ms%idx_salinity > 0) then
         call build_ts_concentration(nx, ny, nz, vcoord%remap_h_old, &
                                     ms%tracers(ms%idx_temperature)%hTr, &
                                     ms%tracers(ms%idx_salinity)%hTr, &
                                     vcoord%remap_conc_t, vcoord%remap_conc_s)
         ! HYCOM = RHO inversion + z*-floor/monotonize deltas (hybrid=.true.);
         ! pure RHO passes hybrid=.false. (bit-identical to RHO-only kernel).
         call vcoord%compute_target_h_rho(vcoord%remap_h_ref, bt_eta, &
                                          vcoord%remap_conc_t, vcoord%remap_conc_s, eos, &
                                          hybrid=(vcoord%coord_type == VCOORD_HYCOM))
      else
         call vcoord%compute_target_h(vcoord%remap_h_ref, bt_eta)
      end if

      ! --- 4b. Grid time-filter (White & Adcroft 2008): relax target toward it
      ! from the old grid by wtd = dt/(τ+dt). Convex blend ⇒ column total
      ! conserved. τ=0 (default) or dt absent ⇒ skipped, bit-identical.
      do_tfilter = vcoord%regrid_time_scale > 0.0_wp .and. present(dt)
      if (do_tfilter) then
         wtd = dt/(vcoord%regrid_time_scale + dt)
         do concurrent(k=1:nz, j=1:ny, i=1:nx)
            vcoord%target_h(i, j, k) = vcoord%remap_h_old(i, j, k) &
                                       + wtd*(vcoord%target_h(i, j, k) - vcoord%remap_h_old(i, j, k))
         end do
      end if

      ! --- 5. Remap every tracer ---
      if (allocated(ms%tracers)) then
         do t = 1, size(ms%tracers)
            if (.not. allocated(ms%tracers(t)%hTr)) cycle
            select case (ms%tracers(t)%budget_id)
            case (TRACER_BUDGET_HEAT)
               call ocean_remap_tracer_field( &
                  nx, ny, nz, vcoord%remap_h_old, vcoord%target_h, ms%tracers(t)%hTr, m, &
                  vcoord%remap_boundary_extrap, vcoord%remap_nonuniform_weights, &
                  budget=ms%heat_budget_remap)
            case (TRACER_BUDGET_SALT)
               call ocean_remap_tracer_field( &
                  nx, ny, nz, vcoord%remap_h_old, vcoord%target_h, ms%tracers(t)%hTr, m, &
                  vcoord%remap_boundary_extrap, vcoord%remap_nonuniform_weights, &
                  budget=ms%salt_budget_remap)
            case default
               call ocean_remap_tracer_field( &
                  nx, ny, nz, vcoord%remap_h_old, vcoord%target_h, ms%tracers(t)%hTr, m, &
                  vcoord%remap_boundary_extrap, vcoord%remap_nonuniform_weights)
            end select
         end do
      end if

      ! --- 6. h_layer = target_h; capture mass-budget delta ---
      do concurrent(k=1:nz, j=1:ny, i=1:nx)
         ms%mass_budget_remap(i, j, k) = ms%mass_budget_remap(i, j, k) &
                                         + (vcoord%target_h(i, j, k) - vcoord%remap_h_old(i, j, k))
         ms%h_layer(i, j, k) = vcoord%target_h(i, j, k)
      end do

      ! --- 7. Recompute bt_eta from new sum (round-off-tight) ---
      do concurrent(j=1:ny, i=1:nx) local(k)
         bt_eta(i, j) = -bt_H_ref(i, j)
         do k = 1, nz
            bt_eta(i, j) = bt_eta(i, j) + ms%h_layer(i, j, k)
         end do
      end do
   end subroutine ocean_apply_ale_remap_centres