porous_fill_stats_resolved Subroutine

   pure subroutine porous_fill_stats_resolved(nx, ny, b, wet_t, &
                                              dmin_u, dmax_u, davg_u, &
                                              dmin_v, dmax_v, davg_v)
      !! Fill the per-face `d_min` / `d_max` / `d_avg` from the RESOLVED
      !! bottom elevation, sampled at three points ALONG each face: the
      !! two end corners and the midpoint.  A corner sample is the mean
      !! of the four cells around it, the midpoint the mean of the two
      !! cells the face separates — all three are genuine points on the
      !! resolved seafloor along the face segment.
      !!
      !! WET GATING.  A corner sample is used ONLY when all four cells in
      !! its stencil are wet; otherwise it falls back to the two-cell face
      !! midpoint.  Without this a wet-wet face with one diagonal LAND
      !! neighbour picks the land elevation up into `d_max` and blocks a
      !! face the grid fully resolves as open (measured: a single 4000 m
      !! land diagonal on a 4000 m column gives ~37% spurious blockage on
      !! the deepest layers).  When both corners are gated out all three
      !! samples coincide, the face statistic is degenerate, and
      !! `porous_update_face_areas` leaves it fully open.
      !!
      !! This is a PROXY for true subgrid statistics (see the module
      !! caveat): it resolves along-face slope but is blind to structure
      !! below the grid scale.  A flat seafloor — or ANY seafloor uniform
      !! along the face — gives `d_min = d_max = d_avg`, the degenerate
      !! case the kernel leaves fully open.  That is the bit-identity
      !! property flat-bottom configurations rely on, and the reason a
      !! ridge running parallel to the face does not become a wall.
      !!
      !! The outermost face columns (`i = 1`, `i = nx+1` for u;
      !! `j = 1`, `j = ny+1` for v) are left at zero: they have no adjacent
      !! cell pair, `porous_update_face_areas` skips them, and their
      !! metric width is already zero.
      !!
      !! Setup-phase host code: plain sequential loops, NOT
      !! `do concurrent` (this runs before `enter_data`, where a DC loop
      !! would round-trip the unmapped arrays through the device).
      integer, intent(in) :: nx, ny
         !! `grid%nx_total`, `grid%ny_total`.
      real(wp), intent(in) :: b(nx, ny)
         !! Bottom elevation at cell centres (m, positive up).
      real(wp), intent(in) :: wet_t(nx, ny)
         !! T-cell wet (1) / land (0) mask (`ocean_metrics_t%wet_T`).
         !! All-ones on a run with no land, where every corner sample is
         !! kept and the statistic is the un-gated one.
      real(wp), intent(out) :: dmin_u(nx + 1, ny), dmax_u(nx + 1, ny), davg_u(nx + 1, ny)
         !! u-face along-face deepest / shallowest / mean height (m).
      real(wp), intent(out) :: dmin_v(nx, ny + 1), dmax_v(nx, ny + 1), davg_v(nx, ny + 1)
         !! v-face along-face deepest / shallowest / mean height (m).

      integer :: i, j
      real(wp) :: s_lo, s_mid, s_hi

      ! ---- u-faces: face i separates cells (i-1,j) and (i,j) ----
      ! Samples run south -> north along the face segment.
      dmin_u = 0.0_wp
      dmax_u = 0.0_wp
      davg_u = 0.0_wp
      do j = 2, ny - 1
         do i = 2, nx
            s_mid = 0.5_wp*(b(i - 1, j) + b(i, j))
            if (corner_is_wet(wet_t(i - 1, j - 1), wet_t(i, j - 1), &
                              wet_t(i - 1, j), wet_t(i, j))) then
               s_lo = 0.25_wp*(b(i - 1, j - 1) + b(i, j - 1) + b(i - 1, j) + b(i, j))
            else
               s_lo = s_mid
            end if
            if (corner_is_wet(wet_t(i - 1, j + 1), wet_t(i, j + 1), &
                              wet_t(i - 1, j), wet_t(i, j))) then
               s_hi = 0.25_wp*(b(i - 1, j + 1) + b(i, j + 1) + b(i - 1, j) + b(i, j))
            else
               s_hi = s_mid
            end if
            dmin_u(i, j) = min(s_lo, s_mid, s_hi)
            dmax_u(i, j) = max(s_lo, s_mid, s_hi)
            ! Simpson (1,4,1)/6 — the consistent along-face MEAN of three
            ! evenly spaced samples.  A plain (1,1,1)/3 average would
            ! under-weight the midpoint and bias `d_avg` toward the ends.
            !
            ! Bracketed because the invariant `d_min <= d_avg <= d_max` is
            ! exact only in exact arithmetic: the weights are a convex
            ! combination, but three near-equal samples can round the sum
            ! an ulp outside the range, and `m = (d_avg-d_min)/(d_max-d_min)`
            ! must stay in [0,1].  The bracket moves nothing else.
            davg_u(i, j) = min(dmax_u(i, j), max(dmin_u(i, j), &
                                                 (s_lo + 4.0_wp*s_mid + s_hi)/6.0_wp))
         end do
      end do
      ! Outermost rows have no cross-face neighbour for the corner
      ! samples: fall back to the two-cell midpoint (a flat face, so the
      ! fit degenerates to the open/closed step and blocks nothing extra).
      do i = 2, nx
         s_mid = 0.5_wp*(b(i - 1, 1) + b(i, 1))
         dmin_u(i, 1) = s_mid
         dmax_u(i, 1) = s_mid
         davg_u(i, 1) = s_mid
         s_mid = 0.5_wp*(b(i - 1, ny) + b(i, ny))
         dmin_u(i, ny) = s_mid
         dmax_u(i, ny) = s_mid
         davg_u(i, ny) = s_mid
      end do

      ! ---- v-faces: face j separates cells (i,j-1) and (i,j) ----
      ! Samples run west -> east along the face segment.
      dmin_v = 0.0_wp
      dmax_v = 0.0_wp
      davg_v = 0.0_wp
      do j = 2, ny
         do i = 2, nx - 1
            s_mid = 0.5_wp*(b(i, j - 1) + b(i, j))
            if (corner_is_wet(wet_t(i - 1, j - 1), wet_t(i - 1, j), &
                              wet_t(i, j - 1), wet_t(i, j))) then
               s_lo = 0.25_wp*(b(i - 1, j - 1) + b(i - 1, j) + b(i, j - 1) + b(i, j))
            else
               s_lo = s_mid
            end if
            if (corner_is_wet(wet_t(i + 1, j - 1), wet_t(i + 1, j), &
                              wet_t(i, j - 1), wet_t(i, j))) then
               s_hi = 0.25_wp*(b(i + 1, j - 1) + b(i + 1, j) + b(i, j - 1) + b(i, j))
            else
               s_hi = s_mid
            end if
            dmin_v(i, j) = min(s_lo, s_mid, s_hi)
            dmax_v(i, j) = max(s_lo, s_mid, s_hi)
            ! Bracketed — see the zonal twin.
            davg_v(i, j) = min(dmax_v(i, j), max(dmin_v(i, j), &
                                                 (s_lo + 4.0_wp*s_mid + s_hi)/6.0_wp))
         end do
      end do
      do j = 2, ny
         s_mid = 0.5_wp*(b(1, j - 1) + b(1, j))
         dmin_v(1, j) = s_mid
         dmax_v(1, j) = s_mid
         davg_v(1, j) = s_mid
         s_mid = 0.5_wp*(b(nx, j - 1) + b(nx, j))
         dmin_v(nx, j) = s_mid
         dmax_v(nx, j) = s_mid
         davg_v(nx, j) = s_mid
      end do
   end subroutine porous_fill_stats_resolved