derive_bt_from_layers Subroutine

   pure subroutine derive_bt_from_layers(grid, bt_work, ms, metrics)
      !! Populate `bt_eta`, `bt_ubt`, `bt_vbt` from the current multilayer
      !! state. `bt_H_ref` must already be set.
      !!   bt_eta = Σ_k h_layer − H_ref;  bt_ubt = Σ_k(u·h_face)/Σ_k h_face.
      !! Face thickness averages the two abutting columns (wall faces use the
      !! single cell). With `use_upstream_h_face`, the interior face-h is the
      !! first-order upwind pick — consistent with `compute_h_face_upstream`.
      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 — and required on purpose.  Under `&vcoord_nml
         !! zfixed_closed_faces` a CLOSED layer carries no transport, so
         !! it must not dilute the face mean either: the weight becomes
         !! `h_face·open` and `ubt` is the OPEN-column depth mean.  That
         !! is not a refinement, it is a consistency requirement: the
         !! barotropic substep transports on `ubt·FA·dy_cu_bt` with
         !! `dy_cu_bt` narrowed by the open fraction, so the `ubt` the
         !! fast loop integrates ALREADY means "the open-column mean".
         !! Deriving `ubt_at_n` from the full column would make the
         !! fold's `Δu = ubt_end − ubt_at_n − dt·F_bt` a difference
         !! between two different quantities — diluted by `0.5·h_live`
         !! per one-sided-filler layer, which at a partial face is not
         !! small.
         !!
         !! This dummy was OPTIONAL for one release and that is exactly
         !! how the `pred_corr` Coriolis-reference defect shipped: one
         !! call site in `set_cor_ref_velocity` omitted it, silently took
         !! the full-column branch, and the un-cancelled `f·(1−φ)·v̄`
         !! forced every barotropic substep.  An optional argument that
         !! silently changes the physics is a defect CLASS, not a defect;
         !! passing it is now mandatory so a new call site cannot quietly
         !! take the wrong branch.
         !!
         !! `metrics%use_closed_faces = .false.` (the default) ⇒ the
         !! ORIGINAL loops run, textually unchanged (byte-identical).  The
         !! open branch is written out in full rather than folded into the
         !! original with a runtime `if` so that the default-path kernel
         !! never NAMES `open_u`/`open_v` at all — with the knob off those
         !! are `(1,1,1)` placeholders, and a placeholder indexed inside a
         !! `do concurrent` is exactly the kind of thing that works on the
         !! host and faults under `mem:separate`.

      integer :: i, j, k, nx, ny, nz, nx_face, ny_face, il, ir, jl, jr
      logical :: use_upstream, use_open
      real(wp) :: total_h, hu_sum, h_face_sum, h_face, hv_sum
      real(wp) :: href_l, href_r, e_prev
      integer :: scheme

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      nx_face = size(ms%u_face_x_layer, 1)
      ny_face = size(ms%v_face_y_layer, 2)
      use_upstream = bt_work%use_upstream_h_face
      use_open = metrics%use_closed_faces
      ! Gated exactly like `compute_bt_rem_from_visc_rem`'s av_rem (see that
      ! routine's `av_rem_scheme` docstring, fix(bt) d885304234): HYBRID's
      ! shelf test is built for a genuine z_fixed/zstar partial-step bed, the
      ! only case `metrics%use_closed_faces` represents.  Off that family
      ! (hycom, eulerian_z, rho, sigma, lagrangian, ...) `href`/`e_prev`'s
      ! "shelf" has no such meaning, and HYBRID measurably corrupts the
      ! barotropic coupling there -- see FRHAT_HYBRID's gate in
      ! `face_depth_mean_u` for the full writeup and the two matrix
      ! witnesses (`decomp_hycom_visc_rem_chain`, `budget_eulerian_z_bc_pgf`)
      ! that found it.
      scheme = merge(bt_work%frhat_scheme, FRHAT_ARITHMETIC, use_open)

      do concurrent(j=1:ny, i=1:nx) local(k, total_h)
         total_h = 0.0_wp
         do k = 1, nz
            total_h = total_h + ms%h_layer(i, j, k)
         end do
         bt_work%bt_eta(i, j) = total_h - bt_work%bt_H_ref(i, j)
      end do

      ! `il`/`ir` are the two columns a u-face abuts (equal at a wall, where
      ! `frhat_h_face_step` degenerates to the single available column —
      ! see that routine's docstring).  `use_upstream` (no MOM6 frhat
      ! counterpart — a roundabout-only alternative orthogonal to
      ! centred/hybrid) bypasses the frhat sweep entirely and is untouched.
      if (use_open) then
         do concurrent(j=1:ny, i=1:nx_face) &
            local(k, il, ir, href_l, href_r, e_prev, hu_sum, h_face_sum, h_face)
            il = max(1, i - 1)
            ir = min(nx, i)
            href_l = bt_work%bt_H_ref(il, j)
            href_r = bt_work%bt_H_ref(ir, j)
            e_prev = -0.5_wp*(href_l + href_r)
            hu_sum = 0.0_wp
            h_face_sum = 0.0_wp
            do k = 1, nz
               if (use_upstream .and. i > 1 .and. i < nx_face) then
                  if (ms%u_face_x_layer(i, j, k) >= 0.0_wp) then
                     h_face = ms%h_layer(i - 1, j, k)
                  else
                     h_face = ms%h_layer(i, j, k)
                  end if
               else
                  call frhat_h_face_step(ms%h_layer(il, j, k), ms%h_layer(ir, j, k), &
                                         href_l, href_r, scheme, e_prev, h_face)
               end if
               h_face = h_face*metrics%open_u(i, j, k)
               hu_sum = hu_sum + ms%u_face_x_layer(i, j, k)*h_face
               h_face_sum = h_face_sum + h_face
            end do
            if (h_face_sum > 0.0_wp) then
               bt_work%bt_ubt(i, j) = hu_sum/h_face_sum
            else
               bt_work%bt_ubt(i, j) = 0.0_wp
            end if
         end do
      else
         do concurrent(j=1:ny, i=1:nx_face) &
            local(k, il, ir, href_l, href_r, e_prev, hu_sum, h_face_sum, h_face)
            il = max(1, i - 1)
            ir = min(nx, i)
            href_l = bt_work%bt_H_ref(il, j)
            href_r = bt_work%bt_H_ref(ir, j)
            e_prev = -0.5_wp*(href_l + href_r)
            hu_sum = 0.0_wp
            h_face_sum = 0.0_wp
            do k = 1, nz
               if (use_upstream .and. i > 1 .and. i < nx_face) then
                  if (ms%u_face_x_layer(i, j, k) >= 0.0_wp) then
                     h_face = ms%h_layer(i - 1, j, k)
                  else
                     h_face = ms%h_layer(i, j, k)
                  end if
               else
                  call frhat_h_face_step(ms%h_layer(il, j, k), ms%h_layer(ir, j, k), &
                                         href_l, href_r, scheme, e_prev, h_face)
               end if
               hu_sum = hu_sum + ms%u_face_x_layer(i, j, k)*h_face
               h_face_sum = h_face_sum + h_face
            end do
            if (h_face_sum > 0.0_wp) then
               bt_work%bt_ubt(i, j) = hu_sum/h_face_sum
            else
               bt_work%bt_ubt(i, j) = 0.0_wp
            end if
         end do
      end if

      if (use_open) then
         do concurrent(j=1:ny_face, i=1:nx) &
            local(k, jl, jr, href_l, href_r, e_prev, hv_sum, h_face_sum, h_face)
            jl = max(1, j - 1)
            jr = min(ny, j)
            href_l = bt_work%bt_H_ref(i, jl)
            href_r = bt_work%bt_H_ref(i, jr)
            e_prev = -0.5_wp*(href_l + href_r)
            hv_sum = 0.0_wp
            h_face_sum = 0.0_wp
            do k = 1, nz
               if (use_upstream .and. j > 1 .and. j < ny_face) then
                  if (ms%v_face_y_layer(i, j, k) >= 0.0_wp) then
                     h_face = ms%h_layer(i, j - 1, k)
                  else
                     h_face = ms%h_layer(i, j, k)
                  end if
               else
                  call frhat_h_face_step(ms%h_layer(i, jl, k), ms%h_layer(i, jr, k), &
                                         href_l, href_r, scheme, e_prev, h_face)
               end if
               h_face = h_face*metrics%open_v(i, j, k)
               hv_sum = hv_sum + ms%v_face_y_layer(i, j, k)*h_face
               h_face_sum = h_face_sum + h_face
            end do
            if (h_face_sum > 0.0_wp) then
               bt_work%bt_vbt(i, j) = hv_sum/h_face_sum
            else
               bt_work%bt_vbt(i, j) = 0.0_wp
            end if
         end do
      else
         do concurrent(j=1:ny_face, i=1:nx) &
            local(k, jl, jr, href_l, href_r, e_prev, hv_sum, h_face_sum, h_face)
            jl = max(1, j - 1)
            jr = min(ny, j)
            href_l = bt_work%bt_H_ref(i, jl)
            href_r = bt_work%bt_H_ref(i, jr)
            e_prev = -0.5_wp*(href_l + href_r)
            hv_sum = 0.0_wp
            h_face_sum = 0.0_wp
            do k = 1, nz
               if (use_upstream .and. j > 1 .and. j < ny_face) then
                  if (ms%v_face_y_layer(i, j, k) >= 0.0_wp) then
                     h_face = ms%h_layer(i, j - 1, k)
                  else
                     h_face = ms%h_layer(i, j, k)
                  end if
               else
                  call frhat_h_face_step(ms%h_layer(i, jl, k), ms%h_layer(i, jr, k), &
                                         href_l, href_r, scheme, e_prev, h_face)
               end if
               hv_sum = hv_sum + ms%v_face_y_layer(i, j, k)*h_face
               h_face_sum = h_face_sum + h_face
            end do
            if (h_face_sum > 0.0_wp) then
               bt_work%bt_vbt(i, j) = hv_sum/h_face_sum
            else
               bt_work%bt_vbt(i, j) = 0.0_wp
            end if
         end do
      end if
   end subroutine derive_bt_from_layers