ocean_sigma_stiffness_worst Subroutine

   pure subroutine ocean_sigma_stiffness_worst(nx, ny, i0, i1, j0, j1, column, wet, &
                                               rx0_max, i_at, j_at, is_x, &
                                               h_thin, h_thick, n_over, n_face)
      !! Worst (largest) `ocean_sigma_stiffness` over every face joining
      !! two WET columns inside `[i0,i1] x [j0,j1]`, with its location,
      !! its two column thicknesses, and how many faces are over
      !! `SIGMA_STIFFNESS_LIMIT`.
      !!
      !! Only INTERIOR-to-INTERIOR faces are scanned (the loops stop one
      !! short of `i1`/`j1`), so the number never depends on what a ghost
      !! ring happens to hold — which is what makes it the same under any
      !! decomposition and safe to quote in a configure message.
      !!
      !! A land column is excluded by `wet`, not by a thickness test: a
      !! land T-cell holds `h_layer = H_VANISHED` per the land-state
      !! contract, so its "column" is `nz*H_VANISHED` and would otherwise
      !! read as a near-vanishing neighbour at every coastline and make
      !! `rx0 -> 1` everywhere. A coastline is a WALL, not a stiff face:
      !! the metrics are zeroed there and no pressure gradient is taken.
      integer, intent(in) :: nx
         !! First dimension of `column`/`wet` (ghosts included).
      integer, intent(in) :: ny
         !! Second dimension.
      integer, intent(in) :: i0
         !! First physical index in x.
      integer, intent(in) :: i1
         !! Last physical index in x.
      integer, intent(in) :: j0
         !! First physical index in y.
      integer, intent(in) :: j1
         !! Last physical index in y.
      real(wp), intent(in) :: column(nx, ny)
         !! Column thickness (m) the vertical coordinate divides — the
         !! WATER column `b - z_draft` under an ice shelf.
      real(wp), intent(in) :: wet(nx, ny)
         !! Static wet (1) / land (0) T-cell mask.
      real(wp), intent(out) :: rx0_max
         !! Largest stiffness found; `0` if no wet-wet face exists.
      integer, intent(out) :: i_at
         !! `i` of the thin side of the worst face.
      integer, intent(out) :: j_at
         !! `j` of the thin side of the worst face.
      logical, intent(out) :: is_x
         !! `.true.` if the worst face is an x (east) face.
      real(wp), intent(out) :: h_thin
         !! Thinner column of the worst face (m).
      real(wp), intent(out) :: h_thick
         !! Thicker column of the worst face (m).
      integer, intent(out) :: n_over
         !! Wet-wet faces with `rx0 > SIGMA_STIFFNESS_LIMIT`.
      integer, intent(out) :: n_face
         !! Wet-wet faces scanned (the denominator for `n_over`).
      integer :: i, j
      real(wp) :: rx0, ha, hb

      rx0_max = 0.0_wp
      i_at = i0
      j_at = j0
      is_x = .true.
      h_thin = 0.0_wp
      h_thick = 0.0_wp
      n_over = 0
      n_face = 0

      do j = j0, j1
         do i = i0, i1
            if (wet(i, j) <= 0.5_wp) cycle
            ha = column(i, j)
            if (ha <= 0.0_wp) cycle
            if (i < i1) then
               if (wet(i + 1, j) > 0.5_wp) then
                  hb = column(i + 1, j)
                  if (hb > 0.0_wp) then
                     rx0 = ocean_sigma_stiffness(ha, hb)
                     n_face = n_face + 1
                     if (rx0 > SIGMA_STIFFNESS_LIMIT) n_over = n_over + 1
                     if (rx0 > rx0_max) then
                        rx0_max = rx0
                        i_at = i
                        j_at = j
                        is_x = .true.
                        h_thin = min(ha, hb)
                        h_thick = max(ha, hb)
                     end if
                  end if
               end if
            end if
            if (j < j1) then
               if (wet(i, j + 1) > 0.5_wp) then
                  hb = column(i, j + 1)
                  if (hb > 0.0_wp) then
                     rx0 = ocean_sigma_stiffness(ha, hb)
                     n_face = n_face + 1
                     if (rx0 > SIGMA_STIFFNESS_LIMIT) n_over = n_over + 1
                     if (rx0 > rx0_max) then
                        rx0_max = rx0
                        i_at = i
                        j_at = j
                        is_x = .false.
                        h_thin = min(ha, hb)
                        h_thick = max(ha, hb)
                     end if
                  end if
               end if
            end if
         end do
      end do
   end subroutine ocean_sigma_stiffness_worst