subtract_fast_cor_ref Subroutine

   pure subroutine subtract_fast_cor_ref(grid, metrics, bt_work, f_corner, &
                                         bc_w, bc_e, bc_s, bc_n, &
                                         has_w, has_e, has_s, has_n)
      !! Subtract the fast-loop Coriolis + vector-invariant advection,
      !! evaluated at the reference barotropic velocity
      !! `bt_work%cor_ref_u`/`cor_ref_v` (filled by
      !! `set_cor_ref_velocity`), from the substep forcing
      !! `F_bt_u_fast`/`F_bt_v_fast`.
      !!
      !! The reference velocity is NOT free: it must be the depth mean
      !! of the same layer velocity the slow `cor%pv_flux_*` inside
      !! `F_bt_u/v` was evaluated on, or the difference survives as a
      !! near-constant per-substep forcing.  Under `ssp_rk2` that is the
      !! stage-entry `bt_ubt/bt_vbt` (what this routine used to read
      !! directly); under `pred_corr` it is the depth mean of `u_av/v_av`.
      !! See `set_cor_ref_velocity` and the `cor_ref_u` docstring in
      !! `rdb_barotropic_workstate`.
      !!
      !! Why: `F_bt_u` is the depth mean of ALL slow layer tendencies —
      !! including the layer Coriolis-advection (`cor%pv_flux_*`), whose
      !! depth mean is ≈ the fast solver's own `(ζ+f)·v − ∇KE` at the
      !! stage-entry state.  The substep then integrates its own LIVE
      !! `(ζ+f)·v − ∇KE` on top, so without this subtraction the
      !! barotropic Coriolis/advection is integrated TWICE — the exact
      !! analogue of the PGF double-count `set_fast_forcing_eta_pf`
      !! already guards against by shedding the free-surface term the
      !! slow PGF carries ("√(gH) inflates to √(2gH)").  The Coriolis double-count is what pumps
      !! the exponential wall/corner barotropic mode on shelf rims under
      !! `VCOORD_LAGRANGIAN` (600² double-gyre h-guard trap).  MOM6
      !! removes it with a reference Coriolis/advection (`Cor_ref_u/v`)
      !! subtracted inside its barotropic substep loop; this is that
      !! subtraction, folded into the forcing so the substep kernel is
      !! untouched.
      !!
      !! After this, at τ=0 the substep's net Coriolis/advection
      !! contribution is zero and only the ANOMALY that develops over
      !! the substeps is integrated — so `Δu = ubt_end − ubt_at_n −
      !! dt·F_bt_u` hands the layers the genuine fast anomaly instead
      !! of an extra `dt·f·v̄` rotation per stage.
      !!
      !! Uses `bt_work%bt_zeta_corner` and `bt_work%bt_ke_centre` as
      !! scratch — both are per-substep scratch the substep recomputes
      !! from scratch before reading (Pass 1 precedes Pass 2b/2c).
      !! ζ/KE formulas and wall closures mirror the substep exactly
      !! (`rdb_barotropic_substep` Pass 1 + the corner-ζ closure), so
      !! the τ=0 cancellation holds at walls too — where the unstable
      !! mode lives.  bt_halo = 0 index convention (the wide-halo
      !! march-in path receives the same interior-computed forcing via
      !! `copy_in`; its seam rows differ at O(ghost) — acceptable, the
      !! halo exchange owns them).
      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) :: f_corner(grid%nx_total + 1, grid%ny_total + 1)
      integer, intent(in) :: bc_w, bc_e, bc_s, bc_n
         !! Per-edge OBC tags (OBC_WALL when no bc present).
      logical, intent(in) :: has_w, has_e, has_s, has_n
         !! Physical-edge flags (.false. at an MPI seam).

      integer :: i, j, nx, ny
      real(wp) :: zeta_at_u, f_at_u, v_at_u, ke_grad_x
      real(wp) :: zeta_at_v, f_at_v, u_at_v, ke_grad_y
      real(wp) :: w_nl
         !! Mirror of the substep's live-nonlinear weight
         !! (`bt_work%substep_zeta_ke`): the reference MUST subtract
         !! exactly what the substep re-integrates live — full
         !! `(ζ+f)·v − ∇KE` when live (1), planetary `f·v̄` only when
         !! the substep is planetary-only (0).

      nx = grid%nx_total
      ny = grid%ny_total
      w_nl = merge(1.0_wp, 0.0_wp, bt_work%substep_zeta_ke)

      ! ---- ζ at corners from the stage-entry bt state (substep Pass-1
      ! formula + closure).  Corner ring + physical-wall lines → 0.
      do concurrent(j=1:ny + 1, i=1:nx + 1)
         if (i >= 2 .and. i <= nx .and. j >= 2 .and. j <= ny) then
            bt_work%bt_zeta_corner(i, j) = &
               ((bt_work%cor_ref_v(i, j)*metrics%dyCv(i, j) &
                 - bt_work%cor_ref_v(i - 1, j)*metrics%dyCv(i - 1, j)) - &
                (bt_work%cor_ref_u(i, j)*metrics%dxCu(i, j) &
                 - bt_work%cor_ref_u(i, j - 1)*metrics%dxCu(i, j - 1)))* &
               metrics%iareaBu(i, j)
         else
            bt_work%bt_zeta_corner(i, j) = 0.0_wp
         end if
         if (bc_w /= OBC_PERIODIC .and. has_w .and. i == grid%nghost + 1) then
            bt_work%bt_zeta_corner(i, j) = 0.0_wp
         end if
         if (bc_e /= OBC_PERIODIC .and. has_e .and. i == grid%nghost + grid%nx_phys + 1) then
            bt_work%bt_zeta_corner(i, j) = 0.0_wp
         end if
         if (bc_s /= OBC_PERIODIC .and. has_s .and. j == grid%nghost + 1) then
            bt_work%bt_zeta_corner(i, j) = 0.0_wp
         end if
         ! The tripolar fold line is a seam (interior corners), not a wall.
         if (bc_n /= OBC_PERIODIC .and. bc_n /= OBC_TRIPOLAR_FOLD .and. has_n .and. &
             j == grid%nghost + grid%ny_phys + 1) then
            bt_work%bt_zeta_corner(i, j) = 0.0_wp
         end if
      end do

      ! ---- KE at centres (substep Pass-1 formula). ----
      do concurrent(j=1:ny, i=1:nx)
         bt_work%bt_ke_centre(i, j) = 0.25_wp*metrics%iareaT(i, j)*( &
                                      metrics%areaCu(i, j)*bt_work%cor_ref_u(i, j)**2 + &
                                      metrics%areaCu(i + 1, j)*bt_work%cor_ref_u(i + 1, j)**2 + &
                                      metrics%areaCv(i, j)*bt_work%cor_ref_v(i, j)**2 + &
                                      metrics%areaCv(i, j + 1)*bt_work%cor_ref_v(i, j + 1)**2)
      end do

      ! ---- Subtract the u-face reference (substep Pass-2b operand). ----
      do concurrent(j=1:ny, i=2:nx) local(zeta_at_u, f_at_u, v_at_u, ke_grad_x)
         zeta_at_u = w_nl*0.5_wp*(bt_work%bt_zeta_corner(i, j) + bt_work%bt_zeta_corner(i, j + 1))
         f_at_u = 0.5_wp*(f_corner(i, j) + f_corner(i, j + 1))
         if (j > 1 .and. j < ny) then
            v_at_u = 0.25_wp*(bt_work%cor_ref_v(i - 1, j) + bt_work%cor_ref_v(i - 1, j + 1) + &
                              bt_work%cor_ref_v(i, j) + bt_work%cor_ref_v(i, j + 1))
         else if (j == 1) then
            v_at_u = 0.5_wp*(bt_work%cor_ref_v(i - 1, j + 1) + bt_work%cor_ref_v(i, j + 1))
         else
            v_at_u = 0.5_wp*(bt_work%cor_ref_v(i - 1, j) + bt_work%cor_ref_v(i, j))
         end if
         ke_grad_x = w_nl*(bt_work%bt_ke_centre(i, j) - bt_work%bt_ke_centre(i - 1, j))*metrics%idxCu(i, j)
         bt_work%F_bt_u_fast(i, j) = bt_work%F_bt_u_fast(i, j) - &
                                     ((zeta_at_u + f_at_u)*v_at_u - ke_grad_x)
      end do

      ! ---- Subtract the v-face reference (substep Pass-2c operand). ----
      do concurrent(j=2:ny, i=1:nx) local(zeta_at_v, f_at_v, u_at_v, ke_grad_y)
         zeta_at_v = w_nl*0.5_wp*(bt_work%bt_zeta_corner(i, j) + bt_work%bt_zeta_corner(i + 1, j))
         f_at_v = 0.5_wp*(f_corner(i, j) + f_corner(i + 1, j))
         if (i > 1 .and. i < nx) then
            u_at_v = 0.25_wp*(bt_work%cor_ref_u(i, j - 1) + bt_work%cor_ref_u(i + 1, j - 1) + &
                              bt_work%cor_ref_u(i, j) + bt_work%cor_ref_u(i + 1, j))
         else if (i == 1) then
            u_at_v = 0.5_wp*(bt_work%cor_ref_u(i + 1, j - 1) + bt_work%cor_ref_u(i + 1, j))
         else
            u_at_v = 0.5_wp*(bt_work%cor_ref_u(i, j - 1) + bt_work%cor_ref_u(i, j))
         end if
         ke_grad_y = w_nl*(bt_work%bt_ke_centre(i, j) - bt_work%bt_ke_centre(i, j - 1))*metrics%idyCv(i, j)
         bt_work%F_bt_v_fast(i, j) = bt_work%F_bt_v_fast(i, j) - &
                                     (-(zeta_at_v + f_at_v)*u_at_v - ke_grad_y)
      end do
   end subroutine subtract_fast_cor_ref