drain_avail_limit Subroutine

   pure subroutine drain_avail_limit(nx, ny, nz, areaT, h_min, h_end, &
                                     uhtr, vhtr, scratch)
      !! Conservative upfront availability limiter on the accumulated window
      !! transports `uhtr/vhtr`, applied BEFORE drain_reconstruct_hprev.
      !!
      !! Guarantees, for every cell, that the reconstructed window-start
      !! volume stays positive (≥ areaT·h_min):
      !!
      !!   vol(i,j,k) = areaT·h_end + (uhtr(i+1)-uhtr(i)) + (vhtr(j+1)-vhtr(j))
      !!              ≥ areaT·h_min
      !!
      !! so the non-conservative max(0,·) clamp + vanishing-layer hatch in
      !! drain_reconstruct_hprev never fire.  Without this, the combined
      !! Fox-Kemper (2008,2011) overturning + resolved/barotropic transport
      !! can overdraw a thin z* surface layer (vol < 0), the clamp destroys
      !! volume and the hatch re-inflates hprev WITHOUT matching tracer mass
      !! → ~5%/day tracer leak.  MOM6 avoids this by sizing its availability
      !! cap against the combined transport; this is the windowed-drain
      !! analogue (general, protects against any overdraw source).
      !!
      !! Mechanism (Zalesak/FCT-style inflow scaling, GPU-uniform fixed
      !! budget): the cell deficit is covered by shrinking the cell's
      !! INFLOW-face transports.  Reducing an inflow magnitude raises the
      !! receiving cell's vol AND the donor neighbour's vol (it is that
      !! neighbour's outflow), so the iteration is monotone in every cell's
      !! vol and converges.  Scaling is applied multiplicatively to the
      !! transport, so the limited uhtr/vhtr are used CONSISTENTLY for both
      !! the reconstruction and the sub-cycle ⇒ the drain still telescopes
      !! exactly onto h_end and Σ(areaT·hTr) is conserved to round-off.
      !!
      !! Inert when every vol > areaT·h_min already (scale = 1 everywhere),
      !! so ratio=1 / FK-off / non-overdrawing windows are bit-identical.
      !!
      !! ONE FCT pass — the caller loops it `AVAIL_MAX_PASS` times with a
      !! seam re-wrap of `scratch` and the faces between passes so the
      !! constraint diffuses consistently across periodic/fold ghosts.
      integer, intent(in) :: nx, ny, nz
      real(wp), intent(in) :: areaT(nx, ny)
      real(wp), intent(in) :: h_min
      real(wp), intent(in) :: h_end(nx, ny, nz)
      real(wp), intent(inout) :: uhtr(nx + 1, ny, nz)
      real(wp), intent(inout) :: vhtr(nx, ny + 1, nz)
      real(wp), intent(inout) :: scratch(nx, ny, nz)
         !! Per-cell inflow scale factor in [0,1] (centre-shaped slot).
      integer :: i, j, k
      real(wp) :: vol, vmin, inflow, deficit, sfac
      ! ---- Per-cell inflow scale factor ----
      do concurrent(k=1:nz, j=1:ny, i=1:nx) &
         local(vol, vmin, inflow, deficit)
         vmin = areaT(i, j)*h_min
         vol = areaT(i, j)*h_end(i, j, k) &
               + (uhtr(i + 1, j, k) - uhtr(i, j, k)) &
               + (vhtr(i, j + 1, k) - vhtr(i, j, k))
         ! Inflow into this cell across its four faces (volume, ≥ 0):
         !   west  face uhtr(i)   inflow if > 0
         !   east  face uhtr(i+1) inflow if < 0
         !   south face vhtr(j)   inflow if > 0
         !   north face vhtr(j+1) inflow if < 0
         inflow = max(0.0_wp, uhtr(i, j, k)) &
                  + max(0.0_wp, -uhtr(i + 1, j, k)) &
                  + max(0.0_wp, vhtr(i, j, k)) &
                  + max(0.0_wp, -vhtr(i, j + 1, k))
         if (vol < vmin .and. inflow > 0.0_wp) then
            deficit = vmin - vol
            scratch(i, j, k) = max(0.0_wp, (inflow - deficit)/inflow)
         else
            scratch(i, j, k) = 1.0_wp
         end if
      end do
   end subroutine drain_avail_limit