face_depth_mean_u Subroutine

   pure subroutine face_depth_mean_u(grid, F_3d, h_layer, F_mean_2d, nz, metrics, href, scheme)
      !! Depth-average a u-face 3D field, weighted by the per-layer face
      !! thickness `frhat_h_face_step` returns (`&ocean_bt_nml
      !! frhat_scheme`; `FRHAT_ARITHMETIC` = the plain mean of the two
      !! abutting cell columns' `h_layer` values, today's default; the
      !! argument is otherwise unchanged).  Writes to a 2D field at the
      !! same u-face shape.  Wall faces (i=1, nx+1) fall back to the
      !! single available cell (`frhat_h_face_step` degenerates there for
      !! either scheme — see its docstring).
      type(hgrid_t), intent(in) :: grid
      ! assumed-shape-ok: face arrays have shape (nx+1,ny,nz) / (nx,ny+1,nz);
      ! a single (nx,ny,nz) explicit-shape triplet would mis-bound the face axis.
      ! size() is used to derive loop bounds from the actual face dimension.
      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(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.  Under `&vcoord_nml zfixed_closed_faces` every
         !! depth mean in the split chain MUST use the same weights
         !! `derive_bt_from_layers` and `apply_bt_correction` use —
         !! `h_face·open` — or the `dt·F_bt` the fold subtracts back out
         !! is not the quantity the fast loop integrated, and the
         !! difference survives as a permanent per-face bias.
         !!
         !! Not optional, deliberately: this dummy WAS optional, and the
         !! one call site that omitted it (`set_cor_ref_velocity`) turned
         !! the `pred_corr` Coriolis reference into `φ·ū_open` and pumped
         !! the ISOMIP+ Ocean0 barotropic mode by ×22.7 in energy over 30
         !! days.  A caller that has no mask still has an
         !! `ocean_metrics_t` to hand; `use_closed_faces = .false.` ⇒ the
         !! ORIGINAL loop runs, textually unchanged (byte-identical).
      real(wp), intent(in) :: href(:, :)  ! assumed-shape-ok: cell-centred (nx,ny); size() unused, indexed by column
         !! Reference column depth (m), `barotropic_workstate_t%bt_H_ref` —
         !! REQUIRED for the same reason `metrics` is: every call site has
         !! a `bt_work` to hand, and an omitted/mismatched `href` would
         !! silently change `FRHAT_HYBRID`'s answer. Unused under
         !! `FRHAT_ARITHMETIC` but still dereferenced (no branch to skip
         !! it), so it must be a real, fully-sized array on every call,
         !! never a placeholder.
      integer, intent(in) :: scheme
         !! A `FRHAT_*` constant (`rdb_constants`; `&ocean_bt_nml
         !! frhat_scheme`, via `bt_work%frhat_scheme`).
      integer :: i, j, k, nu, ny, nx_cells, il, ir
      real(wp) :: h_face, num, denom, 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; off that family this is FRHAT_ARITHMETIC
      ! regardless of what the caller passed as `scheme`.
      eff_scheme = merge(scheme, FRHAT_ARITHMETIC, use_open)

      if (use_open) then
         do concurrent(j=1:ny, i=1:nu) local(k, il, ir, href_l, href_r, e_prev, h_face, num, denom)
            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)
               h_face = h_face*metrics%open_u(i, j, k)
               num = num + F_3d(i, j, k)*h_face
               denom = denom + h_face
            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, num, denom)
            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)
               num = num + F_3d(i, j, k)*h_face
               denom = denom + h_face
            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_u