set_fast_forcing_eta_pf Subroutine

   pure subroutine set_fast_forcing_eta_pf(grid, metrics, bt_work, nx, ny, g_pf, eta_seam, use_seam)
      !! The barotropic substep's frozen forcing under the MOM6 split
      !! (`&ocean_bt_nml bc_pgf_forcing`, default on):
      !!
      !!     F_bt_u_fast = F_bt_u + g_pf·(η_PF(i) − η_PF(i−1))·idxCu
      !!
      !! `F_bt` holds the depth mean of the FULL slow layer PGF.  The
      !! substep integrates `−g_bt·∇η` live, so the one thing the forcing
      !! must shed is the free-surface term the slow PGF itself carries,
      !! `−g_pf·∇η_PF` (`pgf_free_surface_gravity`: 0 for the
      !! surface-relative MONT/FV_LITE/FV_WRIGHT forms, ≈ `g_bt` for
      !! FV_MOM6/GPRIME), evaluated at `η_PF = bt_eta` at stage entry —
      !! the free surface the slow PGF of this stage was built on
      !! (`derive_bt_from_layers` fills it from the same `h_layer`).  The
      !! substep then sees `⟨PGF_bc⟩ − g_bt·∇η`: the depth-mean BAROCLINIC
      !! pressure gradient plus its own free surface.  MOM6: `BT_force`
      !! (Σ wt·bc_accel, PFu included) with `btloop_find_PF`'s
      !! `−gtot·∇(η − eta_PF)`.
      !!
      !! The legacy split (`F_bt − ⟨PGF⟩`) shed the WHOLE depth-mean PGF —
      !! its baroclinic part too, which is the bottom-pressure gradient of
      !! a sloping density field (JEBAR).  Because `apply_bt_correction`
      !! subtracts `dt·F_bt` from every layer, the layers lost it as well:
      !! the depth mean ended every stage at `u_bt^end`, which never felt
      !! it.
      !!
      !! `eta_seam`/`use_seam`: when the surface-pressure load ALSO enters
      !! the slow PGF (`&ocean_pgf_nml p_top_in_bc` with the psurf seam on),
      !! `⟨PGF⟩` carries `−∇p_surf/ρ₀ = +g·∇η_ib` and the substep carries
      !! it again through `eta_forcing`; shedding `g_pf·∇η_ib` here counts
      !! it once.  Otherwise pass any full-size array and `.false.`.
      !!
      !! Array-edge faces (`i = 1`, `nx + 1`; `j = 1`, `ny + 1`) have no
      !! η on one side; the substep never reads their forcing (they are
      !! overwritten by the boundary dispatch / halo), so they take `F_bt`.
      integer, intent(in) :: nx, ny
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
      type(barotropic_workstate_t), intent(inout) :: bt_work
      real(wp), intent(in) :: g_pf
      real(wp), intent(in) :: eta_seam(nx, ny)
      logical, intent(in) :: use_seam
      integer :: i, j, nu, nv
      real(wp) :: d_eta

      if (.false.) nu = grid%nx_total
      nu = nx + 1
      nv = ny + 1
      do concurrent(j=1:ny, i=1:nu) local(d_eta)
         if (i == 1 .or. i == nu) then
            bt_work%F_bt_u_fast(i, j) = bt_work%F_bt_u(i, j)
         else
            d_eta = bt_work%bt_eta(i, j) - bt_work%bt_eta(i - 1, j)
            if (use_seam) d_eta = d_eta - (eta_seam(i, j) - eta_seam(i - 1, j))
            bt_work%F_bt_u_fast(i, j) = bt_work%F_bt_u(i, j) + g_pf*d_eta*metrics%idxCu(i, j)
         end if
      end do
      do concurrent(j=1:nv, i=1:nx) local(d_eta)
         if (j == 1 .or. j == nv) then
            bt_work%F_bt_v_fast(i, j) = bt_work%F_bt_v(i, j)
         else
            d_eta = bt_work%bt_eta(i, j) - bt_work%bt_eta(i, j - 1)
            if (use_seam) d_eta = d_eta - (eta_seam(i, j) - eta_seam(i, j - 1))
            bt_work%F_bt_v_fast(i, j) = bt_work%F_bt_v(i, j) + g_pf*d_eta*metrics%idyCv(i, j)
         end if
      end do
   end subroutine set_fast_forcing_eta_pf