compute_max_cfl Function

   function compute_max_cfl(u_face, v_face, idxT, idyT, dt, nghost, h_layer, vanish_tol) &
      result(max_cfl)
      !! max (|u_c|·dt·idxT + |v_c|·dt·idyT) over PHYSICAL cells (ghosts
      !! excluded). idxT/idyT are metric inverses (= 1/dx,1/dy on uniform).
      !!
      !! Phase-3 optional gate: when `h_layer` + `vanish_tol` are BOTH
      !! present, cells where `h_layer(i,j,k) <= vanish_tol` are skipped.
      !! A vanished cell carries no real momentum; its face-averaged
      !! velocity spike should not trigger a CFL panic.  When absent
      !! (default) ⇒ un-gated path runs verbatim ⇒ bit-identical.
      !!
      !! Two-loop form (un-gated / gated dispatched externally before
      !! calling) preferred over referencing an absent optional inside the
      !! `!$acc parallel loop` region — keeps the parallel body clean.
      real(wp), intent(in) :: u_face(:, :, :), v_face(:, :, :)
      real(wp), intent(in) :: idxT(:, :), idyT(:, :)
      real(wp), intent(in) :: dt
      integer, intent(in) :: nghost
      real(wp), intent(in), optional :: h_layer(:, :, :)
         !! Centre-cell thickness (m). Required together with `vanish_tol`.
      real(wp), intent(in), optional :: vanish_tol
         !! Cells with h_layer <= vanish_tol are excluded. Required
         !! together with `h_layer`.
      real(wp) :: max_cfl
      real(wp) :: acc
      integer :: i, j, k, nx, ny, nz, i_lo, i_hi, j_lo, j_hi
      logical :: gate

      nx = min(size(u_face, 1) - 1, size(v_face, 1), size(idxT, 1))
      ny = min(size(u_face, 2), size(v_face, 2) - 1, size(idyT, 2))
      nz = min(size(u_face, 3), size(v_face, 3))
      i_lo = nghost + 1
      i_hi = nx - nghost
      j_lo = nghost + 1
      j_hi = ny - nghost
      gate = present(h_layer) .and. present(vanish_tol)

      if (gate) then
         ! Gated path: skip cells where h_layer <= vanish_tol.
         ! Separate loop so the compiler keeps the !$acc parallel region clean
         ! (no absent-optional reference inside the parallel body).
         acc = 0.0_wp
         !$acc parallel loop collapse(3) reduction(max:acc) &
         !$acc&         present(u_face, v_face, h_layer, idxT, idyT)
         do k = 1, nz
            do j = j_lo, j_hi
               do i = i_lo, i_hi
                  if (h_layer(i, j, k) > vanish_tol) then
                     acc = max(acc, cfl_cell_value( &
                               u_face(i, j, k), u_face(i + 1, j, k), &
                               v_face(i, j, k), v_face(i, j + 1, k), &
                               idxT(i, j), idyT(i, j), dt))
                  end if
               end do
            end do
         end do
      else
         ! Un-gated path: all layers contribute.  Byte-identical to the
         ! pre-Phase-3 implementation.
         acc = 0.0_wp
         !$acc parallel loop collapse(3) reduction(max:acc) &
         !$acc&         present(u_face, v_face, idxT, idyT)
         do k = 1, nz
            do j = j_lo, j_hi
               do i = i_lo, i_hi
                  acc = max(acc, cfl_cell_value( &
                            u_face(i, j, k), u_face(i + 1, j, k), &
                            v_face(i, j, k), v_face(i, j + 1, k), &
                            idxT(i, j), idyT(i, j), dt))
               end do
            end do
         end do
      end if
      max_cfl = acc
   end function compute_max_cfl