face_depth_mean_rem_u Subroutine

   pure subroutine face_depth_mean_rem_u(grid, F_3d, h_layer, rem, F_mean_2d, nz, metrics, n_inner, &
                                         href, scheme)
      !! `face_depth_mean_u` with MOM6 `wt_u` weighting (`&ocean_bt_nml
      !! forcing_visc_rem`): the weight is
      !! `h_face·visc_rem(k)` instead of `h_face`, so layers the implicit
      !! vertical-friction solve will immediately damp (grounded sliver
      !! stacks under the BBL glue, visc_rem → 0) contribute nothing to
      !! the barotropic forcing.  Without this the spurious grounded-layer
      !! PGF's depth-mean drives the fast loop ballistically even after
      !! the layer velocities themselves are glued (PGF_BUG.md §9).
      !! Denominator falls back to zero-output on an all-remnant-zero
      !! column (the substep should not force an immobilized column).
      !!
      !! `rem` is run through MOM6's exact `wt_u` floor before it weights
      !! anything: `vr = min(rem, 1)`,
      !! `vr = max(vr, 1 - 0.5·Instep/(vr + subroundoff))`,
      !! `vr = max(vr, 0)`, `Instep = 1/n_inner` — NOT roundabout's old
      !! plain `[0,1]` clamp, which let a near-zero `visc_rem` on a
      !! many-substep column weight the forcing far closer to zero than
      !! MOM6 ever lets it (the floor's whole job is to keep the
      !! `Instep`-th root finite on exactly this kind of thin cell; see
      !! the plan's "Thin-cell floors" risk note).  `ieee_is_finite`-
      !! guarded per the CLAUDE.md NaN-clamp gotcha: a non-finite `rem`
      !! passes through unmasked rather than being laundered to 0 or 1.
      type(hgrid_t), intent(in) :: grid
      ! assumed-shape-ok: face arrays have shape (nx+1,ny,nz); a single (nx,ny,nz)
      ! explicit-shape triplet would mis-bound the face axis.
      real(wp), intent(in) :: F_3d(:, :, :)
      real(wp), intent(in) :: h_layer(:, :, :)  ! assumed-shape-ok: face-sized array; size() derives loop bounds
      real(wp), intent(in) :: rem(:, :, :)      ! assumed-shape-ok: face-sized array; size() derives loop bounds
      real(wp), intent(out) :: F_mean_2d(:, :)  ! assumed-shape-ok: face-sized array; size() derives loop bounds
      integer, intent(in) :: nz
      type(ocean_metrics_t), intent(in) :: metrics
         !! REQUIRED (see `face_depth_mean_u`): the weight becomes
         !! `h_face·visc_rem·open`.  It must match
         !! `derive_bt_from_layers` and `apply_bt_correction` or the
         !! `dt·F_bt` the fold subtracts is not what the fast loop
         !! integrated.  Note a CLOSED layer's `visc_rem` is ~1, not 0 —
         !! the closed-face vdiff decoupling leaves it uncoupled, so
         !! `visc_rem` alone does NOT stand in for the mask here.
         !! Knob off ⇒ the ORIGINAL loop, byte-identical.
      integer, intent(in) :: n_inner
         !! Barotropic substep count (MOM6 `nstep`); `Instep = 1/n_inner`
         !! in the `wt_u` floor.  `max(n_inner, 1)` guards the unsplit
         !! (`n_inner = 0`) configuration.
      real(wp), intent(in) :: href(:, :)  ! assumed-shape-ok: cell-centred (nx,ny); see face_depth_mean_u
      integer, intent(in) :: scheme
         !! A `FRHAT_*` constant.  See `face_depth_mean_u`.
      integer :: i, j, k, nu, ny, nx_cells, il, ir
      real(wp) :: h_face, wt, num, denom, vr, instep, href_l, href_r, e_prev
      logical :: use_open
      integer :: eff_scheme

      nu = size(F_3d, 1)
      ny = size(F_3d, 2)
      nx_cells = grid%nx_total

      use_open = metrics%use_closed_faces
      ! FRHAT_HYBRID only under closed faces -- see `derive_bt_from_layers`'
      ! matching comment.
      eff_scheme = merge(scheme, FRHAT_ARITHMETIC, use_open)
      instep = 1.0_wp/real(max(n_inner, 1), wp)

      if (use_open) then
         do concurrent(j=1:ny, i=1:nu) &
            local(k, il, ir, href_l, href_r, e_prev, h_face, wt, num, denom, vr)
            il = max(1, i - 1)
            ir = min(nx_cells, i)
            href_l = href(il, j)
            href_r = href(ir, j)
            e_prev = -0.5_wp*(href_l + href_r)
            num = 0.0_wp
            denom = 0.0_wp
            do k = 1, nz
               call frhat_h_face_step(h_layer(il, j, k), h_layer(ir, j, k), &
                                      href_l, href_r, eff_scheme, e_prev, h_face)
               if (ieee_is_finite(rem(i, j, k))) then
                  ! Clamp to [0, 1] BEFORE MOM6's floor: MOM6's remnant is
                  ! non-negative by construction, but a round-off -1e-17 here
                  ! would make 1 - 0.5*Instep/vr ~ +5e15 and survive max(., 0).
                  vr = max(min(rem(i, j, k), 1.0_wp), 0.0_wp)
                  vr = max(vr, 1.0_wp - 0.5_wp*instep/(vr + VISC_REM_SUBROUNDOFF))
                  vr = max(vr, 0.0_wp)
               else
                  vr = rem(i, j, k)
               end if
               wt = metrics%open_u(i, j, k)*h_face*vr
               num = num + F_3d(i, j, k)*wt
               denom = denom + wt
            end do
            if (denom > 0.0_wp) then
               F_mean_2d(i, j) = num/denom
            else
               F_mean_2d(i, j) = 0.0_wp
            end if
         end do
      else
         do concurrent(j=1:ny, i=1:nu) &
            local(k, il, ir, href_l, href_r, e_prev, h_face, wt, num, denom, vr)
            il = max(1, i - 1)
            ir = min(nx_cells, i)
            href_l = href(il, j)
            href_r = href(ir, j)
            e_prev = -0.5_wp*(href_l + href_r)
            num = 0.0_wp
            denom = 0.0_wp
            do k = 1, nz
               call frhat_h_face_step(h_layer(il, j, k), h_layer(ir, j, k), &
                                      href_l, href_r, eff_scheme, e_prev, h_face)
               if (ieee_is_finite(rem(i, j, k))) then
                  ! Clamp to [0, 1] BEFORE MOM6's floor: MOM6's remnant is
                  ! non-negative by construction, but a round-off -1e-17 here
                  ! would make 1 - 0.5*Instep/vr ~ +5e15 and survive max(., 0).
                  vr = max(min(rem(i, j, k), 1.0_wp), 0.0_wp)
                  vr = max(vr, 1.0_wp - 0.5_wp*instep/(vr + VISC_REM_SUBROUNDOFF))
                  vr = max(vr, 0.0_wp)
               else
                  vr = rem(i, j, k)
               end if
               wt = h_face*vr
               num = num + F_3d(i, j, k)*wt
               denom = denom + wt
            end do
            if (denom > 0.0_wp) then
               F_mean_2d(i, j) = num/denom
            else
               F_mean_2d(i, j) = 0.0_wp
            end if
         end do
      end if
   end subroutine face_depth_mean_rem_u