face_depth_mean_rem_v Subroutine

   pure subroutine face_depth_mean_rem_v(grid, F_3d, h_layer, rem, F_mean_2d, nz, metrics, n_inner, &
                                         href, scheme)
      !! Symmetric v-face counterpart of `face_depth_mean_rem_u` — same
      !! MOM6 `wt_u` floor, `ieee_is_finite`-guarded the same way.
      type(hgrid_t), intent(in) :: grid
      ! assumed-shape-ok: face arrays have shape (nx,ny+1,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`); see `face_depth_mean_rem_u`.
      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, nx, nv, ny_cells, jl, jr
      real(wp) :: h_face, wt, num, denom, vr, instep, href_l, href_r, e_prev
      logical :: use_open
      integer :: eff_scheme

      nx = size(F_3d, 1)
      nv = size(F_3d, 2)
      ny_cells = grid%ny_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:nv, i=1:nx) &
            local(k, jl, jr, href_l, href_r, e_prev, h_face, wt, num, denom, vr)
            jl = max(1, j - 1)
            jr = min(ny_cells, j)
            href_l = href(i, jl)
            href_r = href(i, jr)
            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(i, jl, k), h_layer(i, jr, 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_v(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:nv, i=1:nx) &
            local(k, jl, jr, href_l, href_r, e_prev, h_face, wt, num, denom, vr)
            jl = max(1, j - 1)
            jr = min(ny_cells, j)
            href_l = href(i, jl)
            href_r = href(i, jr)
            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(i, jl, k), h_layer(i, jr, 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_v