continuity_tracer_drain Subroutine

   subroutine continuity_tracer_drain(grid, metrics, this, ms, ratio, bc)
      !! Phase-2 (6b) windowed horizontal tracer-advection drain.
      !!
      !! Spends the accumulated face transports `this%uhtr/vhtr` (built
      !! over `ratio` outer steps in `TR_MODE_ACCUMULATE`) onto the
      !! frozen tracer mass `ms%tracers(:)%hTr` using a fixed-budget
      !! Colella-Woodward swept-average PPM sub-cycle (Adcroft &
      !! Hallberg 2006; MOM6 ADVECT_PPM).  Conservation is structural:
      !!
      !!   hprev = areaT·h_end + div(uhtr)        (≈ window-start thickness)
      !!   Tr_start = hTr_frozen / hprev          (window-start concentration)
      !!
      !! the drain evolves hprev → h_end while moving tracer with each
      !! limited transport portion; at the end hprev == h_end so
      !! Σ(hTr) = Σ(Tr_start·hprev) = Σ(hTr_frozen) to round-off.  Do NOT
      !! seed the concentration from hTr/h_end (that is the drifted value).
      !!
      !! Per-pass re-reconstruction (V2): the PPM parabola is rebuilt from
      !! the CURRENT Tr at the start of every sub-cycle pass.  The
      !! `hup/hlos/min_h` two-test limiter sets the drained transport
      !! `uhh` per pass; the swept-average parabola sets the concentration
      !! multiplying it — the two are orthogonal.
      !!
      !! Single-rank scope: periodic-x / north-fold seams are handled via
      !! the existing ocean periodic/fold ghost wraps.  Multi-rank C-grid
      !! halos ride E1; gated behind the TODO(E1) note below.
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
      type(continuity_t), intent(inout) :: this
      type(multilayer_state_t), intent(inout) :: ms
      integer, intent(in) :: ratio
         !! `dt_tracer_advect_ratio` for this window; sets the
         !! fixed-budget pass count `max_iter = 2·ratio+1` (ratio is an
         !! integer ⇒ ceil(ratio) = ratio).
      type(ocean_bc_state_t), intent(in), optional :: bc

      integer :: nx, ny, nz, nx_phys, ny_phys, nghost
      integer :: max_iter, ipass, it
      integer :: i, j, k
      real(wp) :: cfl_max, denom, fmax
      logical :: per_x, per_y, fold_n

      if (.not. allocated(ms%tracers)) return

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      nx_phys = grid%nx_phys
      ny_phys = grid%ny_phys
      nghost = grid%nghost

      per_x = .false.
      per_y = .false.
      fold_n = .false.
      if (present(bc)) then
         per_x = bc%periodic_x .and. .not. ocean_halo_is_decomposed_x()
         per_y = bc%periodic_y .and. .not. ocean_halo_is_decomposed_y()
         fold_n = bc%north_fold
      end if

      ! ---- Halo of the accumulators (single-rank seams) ----
      ! hprev reads neighbour uhtr/vhtr across the ghost, so wrap them
      ! first.  TODO(E1): the multi-rank C-grid face halo of uhtr/vhtr
      ! rides the MPI halo work; single-rank periodic/fold is in scope.
      call drain_wrap_face_x(this%uhtr, nx, ny, nz, nx_phys, ny_phys, nghost, &
                             per_x, per_y, fold_n)
      call drain_wrap_face_y(this%vhtr, nx, ny, nz, nx_phys, ny_phys, nghost, &
                             per_x, per_y, fold_n)

      ! ---- Conservative availability limiter (Fox-Kemper × windowed-advect) ----
      ! Tighten the accumulated window transports so every cell's
      ! reconstructed window-start volume stays ≥ areaT·h_min — the
      ! combined FK + resolved transport can otherwise overdraw a thin
      ! z* surface layer (vol < 0), tripping the non-conservative
      ! clamp/hatch below and leaking ~5%/day of tracer.  Limited uhtr/vhtr
      ! feed BOTH the reconstruction AND the sub-cycle seed (drain_copy_3d
      ! below), so the drain still telescopes exactly onto h_end and
      ! Σ(areaT·hTr) is conserved.  Inert (bit-identical) when no cell
      ! would otherwise go vol ≤ areaT·h_min (ratio=1 / FK off / benign
      ! windows).  tr_work is reused as the per-cell scale scratch — it is
      ! overwritten by drain_parabola_* in the sub-cycle before any read.
      ! Fixed-budget FCT tighten (GPU-uniform, no data-dependent exit);
      ! AVAIL_LIMIT_PASS bounds the constraint diffusion across the stencil.
      do ipass = 1, AVAIL_LIMIT_PASS
         call drain_avail_limit(nx, ny, nz, metrics%areaT, this%h_min, &
                                ms%h_layer, this%uhtr, this%vhtr, this%tr_work)
         call drain_wrap_centre(this%tr_work, nx, ny, nz, nx_phys, ny_phys, &
                                nghost, per_x, per_y, fold_n)
         call drain_avail_scale_x(nx, ny, nz, this%tr_work, this%uhtr)
         call drain_avail_scale_y(nx, ny, nz, this%tr_work, this%vhtr)
         call drain_wrap_face_x(this%uhtr, nx, ny, nz, nx_phys, ny_phys, nghost, &
                                per_x, per_y, fold_n)
         call drain_wrap_face_y(this%vhtr, nx, ny, nz, nx_phys, ny_phys, nghost, &
                                per_x, per_y, fold_n)
      end do

      ! ---- Reconstruct hprev = areaT·h_end + div(uhtr) (volume → thickness) ----
      ! + vanishing-layer hatch.  Done on the FROZEN hTr grid: the cell's
      ! window-start thickness, against which the frozen tracer mass is a
      ! consistent concentration.
      call drain_reconstruct_hprev(nx, ny, nz, metrics%areaT, metrics%iareaT, &
                                   ms%h_layer, this%uhtr, this%vhtr, this%hprev_work)
      call drain_wrap_centre(this%hprev_work, nx, ny, nz, nx_phys, ny_phys, nghost, &
                             per_x, per_y, fold_n)

      ! Undo the per-stage concentration hold applied by
      ! `continuity_tracer_step_split`: `hTr` currently carries the frozen
      ! window-start concentration weighted by the CURRENT (window-end)
      ! thickness, and the sub-cycle below needs it weighted by the
      ! window-START thickness `hprev` it just reconstructed.  This restores
      ! exactly the content the drain consumed before the concentration hold
      ! existed, so the drain — and everything it produces — is unchanged.
      if (this%hTr_holds_conc .and. allocated(ms%tracers)) then
         do it = 1, size(ms%tracers)
            if (.not. allocated(ms%tracers(it)%hTr)) cycle
            if (.not. ms%tracers(it)%do_horizontal_advection) cycle
            ! Post-`rk2_average` ⇒ POST_AVERAGE weight (see the parameter's
            ! comment).  Together with the in-stage holds recorded above, the
            ! hold/un-hold pair cancels exactly in the accumulator.
            select case (ms%tracers(it)%budget_id)
            case (TRACER_BUDGET_HEAT)
               call drain_rescale_hTr_budget(nx, ny, nz, this%h_win_start, ms%h_layer, &
                                             DRAIN_BUDGET_POST_AVERAGE_WEIGHT, &
                                             ms%tracers(it)%hTr, ms%heat_budget_horiz_adv)
            case (TRACER_BUDGET_SALT)
               call drain_rescale_hTr_budget(nx, ny, nz, this%h_win_start, ms%h_layer, &
                                             DRAIN_BUDGET_POST_AVERAGE_WEIGHT, &
                                             ms%tracers(it)%hTr, ms%salt_budget_horiz_adv)
            case default
               call drain_rescale_hTr(nx, ny, nz, this%h_win_start, ms%h_layer, &
                                      ms%tracers(it)%hTr)
            end select
         end do
         this%hTr_holds_conc = .false.
      end if

      ! ---- Seam-wrap the frozen tracer mass before the first reconstruction ----
      ! The PPM parabola (drain_parabola_*) reads the seam-adjacent cells'
      ! tracer concentration Tr = hTr/hprev at the ghost band (a ±2 stencil).
      ! hprev is already periodic/fold-consistent (wrapped above), but the
      ! frozen hTr the drain inherits is NOT guaranteed periodic in its ghost
      ! halo (the dynamics wrap h_layer, not necessarily the windowed hTr), so
      ! a seam-physical cell's reconstruction read a different Tr than its
      ! translated interior image — breaking bit-exact translation invariance
      ! (the per-pass updates re-wrap hTr at the end, but the FIRST pass's
      ! parabola already consumed the stale ghosts).  Wrap every advected
      ! tracer's hTr once here so pass 1 sees the same periodic/fold ghosts
      ! the later passes do.  No-op on non-periodic walls (bit-identical).
      do it = 1, size(ms%tracers)
         if (.not. allocated(ms%tracers(it)%hTr)) cycle
         if (.not. ms%tracers(it)%do_horizontal_advection) cycle
         call drain_wrap_centre(ms%tracers(it)%hTr, nx, ny, nz, nx_phys, ny_phys, &
                                nghost, per_x, per_y, fold_n)
      end do

      ! ---- Seed remaining transport from the accumulators ----
      call drain_copy_3d(nx + 1, ny, nz, this%uhtr, this%uhr_x)
      call drain_copy_3d(nx, ny + 1, nz, this%vhtr, this%uhr_y)

      ! ---- Fixed-budget sub-cycle, budget sized to the ACTUAL courant ----
      ! The worst case is `2*ratio+1` (per-step CFL ~ 1), but a real run's
      ! accumulated tracer courant is much smaller, so most of those passes are
      ! exact no-ops (once a face's transport is drained the hup/hlos limiter
      ! yields uhh=0 and the pass changes nothing).  Reduce the max per-cell
      ! accumulated courant ONCE (grid-uniform — every column runs the same
      ! `max_iter`, so no warp divergence; no per-pass early-exit) and size the
      ! fixed loop with the SAME `2*ceil(C)+1` form, capped at the worst case.
      ! Dropping the no-op passes is BIT-IDENTICAL (verified vs `2*ratio+1`).
      cfl_max = 0.0_wp
      do concurrent(k=1:nz, j=1:ny, i=1:nx) reduce(max:cfl_max)
         denom = metrics%areaT(i, j)*max(this%hprev_work(i, j, k), this%h_min)
         fmax = max(abs(this%uhr_x(i, j, k)), abs(this%uhr_x(i + 1, j, k)), &
                    abs(this%uhr_y(i, j, k)), abs(this%uhr_y(i, j + 1, k)))
         cfl_max = max(cfl_max, fmax/denom)
      end do
      ! Passes needed = ceiling(cfl_max / per-pass capacity).  The drain_limit_x
      ! two-test limiter guarantees each pass moves at least 0.5·hup (the
      ! `max(0.5*hup, hup-hlos)` floor), i.e. worst-case capacity 0.5 CFL/pass,
      ! so ceiling(2·cfl_max) passes are provably sufficient (and exact in the
      ! divergent worst case).  This is much tighter than the old
      ! `2*ceiling(cfl_max)+1`: at the default ratio = 1 a stable step has
      ! cfl < 0.5 ⇒ ONE pass (was 3 — two full-grid no-op passes).  Still
      ! grid-uniform (one reduction, no per-pass divergence) and capped at the
      ! worst case `2*ratio+1`.
      max_iter = ceiling(2.0_wp*cfl_max)
      if (max_iter > 2*ratio + 1) max_iter = 2*ratio + 1
      if (max_iter < 1) max_iter = 1
      do ipass = 1, max_iter
         ! Zonal sub-pass --------------------------------------------------
         call drain_limit_x(nx, ny, nz, metrics%areaT, this%h_min, &
                            this%uhr_x, this%hprev_work, this%uhh_x)
         ! Wrap the limited per-pass transport so the swept flux at the
         ! periodic seam / north fold reads the SAME uhh from both sides of
         ! the wrapped face (drain_limit_x writes only interior faces 2..nx
         ! + zeros the array edges; without the wrap the two images of a
         ! seam face disagree ⇒ the flux divergence does not telescope and
         ! Σ(areaT·hTr) leaks at the seam).  No-op on non-periodic walls.
         call drain_wrap_face_x(this%uhh_x, nx, ny, nz, nx_phys, ny_phys, nghost, &
                                per_x, per_y, fold_n)
         do it = 1, size(ms%tracers)
            if (.not. allocated(ms%tracers(it)%hTr)) cycle
            if (.not. ms%tracers(it)%do_horizontal_advection) cycle
            if (this%tracer_recon == TRACER_RECON_PPM) then
               call drain_parabola_x(nx, ny, nz, metrics%wet_T, ms%tracers(it)%hTr, this%hprev_work, &
                                     this%tr_work, this%pal, this%par, this%pa6)
               ! Wrap the parabola coefficients so a seam-face donor that lands
               ! in a ghost column carries the SAME (full-PPM) reconstruction as
               ! its physical image — the ghost band is otherwise PCM, which
               ! makes the two seam-face flux evaluations disagree (seam leak).
               ! Batch the 3 independent parabola-coeff wraps async on queue 1,
               ! sync once before the swept flux reads them.
               call drain_wrap_centre(this%pal, nx, ny, nz, nx_phys, ny_phys, nghost, &
                                      per_x, per_y, fold_n, no_wait=.true.)
               call drain_wrap_centre(this%par, nx, ny, nz, nx_phys, ny_phys, nghost, &
                                      per_x, per_y, fold_n, no_wait=.true.)
               call drain_wrap_centre(this%pa6, nx, ny, nz, nx_phys, ny_phys, nghost, &
                                      per_x, per_y, fold_n, no_wait=.true.)
               !$acc wait(1)
               call drain_swept_flux_x(nx, ny, nz, metrics%areaT, this%uhh_x, &
                                       this%hprev_work, this%pal, this%par, this%pa6, &
                                       this%tr_flux_x)
            else
               ! Q6 WENO drain: reconstruct the swept-average donor-edge
               ! concentration with the WENO ladder instead of the CW
               ! parabola.  Build Tr = hTr/hprev, wrap it so seam-adjacent
               ! stencils read consistent ghosts (the WENO reach is wider
               ! than PPM's — the wrap must cover the widest rung), then
               ! evaluate the swept face flux directly from Tr.
               call drain_fill_conc(nx, ny, nz, ms%tracers(it)%hTr, this%hprev_work, this%tr_work)
               call drain_wrap_centre(this%tr_work, nx, ny, nz, nx_phys, ny_phys, nghost, &
                                      per_x, per_y, fold_n)
               call drain_swept_flux_x_weno(nx, ny, nz, nghost, per_x, metrics%areaT, &
                                            this%uhh_x, this%hprev_work, this%tr_work, &
                                            metrics%wet_T, this%tracer_recon, this%tr_flux_x)
            end if
            select case (ms%tracers(it)%budget_id)
            case (TRACER_BUDGET_HEAT)
               call drain_update_tracer_x_budget(nx, ny, nz, metrics%iareaT, this%tr_flux_x, &
                                                 DRAIN_BUDGET_POST_AVERAGE_WEIGHT, &
                                                 ms%tracers(it)%hTr, ms%heat_budget_horiz_adv)
            case (TRACER_BUDGET_SALT)
               call drain_update_tracer_x_budget(nx, ny, nz, metrics%iareaT, this%tr_flux_x, &
                                                 DRAIN_BUDGET_POST_AVERAGE_WEIGHT, &
                                                 ms%tracers(it)%hTr, ms%salt_budget_horiz_adv)
            case default
               call drain_update_tracer_x(nx, ny, nz, metrics%iareaT, this%tr_flux_x, &
                                          ms%tracers(it)%hTr)
            end select
            call drain_wrap_centre(ms%tracers(it)%hTr, nx, ny, nz, nx_phys, ny_phys, &
                                   nghost, per_x, per_y, fold_n)
         end do
         ! Advance the (shared) thickness + remaining transport once.
         call drain_update_h_x(nx, ny, nz, metrics%iareaT, this%uhh_x, this%hprev_work)
         call drain_subtract_3d(nx + 1, ny, nz, this%uhh_x, this%uhr_x)
         call drain_wrap_centre(this%hprev_work, nx, ny, nz, nx_phys, ny_phys, nghost, &
                                per_x, per_y, fold_n)
         call drain_wrap_face_x(this%uhr_x, nx, ny, nz, nx_phys, ny_phys, nghost, &
                                per_x, per_y, fold_n)

         ! Meridional sub-pass ---------------------------------------------
         call drain_limit_y(nx, ny, nz, metrics%areaT, this%h_min, &
                            this%uhr_y, this%hprev_work, this%uhh_y)
         ! Wrap the limited meridional transport (see the zonal note above).
         call drain_wrap_face_y(this%uhh_y, nx, ny, nz, nx_phys, ny_phys, nghost, &
                                per_x, per_y, fold_n)
         do it = 1, size(ms%tracers)
            if (.not. allocated(ms%tracers(it)%hTr)) cycle
            if (.not. ms%tracers(it)%do_horizontal_advection) cycle
            if (this%tracer_recon == TRACER_RECON_PPM) then
               call drain_parabola_y(nx, ny, nz, metrics%wet_T, ms%tracers(it)%hTr, this%hprev_work, &
                                     this%tr_work, this%pal, this%par, this%pa6)
               ! Wrap the meridional parabola coefficients (see the zonal note).
               ! Batch the 3 independent parabola-coeff wraps async on queue 1,
               ! sync once before the swept flux reads them.
               call drain_wrap_centre(this%pal, nx, ny, nz, nx_phys, ny_phys, nghost, &
                                      per_x, per_y, fold_n, no_wait=.true.)
               call drain_wrap_centre(this%par, nx, ny, nz, nx_phys, ny_phys, nghost, &
                                      per_x, per_y, fold_n, no_wait=.true.)
               call drain_wrap_centre(this%pa6, nx, ny, nz, nx_phys, ny_phys, nghost, &
                                      per_x, per_y, fold_n, no_wait=.true.)
               !$acc wait(1)
               call drain_swept_flux_y(nx, ny, nz, metrics%areaT, this%uhh_y, &
                                       this%hprev_work, this%pal, this%par, this%pa6, &
                                       this%tr_flux_y)
            else
               ! Q6 WENO drain (meridional; see the zonal branch note).
               call drain_fill_conc(nx, ny, nz, ms%tracers(it)%hTr, this%hprev_work, this%tr_work)
               call drain_wrap_centre(this%tr_work, nx, ny, nz, nx_phys, ny_phys, nghost, &
                                      per_x, per_y, fold_n)
               call drain_swept_flux_y_weno(nx, ny, nz, nghost, per_y, metrics%areaT, &
                                            this%uhh_y, this%hprev_work, this%tr_work, &
                                            metrics%wet_T, this%tracer_recon, this%tr_flux_y)
            end if
            select case (ms%tracers(it)%budget_id)
            case (TRACER_BUDGET_HEAT)
               call drain_update_tracer_y_budget(nx, ny, nz, metrics%iareaT, this%tr_flux_y, &
                                                 DRAIN_BUDGET_POST_AVERAGE_WEIGHT, &
                                                 ms%tracers(it)%hTr, ms%heat_budget_horiz_adv)
            case (TRACER_BUDGET_SALT)
               call drain_update_tracer_y_budget(nx, ny, nz, metrics%iareaT, this%tr_flux_y, &
                                                 DRAIN_BUDGET_POST_AVERAGE_WEIGHT, &
                                                 ms%tracers(it)%hTr, ms%salt_budget_horiz_adv)
            case default
               call drain_update_tracer_y(nx, ny, nz, metrics%iareaT, this%tr_flux_y, &
                                          ms%tracers(it)%hTr)
            end select
            call drain_wrap_centre(ms%tracers(it)%hTr, nx, ny, nz, nx_phys, ny_phys, &
                                   nghost, per_x, per_y, fold_n)
         end do
         call drain_update_h_y(nx, ny, nz, metrics%iareaT, this%uhh_y, this%hprev_work)
         call drain_subtract_3d(nx, ny + 1, nz, this%uhh_y, this%uhr_y)
         call drain_wrap_centre(this%hprev_work, nx, ny, nz, nx_phys, ny_phys, nghost, &
                                per_x, per_y, fold_n)
         call drain_wrap_face_y(this%uhr_y, nx, ny, nz, nx_phys, ny_phys, nghost, &
                                per_x, per_y, fold_n)
      end do

      ! ---- Reset accumulators + clock for the next window ----
      call drain_zero_3d(nx + 1, ny, nz, this%uhtr)
      call drain_zero_3d(nx, ny + 1, nz, this%vhtr)
      this%t_dyn_rel_adv = 0.0_wp
   end subroutine continuity_tracer_drain