apply_velocity_truncation Subroutine

   subroutine apply_velocity_truncation(ms, metrics, dt, cfl_trunc, maxvel, ntrunc_step, &
                                        vanish_tol, n_nanzero, clip_cell_metric, &
                                        nan_i, nan_j, nan_k, nan_is_u, nan_dx, nan_visc_cfl, nu_h)
      !! Post-RK2 velocity housekeeping: the advective-CFL truncation
      !! (E7) followed by the absolute `maxvel` cap.  Called once at the
      !! end of each outer step (after the RK2 average, before the ALE
      !! remap) so the carried-forward / remapped velocity field is
      !! bounded.  Replaces the bare `apply_maxvel_clamp` call at both
      !! `ocean_dyn_step` / `ocean_dyn_step_split` sites.
      !!
      !! Order (well-defined, idempotent):
      !!   0. **Vanished-layer zero** (Phase 3, gated `vanish_tol > 0`):
      !!      a u-face where BOTH adjacent centre thicknesses are at or
      !!      below `vanish_tol` is set to 0 and NOT counted as a CFL
      !!      truncation (it is a diagnostic-clean zero, not a real clip).
      !!      A face with at least one massive side is untouched (R1).
      !!      Runs BEFORE the CFL clip so a zeroed face cannot re-trigger
      !!      the clip.  Idempotent with Phase-2 reset (R6).
      !!   1. **CFL clip** (gated `cfl_trunc > 0`): any face whose local
      !!      advective CFL `|u|·dt·idx` exceeds `cfl_trunc` is reset to
      !!      `sign(0.9·cfl_trunc/(dt·idx), u_old)` — i.e. magnitude
      !!      `0.9·cfl_trunc·dx/dt`, sign preserved.  The `0.9`
      !!      relaxation (`CFL_TRUNC_RELAX`) keeps the clipped face below
      !!      threshold so it does not re-trip on round-off next step.
      !!      idx is the stored reciprocal metric (`idxCu`/`idyCv`).
      !!   2. **maxvel cap** (gated `maxvel > 0`, via `apply_maxvel_clamp`):
      !!      the coarser absolute physical backstop.  Runs second, so a
      !!      face clipped to `0.9·cfl_trunc·dx/dt > maxvel` is further
      !!      bounded to `±maxvel`.
      !!
      !! The CFL clip + count run in a single `!$acc parallel loop
      !! reduction(+:n)` per face component — *not* a bare `do concurrent`
      !! + `sum()` over scratch (which silently returns 0 on GPU under
      !! stdpar; see the no-managed-memory gotcha).  `ntrunc_step` is the
      !! host-scalar reduction result (reset to 0 here each call).
      !!
      !! Not conservative — clipping a runaway face loses momentum; when
      !! the truncation fires, conservation is already broken by whatever
      !! produced the runaway.  PointAccel-style storage of the worst
      !! offender's (i,j,k) + acceleration breakdown is a deferred future
      !! extension; v1 ships the count only.
      !!
      !! `vanish_tol` absent / <= 0 ⇒ step 0 skipped ⇒ bit-identical.
      type(multilayer_state_t), intent(inout) :: ms
      type(ocean_metrics_t), intent(in) :: metrics
      real(wp), intent(in) :: dt
      real(wp), intent(in) :: cfl_trunc
      real(wp), intent(in) :: maxvel
      integer, intent(out) :: ntrunc_step
      real(wp), intent(in), optional :: vanish_tol
         !! When present (> 0): zero both-sided-vanished faces before the
         !! CFL clip.  Absent / 0 ⇒ skip ⇒ bit-identical.

      real(wp), parameter :: CFL_TRUNC_RELAX = 0.9_wp
         !! MOM6 relaxation factor: clip to 90% of the threshold so the
         !! clipped face stays sub-threshold and does not re-trip.
      integer :: nx_uface, ny_uface, nx_vface, ny_vface, nz
      integer :: nx_centre, ny_centre
      integer, intent(out), optional :: n_nanzero
         !! Count of non-finite (NaN/Inf) face velocities zeroed by the
         !! step -1 NaN-catch this call.  A non-finite velocity is a
         !! producer bug (a 0/0 upstream); it MUST be caught here so it
         !! can never launder to ±maxvel (nvfortran lowers the
         !! `if(u>hi)…else if(u<lo)` pair to a NaN-skipping min/max clamp
         !! on GPU ⇒ NaN → -maxvel).  Loud counter — 0 on a healthy run.
      logical, intent(in), optional :: clip_cell_metric
         !! When present AND .true.: the CFL clip bounds each face on the
         !! CELL metric `max(idxT(i-1,j), idxT(i,j))` (v: `idyT`) — the SAME
         !! metric the console panic / `compute_max_cfl` uses — instead of
         !! the face metric `idxCu`/`idyCv` (pdc 6991988f).  Guarantees the
         !! panic-visible CFL is bounded even where the face metric is
         !! anomalous at a grounding face (`|u|·dt·idxCu ≤ cfl_trunc` while
         !! `|u|·dt·idxT > cfl_trunc` — the θ-edge escape; measured at
         !! 1024²/dt=800 as MaxCFL 0.689 sailing past a 0.5 ceiling).
         !! Bit-identical on grids where `idxCu == idxT`.  Absent /
         !! .false. ⇒ face metric ⇒ bit-identical.
      integer, intent(out), optional :: nan_i, nan_j, nan_k
         !! Grid location (local, including ghosts) of the FIRST non-finite
         !! face this call caught, in (k,j,i)-ascending scan order —
         !! actionable in place of the old bare count ("producer 0/0
         !! upstream — investigate"; FINDINGS.md's global-tripolar-aquaplanet
         !! debugging session had nothing better to go on).  0 when
         !! `n_nanzero` is absent/0 (nothing to report) or when none of
         !! these outputs were requested (the search is skipped entirely —
         !! see `nan_dx`/`nan_visc_cfl`).
      logical, intent(out), optional :: nan_is_u
         !! `.true.` if the first non-finite face was a u-face
         !! (`u_face_x_layer`), `.false.` if a v-face. Only meaningful
         !! when `nan_i`/`nan_j`/`nan_k` were actually found (`n_nan > 0`
         !! AND the location search ran).
      real(wp), intent(out), optional :: nan_dx
         !! Local cell size (m) at `(nan_i, nan_j)` — `1/max(idxT(i-1,j),
         !! idxT(i,j))` (v-face: idyT analogue), the SAME cell metric
         !! `clip_cell_metric` uses. Cheap (one extra metrics read at the
         !! single located cell, not a scan) — always computed alongside
         !! the location when a location was found.
      real(wp), intent(out), optional :: nan_visc_cfl
         !! Local viscous CFL `nu_h*dt/nan_dx^2` at the located cell —
         !! only computed when `nu_h` is supplied (the caller's
         !! `ocean_horizontal_viscosity_t%nu_h`); 0 otherwise. Mirrors
         !! `rdb_ocean_stability_audit`'s configure-time check, evaluated
         !! HERE at the exact runtime location the NaN was caught, so a
         !! user reading the message can immediately see whether an
         !! under-resolved viscous CFL is the likely producer.
      real(wp), intent(in), optional :: nu_h
         !! Constant horizontal viscosity (m^2/s), for `nan_visc_cfl`
         !! only. Absent ⇒ `nan_visc_cfl` (if requested) is 0.
      integer :: i, j, k, n_u, n_v, n_nan
      real(wp) :: idx_use, idy_use
      logical :: do_clip_cell
      real(wp) :: vt
      logical :: do_vanish_zero
      logical :: want_nan_loc
      integer(int64) :: lin_u_min, lin_v_min, lin
      integer(int64) :: nxu64, nyu64, nxv64, nyv64

      do_clip_cell = .false.
      if (present(clip_cell_metric)) do_clip_cell = clip_cell_metric
      ntrunc_step = 0
      do_vanish_zero = .false.
      vt = 0.0_wp
      if (present(vanish_tol)) then
         if (vanish_tol > 0.0_wp) then
            do_vanish_zero = .true.
            vt = vanish_tol
         end if
      end if

      ! Step -1: NaN/Inf CATCH (loud, unconditional; pdc 8c2fd674).  A
      ! non-finite face velocity — born of a 0/0 in an upstream producer
      ! (implicit vdiff on an all-at-floor column, a per-column BT-fold/PGF
      ! division at a grounding θ-edge) — is NaN-false to the CFL clip's
      ! `abs(u)…>cfl_trunc` (no clip) and LAUNDERS to ±maxvel in
      ! apply_maxvel_clamp (nvfortran lowers `if(u>hi)…else if(u<lo)…` to a
      ! NaN-skipping min/max clamp on GPU).  It then transports mass at
      ! maxvel.  Zero it here — BEFORE clip + maxvel — and count it loudly
      ! (max-reductions are NaN-blind, so the count is explicit).  Split
      ! count/zero (write out of the reduction loop, per the clip fix).
      ! Finite fields ⇒ 0 non-finite ⇒ bit-identical.
      nz = ms%nz_ml
      n_nan = 0
      nx_uface = size(ms%u_face_x_layer, 1)
      ny_uface = size(ms%u_face_x_layer, 2)
      nx_vface = size(ms%v_face_y_layer, 1)
      ny_vface = size(ms%v_face_y_layer, 2)
      !$acc parallel loop collapse(3) reduction(+:n_nan) present(ms%u_face_x_layer)
      do k = 1, nz
         do j = 1, ny_uface
            do i = 1, nx_uface
               if (.not. ieee_is_finite(ms%u_face_x_layer(i, j, k))) n_nan = n_nan + 1
            end do
         end do
      end do
      !$acc parallel loop collapse(3) reduction(+:n_nan) present(ms%v_face_y_layer)
      do k = 1, nz
         do j = 1, ny_vface
            do i = 1, nx_vface
               if (.not. ieee_is_finite(ms%v_face_y_layer(i, j, k))) n_nan = n_nan + 1
            end do
         end do
      end do
      ! Locate the FIRST non-finite face (min-reduction over an encoded
      ! linear index, k-major) — ONLY when n_nan>0 (something is already
      ! broken) AND the caller actually asked for a location (any of
      ! nan_i/nan_j/nan_k/nan_is_u/nan_dx/nan_visc_cfl present). This is
      ! the per-step-path cost gate: a healthy run (n_nan==0, the
      ! overwhelming common case) never runs this block at all, and a
      ! caller not asking for a location (the legacy call sites) pays
      ! nothing either. MUST run BEFORE the zero-out below — the zeroed
      ! array has nothing left to locate.
      want_nan_loc = present(nan_i) .or. present(nan_j) .or. present(nan_k) .or. &
                     present(nan_is_u) .or. present(nan_dx) .or. present(nan_visc_cfl)
      if (present(nan_i)) nan_i = 0
      if (present(nan_j)) nan_j = 0
      if (present(nan_k)) nan_k = 0
      if (present(nan_is_u)) nan_is_u = .true.
      if (present(nan_dx)) nan_dx = 0.0_wp
      if (present(nan_visc_cfl)) nan_visc_cfl = 0.0_wp
      if (n_nan > 0 .and. want_nan_loc) then
         nxu64 = int(nx_uface, int64)
         nyu64 = int(ny_uface, int64)
         nxv64 = int(nx_vface, int64)
         nyv64 = int(ny_vface, int64)
         lin_u_min = huge(1_int64)
         !$acc parallel loop collapse(3) reduction(min:lin_u_min) present(ms%u_face_x_layer)
         do k = 1, nz
            do j = 1, ny_uface
               do i = 1, nx_uface
                  if (.not. ieee_is_finite(ms%u_face_x_layer(i, j, k))) then
                     lin_u_min = min(lin_u_min, &
                                     (int(k - 1, int64)*nyu64 + int(j - 1, int64))*nxu64 + int(i - 1, int64))
                  end if
               end do
            end do
         end do
         lin_v_min = huge(1_int64)
         !$acc parallel loop collapse(3) reduction(min:lin_v_min) present(ms%v_face_y_layer)
         do k = 1, nz
            do j = 1, ny_vface
               do i = 1, nx_vface
                  if (.not. ieee_is_finite(ms%v_face_y_layer(i, j, k))) then
                     lin_v_min = min(lin_v_min, &
                                     (int(k - 1, int64)*nyv64 + int(j - 1, int64))*nxv64 + int(i - 1, int64))
                  end if
               end do
            end do
         end do

         block
            integer :: fi, fj, fk
            logical :: found_u
            real(wp) :: idx_local, dx_local
            integer(int64) :: lin_dec, lin_rem

            ! n_nan>0 guarantees at least one of lin_u_min/lin_v_min is a
            ! real (non-huge) encoded index; a tie prefers u.
            found_u = lin_u_min <= lin_v_min
            if (found_u) then
               lin_dec = lin_u_min
               fi = int(mod(lin_dec, nxu64), kind(fi)) + 1
               lin_rem = lin_dec/nxu64
               fj = int(mod(lin_rem, nyu64), kind(fj)) + 1
               fk = int(lin_rem/nyu64, kind(fk)) + 1
            else
               lin_dec = lin_v_min
               fi = int(mod(lin_dec, nxv64), kind(fi)) + 1
               lin_rem = lin_dec/nxv64
               fj = int(mod(lin_rem, nyv64), kind(fj)) + 1
               fk = int(lin_rem/nyv64, kind(fk)) + 1
            end if

            if (present(nan_i)) nan_i = fi
            if (present(nan_j)) nan_j = fj
            if (present(nan_k)) nan_k = fk
            if (present(nan_is_u)) nan_is_u = found_u

            if (present(nan_dx) .or. present(nan_visc_cfl)) then
               ! Host scalar read of metrics%idxT/idyT: SAFE under
               ! mem:separate despite metrics being device-resident —
               ! these are static, `copyin`-once metric arrays (never
               ! re-written after configure_ocean_metrics), so the host
               ! and device copies stay identical for the life of the
               ! run; no `!$acc update self` needed (contrast with a
               ! per-step state array, which would be stale here).
               nx_centre = size(ms%h_layer, 1)
               ny_centre = size(ms%h_layer, 2)
               if (found_u) then
                  idx_local = max(metrics%idxT(max(fi - 1, 1), fj), &
                                  metrics%idxT(min(fi, nx_centre), fj))
               else
                  idx_local = max(metrics%idyT(fi, max(fj - 1, 1)), &
                                  metrics%idyT(fi, min(fj, ny_centre)))
               end if
               dx_local = 0.0_wp
               if (idx_local > 0.0_wp) dx_local = 1.0_wp/idx_local
               if (present(nan_dx)) nan_dx = dx_local
               if (present(nan_visc_cfl)) then
                  if (present(nu_h) .and. dx_local > 0.0_wp) then
                     nan_visc_cfl = nu_h*dt/dx_local**2
                  end if
               end if
            end if
         end block
      end if

      if (n_nan > 0) then
         do concurrent(k=1:nz, j=1:ny_uface, i=1:nx_uface)
            if (.not. ieee_is_finite(ms%u_face_x_layer(i, j, k))) ms%u_face_x_layer(i, j, k) = 0.0_wp
         end do
         do concurrent(k=1:nz, j=1:ny_vface, i=1:nx_vface)
            if (.not. ieee_is_finite(ms%v_face_y_layer(i, j, k))) ms%v_face_y_layer(i, j, k) = 0.0_wp
         end do
      end if
      if (present(n_nanzero)) n_nanzero = n_nan

      ! Step 0: zero both-sided-vanished faces (Phase 3, BEFORE CFL clip).
      ! A face zeroed here does NOT increment ntrunc_step (it is not a CFL
      ! clip — it is a diagnostic-clean zero of a face with no real mass).
      ! Data-parallel masked write (no reduction) ⇒ do concurrent with the
      ! both-sided-vanished test folded into the loop mask; mirrors
      ! reset_vanished_layer_velocities.  ms arrays are device-resident.
      if (do_vanish_zero) then
         nx_uface = size(ms%u_face_x_layer, 1)
         ny_uface = size(ms%u_face_x_layer, 2)
         nx_vface = size(ms%v_face_y_layer, 1)
         ny_vface = size(ms%v_face_y_layer, 2)
         nx_centre = size(ms%h_layer, 1)
         ny_centre = size(ms%h_layer, 2)
         nz = ms%nz_ml
         ! u-faces: i straddles centres i-1 and i; safe range i=2..nx_uface-1
         ! (both-sided-vanished test as an inner `if`, not a DC mask — masked
         ! headers fail tools/dc_audit.py --strict)
         do concurrent(k=1:nz, j=1:ny_uface, i=2:nx_uface - 1)
            if (max(ms%h_layer(i - 1, j, k), ms%h_layer(i, j, k)) <= vt) then
               ms%u_face_x_layer(i, j, k) = 0.0_wp
            end if
         end do
         ! v-faces: j straddles centres j-1 and j; safe range j=2..ny_vface-1
         do concurrent(k=1:nz, j=2:ny_vface - 1, i=1:nx_vface)
            if (max(ms%h_layer(i, j - 1, k), ms%h_layer(i, j, k)) <= vt) then
               ms%v_face_y_layer(i, j, k) = 0.0_wp
            end if
         end do
      end if

      if (cfl_trunc > 0.0_wp) then
         nx_uface = size(ms%u_face_x_layer, 1)
         ny_uface = size(ms%u_face_x_layer, 2)
         nx_vface = size(ms%v_face_y_layer, 1)
         ny_vface = size(ms%v_face_y_layer, 2)
         nz = ms%nz_ml
         n_u = 0
         n_v = 0

         ! u-faces (Cu): clip in TWO passes — a read-only COUNT reduction,
         ! then a masked-write `do concurrent` CLIP.  A conditional array
         ! WRITE fused into an `!$acc parallel loop reduction` silently
         ! fails to land on the cc70 600² build (the 6 m/s face survived
         ! an unconditionally-true clip — the truncation was partially
         ! inert, "truncations climb but don't save").  Both
         ! passes read the ORIGINAL u ⇒ same faces, same ntrunc_step ⇒
         ! bit-identical to a correct fused loop.
         ! Clip metric: face idxCu by default; the CELL metric
         ! max(idxT(i-1), idxT(i)) — the panic basis — under
         ! clip_cell_metric.  On idxCu==idxT grids the two are identical
         ! ⇒ bit-identical.
         nx_centre = size(ms%h_layer, 1)
         ny_centre = size(ms%h_layer, 2)
         !$acc parallel loop collapse(3) reduction(+:n_u) private(idx_use) &
         !$acc   present(ms%u_face_x_layer, metrics%idxCu, metrics%idxT)
         do k = 1, nz
            do j = 1, ny_uface
               do i = 1, nx_uface
                  if (do_clip_cell) then
                     idx_use = max(metrics%idxT(max(i - 1, 1), j), &
                                   metrics%idxT(min(i, nx_centre), j))
                  else
                     idx_use = metrics%idxCu(i, j)
                  end if
                  if (abs(ms%u_face_x_layer(i, j, k))*dt*idx_use > cfl_trunc) then
                     n_u = n_u + 1
                  end if
               end do
            end do
         end do
         do concurrent(k=1:nz, j=1:ny_uface, i=1:nx_uface) local(idx_use)
            if (do_clip_cell) then
               idx_use = max(metrics%idxT(max(i - 1, 1), j), &
                             metrics%idxT(min(i, nx_centre), j))
            else
               idx_use = metrics%idxCu(i, j)
            end if
            if (abs(ms%u_face_x_layer(i, j, k))*dt*idx_use > cfl_trunc) then
               ms%u_face_x_layer(i, j, k) = sign( &
                                            CFL_TRUNC_RELAX*cfl_trunc/(dt*idx_use), &
                                            ms%u_face_x_layer(i, j, k))
            end if
         end do

         ! v-faces (Cv): same two-pass split as the u-faces above.
         !$acc parallel loop collapse(3) reduction(+:n_v) private(idy_use) &
         !$acc   present(ms%v_face_y_layer, metrics%idyCv, metrics%idyT)
         do k = 1, nz
            do j = 1, ny_vface
               do i = 1, nx_vface
                  if (do_clip_cell) then
                     idy_use = max(metrics%idyT(i, max(j - 1, 1)), &
                                   metrics%idyT(i, min(j, ny_centre)))
                  else
                     idy_use = metrics%idyCv(i, j)
                  end if
                  if (abs(ms%v_face_y_layer(i, j, k))*dt*idy_use > cfl_trunc) then
                     n_v = n_v + 1
                  end if
               end do
            end do
         end do
         do concurrent(k=1:nz, j=1:ny_vface, i=1:nx_vface) local(idy_use)
            if (do_clip_cell) then
               idy_use = max(metrics%idyT(i, max(j - 1, 1)), &
                             metrics%idyT(i, min(j, ny_centre)))
            else
               idy_use = metrics%idyCv(i, j)
            end if
            if (abs(ms%v_face_y_layer(i, j, k))*dt*idy_use > cfl_trunc) then
               ms%v_face_y_layer(i, j, k) = sign( &
                                            CFL_TRUNC_RELAX*cfl_trunc/(dt*idy_use), &
                                            ms%v_face_y_layer(i, j, k))
            end if
         end do

         ntrunc_step = n_u + n_v
      end if

      ! Absolute physical backstop (MOM6 MAXVEL); no-op when maxvel <= 0.
      call apply_maxvel_clamp(ms, maxvel)
   end subroutine apply_velocity_truncation