bt_cfl_dt_wet Subroutine

   pure subroutine bt_cfl_dt_wet(nx, ny, i0, i1, j0, j1, b, wet, dxT, dyT, &
                                 cfl_safety, dt_bt, h_at, l_at, n_wet)
      !! Per-WET-CELL external-gravity-wave CFL limit (MOM6 `set_dtbt`):
      !!
      !!     dt_bt = min over wet (i,j) of  cfl_safety * l(i,j) / c(i,j),
      !!     l(i,j) = 1/sqrt(1/dxT(i,j)^2 + 1/dyT(i,j)^2),
      !!     c(i,j) = sqrt(g * max(b(i,j), 1 m)),
      !!
      !! i.e. the LOCAL depth with the LOCAL cell size, over OCEAN only.
      !! MOM6 evaluates the same quantity as `gtot*dt^2*(1/dx^2+1/dy^2)`
      !! per wet point.
      !!
      !! **Why per point.** The former estimate combined the deepest depth
      !! ANYWHERE with the smallest cell ANYWHERE — on the 1° tripolar grid
      !! that was 6000 m of ocean against a 362 m LAND cell at a
      !! land-locked bipole, and bought 1930 substeps where ~32 suffice.
      !! A land cell carries no gravity wave, and a small cell over a
      !! shallow shelf does not see the abyssal wave speed.
      !!
      !! **Bit-identity where the two agree.** Each point's value is
      !! evaluated with EXACTLY the arithmetic the global-extremes estimate
      !! used (`1/sqrt(inv_l2)`, `sqrt(g*max(b,1))`, `safety*l/c`, in that
      !! order).  Every step is monotone under round-to-nearest, so where
      !! the deepest wet column and the smallest wet cell COINCIDE (a
      !! flat-bottomed uniform grid, a wall basin with no land) the minimum
      !! is attained at that point and equals the old number bit-for-bit.
      !!
      !! No wet cell in the window ⇒ `dt_bt = huge`, `n_wet = 0` — the
      !! identity of the cross-rank `min` reduction (a rank that is all
      !! land does not constrain the others).
      !!
      !! Host-side, configure time: plain loops, no `do concurrent`.
      integer, intent(in) :: nx
         !! First extent of the centre arrays (ghosts included).
      integer, intent(in) :: ny
         !! Second extent of the centre arrays (ghosts included).
      integer, intent(in) :: i0
         !! First physical i.
      integer, intent(in) :: i1
         !! Last physical i.
      integer, intent(in) :: j0
         !! First physical j.
      integer, intent(in) :: j1
         !! Last physical j.
      real(wp), intent(in) :: b(nx, ny)
         !! Bed depth (m, positive down).
      real(wp), intent(in) :: wet(nx, ny)
         !! Static wet (1) / land (0) T-cell mask.
      real(wp), intent(in) :: dxT(nx, ny)
         !! T-cell x length (m).
      real(wp), intent(in) :: dyT(nx, ny)
         !! T-cell y length (m).
      real(wp), intent(in) :: cfl_safety
         !! `&ocean_bt_nml cfl_bt_safety`.
      real(wp), intent(out) :: dt_bt
         !! Smallest per-wet-cell safe substep (s); `huge` if no wet cell.
      real(wp), intent(out) :: h_at
         !! Bed depth at the limiting cell (m); 0 if no wet cell.
      real(wp), intent(out) :: l_at
         !! 2-D CFL length at the limiting cell (m); 0 if no wet cell.
      integer, intent(out) :: n_wet
         !! Wet cells scanned.
      integer :: i, j
      real(wp) :: inv_l2, l_ij, c_ij, dt_ij

      dt_bt = huge(1.0_wp)
      h_at = 0.0_wp
      l_at = 0.0_wp
      n_wet = 0
      do j = j0, j1
         do i = i0, i1
            if (wet(i, j) <= 0.5_wp) cycle
            n_wet = n_wet + 1
            inv_l2 = 1.0_wp/dxT(i, j)**2 + 1.0_wp/dyT(i, j)**2
            l_ij = 1.0_wp/sqrt(inv_l2)
            c_ij = sqrt(GRAVITY*max(b(i, j), 1.0_wp))
            dt_ij = cfl_safety*l_ij/c_ij
            if (dt_ij < dt_bt) then
               dt_bt = dt_ij
               h_at = b(i, j)
               l_at = l_ij
            end if
         end do
      end do
   end subroutine bt_cfl_dt_wet