compute_h_face_upstream Subroutine

   pure subroutine compute_h_face_upstream(grid, bt_work, ms, metrics)
      !! Per-face upstream column-sum thickness `h_face_up_x/y(I,j) =
      !! Σ_k h_layer(I_upstream,j,k)` used by the BT chain when
      !! `use_upstream_h_face = .true.`. First-order upwind pick by face-velocity
      !! sign (sampled at the top of the outer step). Wall faces use the single
      !! available cell. No-op when the knob is off.
      !!
      !! Under `&vcoord_nml zfixed_closed_faces` the sum is over the OPEN
      !! column instead — see `metrics` and `h_face_upstream_open_impl`.
      type(hgrid_t), intent(in) :: grid
      type(barotropic_workstate_t), intent(inout) :: bt_work
      type(multilayer_state_t), intent(in) :: ms
      type(ocean_metrics_t), intent(in) :: metrics
         !! REQUIRED, for the reason `derive_bt_from_layers` gives: an
         !! optional dummy that silently selects the full-column branch is
         !! how the `pred_corr` Coriolis-reference defect shipped.
         !!
         !! `metrics%use_closed_faces = .false.` (the default) ⇒ the
         !! ORIGINAL full-column loops below run, textually unchanged
         !! (byte-identical), and neither `open_*` nor `dy_cu_bt` is
         !! named.  `.true.` ⇒ the open-column builder
         !! `h_face_upstream_open_impl`, whose docstring derives why the
         !! full-column sum is not merely imprecise there but makes the
         !! barotropic transport disagree with the renormalised layer
         !! transports by a factor up to ~2 at a staircase face.

      integer :: i, j, k, nx, ny, nz, nx_face, ny_face
      real(wp) :: h_sum, h_k

      if (.not. bt_work%use_upstream_h_face) return

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      nx_face = size(bt_work%h_face_up_x, 1)
      ny_face = size(bt_work%h_face_up_y, 2)

      if (metrics%use_closed_faces) then
         ! Two calls, one per porous state: with porous OFF the
         ! `por_face_area_*` arrays are the `(1,1,1)` placeholder and must
         ! NOT reach the callee's explicit-shape dummy (nvfortran builds
         ! the `do concurrent` data clause from the loop bounds, so a
         ! placeholder aborts under `mem:separate` even when the branch
         ! that indexes it is never taken).  The mask itself is the inert
         ! stand-in — right shape, already mapped, `intent(in)` at both
         ! dummies — the same device the `ocean_porous_refresh` call uses.
         if (metrics%use_porous) then
            call h_face_upstream_open_impl(nx, ny, nz, .true., ms%h_layer, &
                                           ms%u_face_x_layer, ms%v_face_y_layer, &
                                           metrics%open_u, metrics%open_v, &
                                           metrics%por_face_area_u, &
                                           metrics%por_face_area_v, &
                                           metrics%dy_cu, metrics%dx_cv, &
                                           metrics%dy_cu_bt, metrics%dx_cv_bt, &
                                           bt_work%h_face_up_x, bt_work%h_face_up_y)
         else
            call h_face_upstream_open_impl(nx, ny, nz, .false., ms%h_layer, &
                                           ms%u_face_x_layer, ms%v_face_y_layer, &
                                           metrics%open_u, metrics%open_v, &
                                           metrics%open_u, metrics%open_v, &
                                           metrics%dy_cu, metrics%dx_cv, &
                                           metrics%dy_cu_bt, metrics%dx_cv_bt, &
                                           bt_work%h_face_up_x, bt_work%h_face_up_y)
         end if
         return
      end if

      ! East-face: upstream pick from u_face_x_layer sign.
      do concurrent(j=1:ny, i=1:nx_face) local(k, h_sum, h_k)
         h_sum = 0.0_wp
         do k = 1, nz
            if (i == 1) then
               h_k = ms%h_layer(1, j, k)
            else if (i == nx_face) then
               h_k = ms%h_layer(nx, j, k)
            else if (ms%u_face_x_layer(i, j, k) >= 0.0_wp) then
               h_k = ms%h_layer(i - 1, j, k)
            else
               h_k = ms%h_layer(i, j, k)
            end if
            h_sum = h_sum + h_k
         end do
         bt_work%h_face_up_x(i, j) = h_sum
      end do

      ! North-face: upstream pick from v_face_y_layer sign.
      do concurrent(j=1:ny_face, i=1:nx) local(k, h_sum, h_k)
         h_sum = 0.0_wp
         do k = 1, nz
            if (j == 1) then
               h_k = ms%h_layer(i, 1, k)
            else if (j == ny_face) then
               h_k = ms%h_layer(i, ny, k)
            else if (ms%v_face_y_layer(i, j, k) >= 0.0_wp) then
               h_k = ms%h_layer(i, j - 1, k)
            else
               h_k = ms%h_layer(i, j, k)
            end if
            h_sum = h_sum + h_k
         end do
         bt_work%h_face_up_y(i, j) = h_sum
      end do
   end subroutine compute_h_face_upstream