pure subroutine ocean_apply_ale_remap_step(grid, vcoord, ms, bt_eta, bt_H_ref, method, eos, dt)
!! Top-level entry the driver calls between outer steps: snapshot h_old once,
!! remap centres (h_layer + tracers) AND faces using the same snapshot, then
!! re-derive bt_eta. Returns early for EULERIAN_Z/LAGRANGIAN.
!! `eos` (optional): required ONLY for VCOORD_RHO (isopycnal density inversion).
!! `dt` (optional, s): only for the grid time-filter (regrid_time_scale > 0);
!! absent or τ=0 (default) ⇒ filter skipped, bit-identical.
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(:, :)
real(wp), intent(in) :: bt_H_ref(:, :) ! assumed-shape-ok: outer-driver allocatable; grid%nx_total available; deferred
integer, intent(in), optional :: method
type(eos_t), intent(in), optional :: eos
real(wp), intent(in), optional :: dt
integer :: t, m, nx, ny, nz, i, j, k
real(wp) :: wtd
logical :: do_tfilter, conserve_ke, bnd_extrap, nonunif
! 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
bnd_extrap = vcoord%remap_boundary_extrap
nonunif = vcoord%remap_nonuniform_weights
nx = grid%nx_total
ny = grid%ny_total
nz = ms%nz_ml
! 1. Column total + target_h (persistent scratch on vcoord — per-call
! allocates here generate H↔D transfers, not being in the present table).
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
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
! 2. Snapshot h_old on-device (host source= reads stale memory after a
! dynamics-side OpenACC kernel). Before target_h for VCOORD_RHO's inversion.
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
! 3. Populate target_h. Geometric coords use compute_target_h; isopycnal
! needs per-layer T/S concentrations + EOS via compute_target_h_rho.
! Concentrations built into persistent scratch (guarded c = hTr/h).
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
! 3b. 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
! `vcoord%remap_h_old` and `vcoord%target_h` are the exact pair every
! column kernel below consumes, and neither is overwritten by steps 4-7
! — so `ocean_remap_scan_preconditions` can assert them from the
! (impure) driver AFTER this call. See `&vcoord_nml
! remap_check_preconditions` and `rdb_ocean_dyn :: ocean_dyn_step_split`.
! 4. Tracer remap (centre cells)
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, &
bnd_extrap, nonunif, 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, &
bnd_extrap, nonunif, 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, &
bnd_extrap, nonunif)
end select
end do
end if
! 5. Face-velocity remap (h_old → target_h). KE-conserving anomaly rescale
! gated on the vcoord knob (default off ⇒ momentum-only, bit-identical).
conserve_ke = vcoord%remap_vel_conserve_ke
if (allocated(ms%u_face_x_layer) .and. allocated(ms%v_face_y_layer)) then
call remap_x_face_velocity(nx, ny, nz, vcoord%remap_h_old, vcoord%target_h, &
ms%u_face_x_layer, m, conserve_ke, &
vcoord%zfixed_closed_faces, bnd_extrap, nonunif)
call remap_y_face_velocity(nx, ny, nz, vcoord%remap_h_old, vcoord%target_h, &
ms%v_face_y_layer, m, conserve_ke, &
vcoord%zfixed_closed_faces, bnd_extrap, nonunif)
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. Re-derive bt_eta
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_step