porous_update_face_areas Subroutine

   pure subroutine porous_update_face_areas(nx, ny, nz, interp, mask_depth, &
                                            b, h_layer, &
                                            dmin_u, dmax_u, davg_u, &
                                            dmin_v, dmax_v, davg_v, &
                                            dy_cu, dx_cv, &
                                            por_u, por_v, dy_cu_bt, dx_cv_bt)
      !! Recompute the layer-averaged open-area fractions from the
      !! CURRENT layer thicknesses.  Interface-height dependent, so this
      !! runs once per RK2 stage on the device.
      !!
      !! Column recurrence (bottom-up, `k=1` is the bed): the face
      !! interface height is interpolated from the two adjacent columns'
      !! running interface heights, the cumulative open area is evaluated
      !! there, and the layer fraction is the increment over the layer
      !! divided by the layer's face thickness.  Everything is a running
      !! scalar — no per-column workspace, so there is nothing to
      !! allocate and nothing to map.
      !!
      !! Three faces are left FULLY OPEN (`por = 1`, BT width untouched):
      !!   * array-edge faces (no adjacent cell pair; metric width is
      !!     already zero),
      !!   * faces whose mean along-face height is at or above
      !!     `mask_depth` — MOM6's `PORBAR_MASKING_DEPTH` gate; porous
      !!     barriers are a deep-sill parameterization and should not
      !!     narrow shelf faces the grid already resolves,
      !!   * DEGENERATE faces, `d_max <= d_min`.  The fit's limit there is
      !!     a step at `d_min`, which on a seafloor uniform ALONG the face
      !!     puts a hard wall at the two-cell mean height — an artifact of
      !!     the resolved-bathymetry proxy, not subgrid information (see
      !!     the module docstring).  Leaving them open also makes a flat
      !!     basin a literal no-op, which is what the bit-identity
      !!     argument for `source="resolved"` rests on.
      !!
      !! A wholly VANISHED column (top interface within `H_VANISHED` of
      !! the bed) gets `por = 0` on every layer AND a zero BT width, so
      !! the barotropic mode and the layers agree that nothing can move
      !! through it.  Handing the BT the un-narrowed width there would let
      !! it transport at full width across a face every layer has blocked.
      integer, intent(in) :: nx, ny, nz
         !! `grid%nx_total`, `grid%ny_total`, number of layers.
      integer, intent(in) :: interp
         !! Interface-at-velocity-point rule, a `POROUS_ETA_*` value.
      real(wp), intent(in) :: mask_depth
         !! Gate height (m, positive up, `<= 0`): faces with
         !! `d_avg >= mask_depth` stay fully open.
      real(wp), intent(in) :: b(nx, ny)
         !! Bottom elevation at cell centres (m, positive up).
      real(wp), intent(in) :: h_layer(nx, ny, nz)
         !! Layer thicknesses (m), `k=1` the bed layer.
      real(wp), intent(in) :: 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(in) :: dmin_v(nx, ny + 1), dmax_v(nx, ny + 1), davg_v(nx, ny + 1)
         !! v-face along-face deepest / shallowest / mean height (m).
      real(wp), intent(in) :: dy_cu(nx + 1, ny)
         !! Un-narrowed open u-face width for transport (m).
      real(wp), intent(in) :: dx_cv(nx, ny + 1)
         !! Un-narrowed open v-face width for transport (m).
      real(wp), intent(out) :: por_u(nx + 1, ny, nz)
         !! u-face layer-averaged open-area fraction (nondim, `[0,1]`).
      real(wp), intent(out) :: por_v(nx, ny + 1, nz)
         !! v-face layer-averaged open-area fraction (nondim, `[0,1]`).
      real(wp), intent(out) :: dy_cu_bt(nx + 1, ny)
         !! u-face width the BAROTROPIC substep transports on (m):
         !! `dy_cu` scaled by the COLUMN-INTEGRATED open fraction.
      real(wp), intent(out) :: dx_cv_bt(nx, ny + 1)
         !! v-face twin (m).

      integer :: i, j, k
      real(wp) :: z_a, z_b, e_lo, e_hi, a_lo, a_hi, dz
      real(wp) :: e_bed, a_bed, por_col

      ! ---- u-faces ----
      do concurrent(j=1:ny, i=1:nx + 1) &
         local(k, z_a, z_b, e_lo, e_hi, a_lo, a_hi, dz, e_bed, a_bed, por_col)
         if (i < 2 .or. i > nx) then
            ! Array-edge faces have no adjacent cell pair; their metric
            ! width is already zero, so leave them fully open.
            do k = 1, nz
               por_u(i, j, k) = 1.0_wp
            end do
            dy_cu_bt(i, j) = dy_cu(i, j)
         else if (davg_u(i, j) >= mask_depth) then
            do k = 1, nz
               por_u(i, j, k) = 1.0_wp
            end do
            dy_cu_bt(i, j) = dy_cu(i, j)
         else if (dmax_u(i, j) <= dmin_u(i, j)) then
            ! Degenerate face: the along-face statistic carries no
            ! information, so block nothing (see the routine docstring).
            do k = 1, nz
               por_u(i, j, k) = 1.0_wp
            end do
            dy_cu_bt(i, j) = dy_cu(i, j)
         else
            z_a = b(i - 1, j)
            z_b = b(i, j)
            e_lo = porous_eta_face(z_a, z_b, interp)
            a_lo = porous_cum_area(dmin_u(i, j), dmax_u(i, j), davg_u(i, j), e_lo)
            e_bed = e_lo
            a_bed = a_lo
            do k = 1, nz
               z_a = z_a + h_layer(i - 1, j, k)
               z_b = z_b + h_layer(i, j, k)
               e_hi = porous_eta_face(z_a, z_b, interp)
               a_hi = porous_cum_area(dmin_u(i, j), dmax_u(i, j), davg_u(i, j), e_hi)
               dz = e_hi - e_lo
               if (dz > H_VANISHED) then
                  por_u(i, j, k) = clamp_fraction((a_hi - a_lo)/dz)
               else
                  ! Vanished layer at this face: no open area to speak of.
                  por_u(i, j, k) = 0.0_wp
               end if
               e_lo = e_hi
               a_lo = a_hi
            end do
            ! Column-integrated open fraction: the SAME layer-average
            ! formula applied over the whole water column, which is
            ! identically the thickness-weighted mean of the per-layer
            ! fractions.  This is what makes the barotropic transport see
            ! the barrier (MOM6 folds the equivalent quantity into its
            ! BT_cont face area).  A wholly VANISHED column has every
            ! layer blocked above, so the BT width goes to zero too — the
            ! two modes must agree about a face nothing can pass through.
            if (e_lo - e_bed > H_VANISHED) then
               por_col = clamp_fraction((a_lo - a_bed)/(e_lo - e_bed))
               dy_cu_bt(i, j) = dy_cu(i, j)*por_col
            else
               dy_cu_bt(i, j) = 0.0_wp
            end if
         end if
      end do

      ! ---- v-faces ----
      do concurrent(j=1:ny + 1, i=1:nx) &
         local(k, z_a, z_b, e_lo, e_hi, a_lo, a_hi, dz, e_bed, a_bed, por_col)
         if (j < 2 .or. j > ny) then
            do k = 1, nz
               por_v(i, j, k) = 1.0_wp
            end do
            dx_cv_bt(i, j) = dx_cv(i, j)
         else if (davg_v(i, j) >= mask_depth) then
            do k = 1, nz
               por_v(i, j, k) = 1.0_wp
            end do
            dx_cv_bt(i, j) = dx_cv(i, j)
         else if (dmax_v(i, j) <= dmin_v(i, j)) then
            ! Degenerate face — see the zonal twin.
            do k = 1, nz
               por_v(i, j, k) = 1.0_wp
            end do
            dx_cv_bt(i, j) = dx_cv(i, j)
         else
            z_a = b(i, j - 1)
            z_b = b(i, j)
            e_lo = porous_eta_face(z_a, z_b, interp)
            a_lo = porous_cum_area(dmin_v(i, j), dmax_v(i, j), davg_v(i, j), e_lo)
            e_bed = e_lo
            a_bed = a_lo
            do k = 1, nz
               z_a = z_a + h_layer(i, j - 1, k)
               z_b = z_b + h_layer(i, j, k)
               e_hi = porous_eta_face(z_a, z_b, interp)
               a_hi = porous_cum_area(dmin_v(i, j), dmax_v(i, j), davg_v(i, j), e_hi)
               dz = e_hi - e_lo
               if (dz > H_VANISHED) then
                  por_v(i, j, k) = clamp_fraction((a_hi - a_lo)/dz)
               else
                  por_v(i, j, k) = 0.0_wp
               end if
               e_lo = e_hi
               a_lo = a_hi
            end do
            ! Column-integrated open fraction — see the zonal twin.
            if (e_lo - e_bed > H_VANISHED) then
               por_col = clamp_fraction((a_lo - a_bed)/(e_lo - e_bed))
               dx_cv_bt(i, j) = dx_cv(i, j)*por_col
            else
               dx_cv_bt(i, j) = 0.0_wp
            end if
         end if
      end do
   end subroutine porous_update_face_areas