print_bt_budget Subroutine

   subroutine print_bt_budget(grid, ms, bt_work, pgf, cor, hv, bd, ss, &
                              t_value, t_unit, stage_tag, header)
      !! Compute + print the per-region BT-mode budget snapshot.
      !!
      !! Power terms `P_<term>` are computed as the cell-centred
      !! u_bt·F_u + v_bt·F_v dot product, then averaged over the
      !! cells in each region.  Units: m²/s³ (acceleration × velocity).
      !! Sign convention: positive = the term ADDS energy to the BT
      !! mode (the wind P_ss should be positive in the gyres; hvisc /
      !! drag should be negative).
      type(hgrid_t), intent(in) :: grid
      type(multilayer_state_t), intent(in) :: ms
      type(barotropic_workstate_t), intent(in) :: bt_work
      type(ocean_pressure_force_t), intent(in) :: pgf
      type(coriolis_adv_t), intent(in) :: cor
      type(ocean_horizontal_viscosity_t), intent(in) :: hv
      type(ocean_bottom_drag_t), intent(in) :: bd
      type(ocean_surface_stress_t), intent(in) :: ss
      real(wp), intent(in) :: t_value
      character(len=*), intent(in) :: t_unit, stage_tag
      logical, intent(in) :: header

      integer :: nx, ny, nz, nghost
      integer :: jp_lo, jp_mid_lo, jp_mid_hi, jp_hi
      integer :: nz_top
      real(wp) :: ke_sub, ke_jet, ke_pol
      real(wp) :: eta_sub_min, eta_sub_max, eta_jet_min, eta_jet_max
      real(wp) :: eta_pol_min, eta_pol_max
      real(wp) :: usurf_sub, usurf_jet, usurf_pol
      real(wp) :: ubed_sub, ubed_jet, ubed_pol
      real(wp) :: vsurf_sub, vsurf_jet, vsurf_pol
      real(wp) :: vbed_sub, vbed_jet, vbed_pol
      real(wp) :: p_pgf_sub, p_pgf_jet, p_pgf_pol
      real(wp) :: p_cor_sub, p_cor_jet, p_cor_pol
      real(wp) :: p_hv_sub, p_hv_jet, p_hv_pol
      real(wp) :: p_bd_sub, p_bd_jet, p_bd_pol
      real(wp) :: p_ss_sub, p_ss_jet, p_ss_pol

      ! ---- Region selection: thirds along j (NS) within the
      ! physical domain.  The double_gyre setup puts the spoon's
      ! east-west margin in the top third (subpolar gyre).
      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      nghost = grid%nghost
      nz_top = nz                            ! Surface layer index
      jp_lo = nghost + 1                     ! First interior j
      jp_hi = nghost + grid%ny_phys          ! Last interior j
      jp_mid_lo = nghost + grid%ny_phys/3 + 1
      jp_mid_hi = nghost + 2*grid%ny_phys/3

      ! ---- η ranges + |u|/|v| extremes per region ----
      call region_eta_uv(grid, ms, bt_work, &
                         jp_lo, jp_mid_lo - 1, &
                         eta_sub_min, eta_sub_max, &
                         usurf_sub, ubed_sub, vsurf_sub, vbed_sub, ke_sub)
      call region_eta_uv(grid, ms, bt_work, &
                         jp_mid_lo, jp_mid_hi, &
                         eta_jet_min, eta_jet_max, &
                         usurf_jet, ubed_jet, vsurf_jet, vbed_jet, ke_jet)
      call region_eta_uv(grid, ms, bt_work, &
                         jp_mid_hi + 1, jp_hi, &
                         eta_pol_min, eta_pol_max, &
                         usurf_pol, ubed_pol, vsurf_pol, vbed_pol, ke_pol)

      ! ---- Per-term BT power per region ----
      call region_power(grid, ms, bt_work, &
                        pgf%dpdx_face%data, pgf%dpdy_face%data, &
                        jp_lo, jp_mid_lo - 1, p_pgf_sub)
      call region_power(grid, ms, bt_work, &
                        pgf%dpdx_face%data, pgf%dpdy_face%data, &
                        jp_mid_lo, jp_mid_hi, p_pgf_jet)
      call region_power(grid, ms, bt_work, &
                        pgf%dpdx_face%data, pgf%dpdy_face%data, &
                        jp_mid_hi + 1, jp_hi, p_pgf_pol)

      call region_power(grid, ms, bt_work, &
                        cor%pv_flux_x%data, cor%pv_flux_y%data, &
                        jp_lo, jp_mid_lo - 1, p_cor_sub)
      call region_power(grid, ms, bt_work, &
                        cor%pv_flux_x%data, cor%pv_flux_y%data, &
                        jp_mid_lo, jp_mid_hi, p_cor_jet)
      call region_power(grid, ms, bt_work, &
                        cor%pv_flux_x%data, cor%pv_flux_y%data, &
                        jp_mid_hi + 1, jp_hi, p_cor_pol)

      call region_power(grid, ms, bt_work, &
                        hv%du_visc%data, hv%dv_visc%data, &
                        jp_lo, jp_mid_lo - 1, p_hv_sub)
      call region_power(grid, ms, bt_work, &
                        hv%du_visc%data, hv%dv_visc%data, &
                        jp_mid_lo, jp_mid_hi, p_hv_jet)
      call region_power(grid, ms, bt_work, &
                        hv%du_visc%data, hv%dv_visc%data, &
                        jp_mid_hi + 1, jp_hi, p_hv_pol)

      ! Bottom drag: explicit-tendency path writes `-r·u·(h_in_bbl/h)`
      ! to `du_drag/dv_drag` directly. (Implicit-rate variant handled by
      ! region_power_drag_implicit when used.)
      call region_power(grid, ms, bt_work, &
                        bd%du_drag%data, bd%dv_drag%data, &
                        jp_lo, jp_mid_lo - 1, p_bd_sub)
      call region_power(grid, ms, bt_work, &
                        bd%du_drag%data, bd%dv_drag%data, &
                        jp_mid_lo, jp_mid_hi, p_bd_jet)
      call region_power(grid, ms, bt_work, &
                        bd%du_drag%data, bd%dv_drag%data, &
                        jp_mid_hi + 1, jp_hi, p_bd_pol)

      call region_power(grid, ms, bt_work, &
                        ss%du_stress%data, ss%dv_stress%data, &
                        jp_lo, jp_mid_lo - 1, p_ss_sub)
      call region_power(grid, ms, bt_work, &
                        ss%du_stress%data, ss%dv_stress%data, &
                        jp_mid_lo, jp_mid_hi, p_ss_jet)
      call region_power(grid, ms, bt_work, &
                        ss%du_stress%data, ss%dv_stress%data, &
                        jp_mid_hi + 1, jp_hi, p_ss_pol)

      ! ---- Print ----
      if (header) then
         write (*, "(a)") "# BT-BUDGET PROBE"
         write (*, "(a)") "# columns: t stage region <KE>[m2/s2] eta_min eta_max "// &
            "|u_surf|max |u_bed|max |v_surf|max |v_bed|max "// &
            "P_pgf P_cor P_hv P_bd P_ss   [m2/s3]"
         write (*, "(a)") "# regions: SUB=subtropical (low j), JET=middle, POL=subpolar (high j)"
      end if

      call emit_row(t_value, t_unit, stage_tag, "SUB", ke_sub, &
                    eta_sub_min, eta_sub_max, &
                    usurf_sub, ubed_sub, vsurf_sub, vbed_sub, &
                    p_pgf_sub, p_cor_sub, p_hv_sub, p_bd_sub, p_ss_sub)
      call emit_row(t_value, t_unit, stage_tag, "JET", ke_jet, &
                    eta_jet_min, eta_jet_max, &
                    usurf_jet, ubed_jet, vsurf_jet, vbed_jet, &
                    p_pgf_jet, p_cor_jet, p_hv_jet, p_bd_jet, p_ss_jet)
      call emit_row(t_value, t_unit, stage_tag, "POL", ke_pol, &
                    eta_pol_min, eta_pol_max, &
                    usurf_pol, ubed_pol, vsurf_pol, vbed_pol, &
                    p_pgf_pol, p_cor_pol, p_hv_pol, p_bd_pol, p_ss_pol)
   end subroutine print_bt_budget