ocean_dyn_step_barotropic Subroutine

   pure subroutine ocean_dyn_step_barotropic(grid, metrics, dyn, cor, ct, bs, dt)
      !! Public only for the unit-test suite (no production module imports it);
      !! ignore when developing production code in other modules.
      !! Unsplit SSP-RK2 (Heun's method) outer step on the barotropic
      !! C-grid state.  Couples continuity-PPM (h-update) and the
      !! Sadourny Coriolis-advection tendency (u, v update) into one
      !! second-order-accurate step.
      !!
      !! SSP-RK2 (Shu-Osher form):
      !!   u^(1)    = u^n + dt * L(u^n)
      !!   u^(n+1)  = 1/2 * u^n + 1/2 * (u^(1) + dt * L(u^(1)))
      !!
      !! In our code:
      !!   1. Save u^n into the *_0 buffers on `bs`.
      !!   2. Stage 1: compute continuity flux divergence and the
      !!      Sadourny momentum tendency at u^n, apply both as a
      !!      forward-Euler step.  State is now u^(1).
      !!   3. Stage 2: compute tendencies at u^(1), apply another
      !!      forward-Euler step.  State is u^(1) + dt * L(u^(1)).
      !!   4. RK2 average: state <- 1/2 * (u_0 + state).
      !!
      !! Errors per step are O(dt^3); compared to plain FE, the
      !! inertial-oscillator amplitude growth drops from O(dt^2) per
      !! step to O(dt^4) per step.  At f*dt = 0.01 over 157 steps the
      !! magnitude error drops from ~0.8% (FE) to ~2e-7 (RK2).
      !!
      !! `dyn` carries diagnostic state (step counter, CFL, KE) for
      !! the outer step; Phase 4a only bumps `outer_step_count`.
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
      type(ocean_dyn_t), intent(inout) :: dyn
      type(coriolis_adv_t), intent(inout) :: cor
      type(continuity_t), intent(inout) :: ct
      type(barotropic_state_t), intent(inout) :: bs
      real(wp), intent(in) :: dt

      integer :: i, j, nx, ny, nx_face, ny_uface, nx_vface, ny_face

      nx = grid%nx_total
      ny = grid%ny_total
      nx_face = size(bs%u_face_x, 1)
      ny_uface = size(bs%u_face_x, 2)
      nx_vface = size(bs%v_face_y, 1)
      ny_face = size(bs%v_face_y, 2)

      ! ---- 1. Save u^n into h0 / u_face_x0 / v_face_y0 ----
      do concurrent(j=1:ny, i=1:nx)
         bs%h0(i, j) = bs%h(i, j)
      end do
      do concurrent(j=1:ny_uface, i=1:nx_face)
         bs%u_face_x0(i, j) = bs%u_face_x(i, j)
      end do
      do concurrent(j=1:ny_face, i=1:nx_vface)
         bs%v_face_y0(i, j) = bs%v_face_y(i, j)
      end do

      ! ---- 2. Stage 1: tendencies at u^n, FE step -> u^(1) ----
      call continuity_compute_fluxes_barotropic(grid, metrics, ct, bs)
      call coriolis_adv_compute_tendencies_barotropic(grid, metrics, cor, bs)
      call continuity_apply_fluxes_barotropic(bs, dt)
      call coriolis_adv_apply_tendencies_barotropic(cor, bs, dt)

      ! ---- 3. Stage 2: tendencies at u^(1), FE step -> u^(1) + dt*L(u^(1)) ----
      call continuity_compute_fluxes_barotropic(grid, metrics, ct, bs)
      call coriolis_adv_compute_tendencies_barotropic(grid, metrics, cor, bs)
      call continuity_apply_fluxes_barotropic(bs, dt)
      call coriolis_adv_apply_tendencies_barotropic(cor, bs, dt)

      ! ---- 4. RK2 average: u^(n+1) = 1/2 * (u^n + (u^(1) + dt*L(u^(1)))) ----
      do concurrent(j=1:ny, i=1:nx)
         bs%h(i, j) = 0.5_wp*(bs%h0(i, j) + bs%h(i, j))
      end do
      do concurrent(j=1:ny_uface, i=1:nx_face)
         bs%u_face_x(i, j) = 0.5_wp*(bs%u_face_x0(i, j) + bs%u_face_x(i, j))
      end do
      do concurrent(j=1:ny_face, i=1:nx_vface)
         bs%v_face_y(i, j) = 0.5_wp*(bs%v_face_y0(i, j) + bs%v_face_y(i, j))
      end do

      dyn%outer_step_count = dyn%outer_step_count + 1
   end subroutine ocean_dyn_step_barotropic