pd_limit_zonal_impl Subroutine

   subroutine pd_limit_zonal_impl(nx, ny, nz, dt, h_lim, iareaT, h_layer, &
                                  mass_flux_x, theta, n_limited, u_cor)
      !! Positive-definite per-donor outflux limiter — zonal (x) pass (P2,
      !! plan §P2 / decision D2).  Scales the per-layer east-face mass fluxes
      !! DOWN so no donor cell loses more thickness than it holds above the
      !! floor `h_lim`: guarantees `h_layer >= h_lim` after
      !! `continuity_apply_zonal`, with ZERO mass created — outfluxes shrink,
      !! thickness is never inflated (the deliberate contrast with MOM6's
      !! `max(h, Angstrom)` injection; the conservative borrow stays the
      !! backstop).  Two device passes:
      !!
      !!   1. per-cell `θ(i,j,k) = min(1, avail/demand)` over the FULL index
      !!      range (incl. ghosts), with `demand = dt·outflow·iareaT` the
      !!      thickness (m) the two x-faces would drain and
      !!      `avail = max(h − h_lim, 0)`.  `outflow` counts only the OUTGOING
      !!      part of each face — east face `i+1` when its flux is positive,
      !!      west face `i` when its flux is negative — because in a Lie pass
      !!      only outgoing faces drain the cell.  θ construction (the
      !!      `merge(avail/max(demand,H_DIV_EPS), 1, demand>avail)` form with
      !!      pure 1/0 armour) mirrors the barotropic wet/dry sweep.
      !!   2. each interior face `i ∈ 2..nx` is multiplied by its UPWIND
      !!      donor's θ (west cell `i−1` when the face flux ≥ 0, else east
      !!      cell `i`) — donor-side only, NO two-cell min: a face's mass
      !!      leaves exactly one donor per direction pass, and that donor's
      !!      OTHER outgoing face is scaled by the SAME θ, so realised total
      !!      outflow ≤ `avail` is already guaranteed by θ_donor alone.
      !!
      !! The θ range spans ghosts so a seam/periodic donor in a ghost column
      !! carries a current factor.  MPI-seam determinism (deferred multi-rank):
      !! two ranks sharing a seam face compute the same θ from the same
      !! exchanged donor `h_layer` and the same face flux, so they scale the
      !! shared face identically.  `n_limited` accumulates the count of scaled
      !! faces (θ_donor < 1, flux ≠ 0) via an `!$acc parallel loop reduction`
      !! (a `do concurrent` + `sum()` would silently read the stale host shadow
      !! on the mem:separate GPU build).
      !!
      !! v1.1 — u_cor re-matching (optional `u_cor`).  `split_scheme="pred_corr"`
      !! captures the MOM6 `u_cor` (the transport-matched velocity the fluxes
      !! correspond to) inside `continuity_zonal_flux` BEFORE this limiter runs;
      !! the `use_state_fluxes` corrector then evaluates at `u_cor` but transports
      !! with the LIMITED `mass_flux_x`.  After a face is scaled by θ its captured
      !! `u_cor` no longer corresponds — `flux = h_face·u·metric` with `h_face`
      !! fixed ⇒ `u` scales linearly with θ, so `u_cor *= th_face` restores the
      !! flux↔velocity match.  Left unfixed this is a flux/velocity
      !! inconsistency (anti-damping class) that surfaces as a late-onset
      !! CFL blow-up on the split-scheme path.  The re-scale runs in
      !! a `do concurrent` guarded by a host `present` flag (the renorm's proven
      !! optional-device-array pattern) BEFORE the mass-flux sweep, so both read
      !! the same UNSCALED-flux donor sign.  Absent ⇒ bit-identical to v1.
      integer, intent(in) :: nx, ny, nz
      real(wp), intent(in) :: dt, h_lim
      real(wp), intent(in) :: iareaT(nx, ny)
      real(wp), intent(in) :: h_layer(nx, ny, nz)
      real(wp), intent(inout) :: mass_flux_x(nx + 1, ny, nz)
      real(wp), intent(inout) :: theta(nx, ny, nz)
      integer, intent(inout) :: n_limited
      real(wp), intent(inout), optional :: u_cor(nx + 1, ny, nz)
         !! MOM6 `u_cor` capture (transport-matched velocity); when present,
         !! re-scaled by the SAME per-face θ as `mass_flux_x` so it stays
         !! consistent with the limited flux the mom6-scheme corrector reads.

      integer :: i, j, k, n_lim
      real(wp) :: out_e, out_w, outflow, demand, avail, th_face
      logical :: do_ucor

      do_ucor = present(u_cor)

      do concurrent(k=1:nz, j=1:ny, i=1:nx) &
         local(out_e, out_w, outflow, demand, avail)
         out_e = max(mass_flux_x(i + 1, j, k), 0.0_wp)
         out_w = max(-mass_flux_x(i, j, k), 0.0_wp)
         outflow = out_e + out_w
         demand = dt*outflow*iareaT(i, j)
         avail = max(h_layer(i, j, k) - h_lim, 0.0_wp)
         theta(i, j, k) = merge(avail/max(demand, H_DIV_EPS), 1.0_wp, demand > avail)
      end do

      ! v1.1: re-match u_cor to the limited flux, reading the UNSCALED-flux
      ! donor sign (runs before the mass-flux sweep below).  Host-flag-guarded
      ! do concurrent = the renorm's proven optional-device-array pattern.
      if (do_ucor) then
         do concurrent(k=1:nz, j=1:ny, i=2:nx) local(th_face)
            if (mass_flux_x(i, j, k) >= 0.0_wp) then
               th_face = theta(i - 1, j, k)
            else
               th_face = theta(i, j, k)
            end if
            u_cor(i, j, k) = u_cor(i, j, k)*th_face
         end do
      end if

      n_lim = 0
      do concurrent(k=1:nz, j=1:ny, i=2:nx) local(th_face) reduce(+:n_lim)
         if (mass_flux_x(i, j, k) >= 0.0_wp) then
            th_face = theta(i - 1, j, k)
         else
            th_face = theta(i, j, k)
         end if
         if (th_face < 1.0_wp .and. mass_flux_x(i, j, k) /= 0.0_wp) n_lim = n_lim + 1
         mass_flux_x(i, j, k) = mass_flux_x(i, j, k)*th_face
      end do
      n_limited = n_limited + n_lim
   end subroutine pd_limit_zonal_impl