coriolis_adv_compute_tendencies_barotropic Subroutine

   pure subroutine coriolis_adv_compute_tendencies_barotropic(grid, metrics, this, bs)
      !! Sadourny (1975) energy-conserving Coriolis + horizontal-
      !! momentum-advection form on the barotropic C-grid state:
      !!
      !!   du/dt = +(zeta + f) * v_at_u_face - d/dx(KE)
      !!   dv/dt = -(zeta + f) * u_at_v_face - d/dy(KE)
      !!
      !! where zeta = dv/dx - du/dy is the relative vorticity at
      !! cell corners (south-west corner of cell (i, j) sits at
      !! position (i-1/2, j-1/2)) and KE = (1/2)*(u^2 + v^2) is the
      !! kinetic energy per unit mass, averaged from the surrounding
      !! face velocities at each cell centre.
      !!
      !! For spatially uniform u, v this reduces to plain Coriolis:
      !! zeta vanishes and d/dx(KE) vanishes, leaving the f*v / -f*u
      !! pair (the Phase 3a kernel).
      !!
      !! Three passes:
      !!   1. zeta at corners (q_corner buffer; reused as "vorticity"
      !!      until Phase 5 adds the q = (zeta+f)/h division for
      !!      layered PV).
      !!   2. KE at cell centres (ke_centre buffer), using
      !!      KE = (1/4)*(u_W^2 + u_E^2 + v_S^2 + v_N^2) for the
      !!      energy-consistent form.
      !!   3. Tendencies at faces, combining the vorticity-flux and
      !!      KE-gradient terms.
      !!
      !! Walls: zeta uses zero fallback at outer corners (where the
      !! 4-point stencil falls off the grid); face tendencies at
      !! domain walls are zeroed (closed-wall, no momentum at the
      !! boundary).
      !!
      !! Loop order: j-then-i for NVHPC GPU coalescing.
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
      type(coriolis_adv_t), intent(inout) :: this
      type(barotropic_state_t), intent(in) :: bs

      integer :: i, j, nx, ny
      real(wp) :: v_at_u, u_at_v, zeta_at_u, zeta_at_v
      real(wp) :: f_at_u, f_at_v, ke_grad_x, ke_grad_y
      real(wp) :: h_vf_SW, h_vf_NW, h_vf_SE, h_vf_NE
      real(wp) :: h_uf_SW, h_uf_NW, h_uf_SE, h_uf_NE
      real(wp) :: vh_sum, uh_sum, h_eff_sum
      real(wp) :: ns

      nx = grid%nx_total
      ny = grid%ny_total
      ! Slip selector (C1): ns=0 ⇒ factor=wet_q (free-slip); ns=1 ⇒
      ! factor=2-wet_q (no-slip image vorticity).  Branchless per corner.
      ns = merge(1.0_wp, 0.0_wp, this%no_slip)

      ! ---- Pass 1: relative vorticity at SW corners (circulation/area) ----
      ! Inline twin of `rdb_rvc_zeta_corner` (shared_module_utilities/
      ! rdb_rel_vort_corner.inc, read by the `vorticity_z` diag): keep in step.
      ! zeta_corner(i, j) sits at position (i-1/2, j-1/2).  Curvilinear
      ! circulation form (design §2):
      !   zeta = ( v(i,j)·dyCv(i,j) - v(i-1,j)·dyCv(i-1,j)
      !          - (u(i,j)·dxCu(i,j) - u(i,j-1)·dxCu(i,j-1)) ) · iareaBu(i,j)
      ! On uniform square metrics this is `(Δv)/dx - (Δu)/dy` bitwise.
      ! Interior corners: i=2..nx, j=2..ny.  Outer corners
      ! (i=1, j=1, i=nx+1, j=ny+1) get zero — no neighbouring cell
      ! across the wall.  The rel-vort is multiplied by the slip factor
      ! `(1-2·ns)·wet_q + 2·ns` (C1, MOM6 free-slip/no-slip); planetary f added
      ! later stays UNMASKED.  Interior land corners zero exactly as the
      ! domain-wall corners do.
      do concurrent(j=2:ny, i=2:nx)
         this%q_corner%data(i, j, 1) = &
            ((1.0_wp - 2.0_wp*ns)*metrics%wet_q(i, j) + 2.0_wp*ns)* &
            ((bs%v_face_y(i, j)*metrics%dyCv(i, j) - bs%v_face_y(i - 1, j)*metrics%dyCv(i - 1, j)) - &
             (bs%u_face_x(i, j)*metrics%dxCu(i, j) - bs%u_face_x(i, j - 1)*metrics%dxCu(i, j - 1)))* &
            metrics%iareaBu(i, j)
      end do
      ! Outer-corner fallback: zero vorticity at the boundary.
      do concurrent(j=1:ny + 1)
         this%q_corner%data(1, j, 1) = 0.0_wp
         this%q_corner%data(nx + 1, j, 1) = 0.0_wp
      end do
      do concurrent(i=1:nx + 1)
         this%q_corner%data(i, 1, 1) = 0.0_wp
         this%q_corner%data(i, ny + 1, 1) = 0.0_wp
      end do

      ! ---- Pass 2: kinetic energy at cell centres (area-weighted) ----
      ! KE_centre(i,j) = 0.25·iareaT·( areaCu(i)·u(i)² + areaCu(i+1)·u(i+1)²
      !                              + areaCv(j)·v(j)² + areaCv(j+1)·v(j+1)² )
      ! On uniform metrics areaCu=areaCv=areaT, iareaT=1/areaT, so this
      ! reduces to the simple 0.25·(u²+u²+v²+v²) form bitwise.
      do concurrent(j=1:ny, i=1:nx)
         this%ke_centre%data(i, j, 1) = 0.25_wp*metrics%iareaT(i, j)*( &
                                        metrics%areaCu(i, j)*bs%u_face_x(i, j)**2 + &
                                        metrics%areaCu(i + 1, j)*bs%u_face_x(i + 1, j)**2 + &
                                        metrics%areaCv(i, j)*bs%v_face_y(i, j)**2 + &
                                        metrics%areaCv(i, j + 1)*bs%v_face_y(i, j + 1)**2)
      end do

      ! ---- Pass 3a: du/dt at interior east faces ----
      ! Thickness-weighted v at the u-face: average `bs%mass_flux_y`
      ! (= v · h_face) over the 4 abutting v-faces, divide by the
      ! sum of the v-face thicknesses.  See the multilayer kernel
      ! comment for the rationale — same form, dropped k axis.
      ! For uniform h this is bit-identical to the simple 4-point
      ! velocity average.
      do concurrent(j=1:ny, i=2:nx) &
         local(f_at_u, h_vf_SW, h_vf_NW, h_vf_SE, h_vf_NE, &
               vh_sum, h_eff_sum)
         h_vf_SW = 0.5_wp*(bs%h(i - 1, max(1, j - 1)) + bs%h(i - 1, j))
         h_vf_NW = 0.5_wp*(bs%h(i - 1, j) + bs%h(i - 1, min(ny, j + 1)))
         h_vf_SE = 0.5_wp*(bs%h(i, max(1, j - 1)) + bs%h(i, j))
         h_vf_NE = 0.5_wp*(bs%h(i, j) + bs%h(i, min(ny, j + 1)))
         vh_sum = (bs%v_face_y(i - 1, j)*h_vf_SW + bs%v_face_y(i - 1, j + 1)*h_vf_NW) + &
                  (bs%v_face_y(i, j)*h_vf_SE + bs%v_face_y(i, j + 1)*h_vf_NE)
         h_eff_sum = (h_vf_SW + h_vf_NW) + (h_vf_SE + h_vf_NE)
         if (h_eff_sum > 0.0_wp) then
            v_at_u = vh_sum/h_eff_sum
         else
            v_at_u = 0.0_wp
         end if
         zeta_at_u = 0.5_wp*(this%q_corner%data(i, j, 1) + &
                             this%q_corner%data(i, j + 1, 1))
         f_at_u = 0.5_wp*(this%f_corner(i, j) + this%f_corner(i, j + 1))
         ke_grad_x = (this%ke_centre%data(i, j, 1) - &
                      this%ke_centre%data(i - 1, j, 1))*metrics%idxCu(i, j)
         this%pv_flux_x%data(i, j, 1) = (zeta_at_u + f_at_u)*v_at_u - ke_grad_x
      end do
      do concurrent(j=1:ny)
         this%pv_flux_x%data(1, j, 1) = 0.0_wp
         this%pv_flux_x%data(nx + 1, j, 1) = 0.0_wp
      end do

      ! ---- Pass 3b: dv/dt at interior north faces ----
      do concurrent(j=2:ny, i=1:nx) &
         local(f_at_v, h_uf_SW, h_uf_NW, h_uf_SE, h_uf_NE, &
               uh_sum, h_eff_sum)
         h_uf_SW = 0.5_wp*(bs%h(max(1, i - 1), j - 1) + bs%h(i, j - 1))
         h_uf_SE = 0.5_wp*(bs%h(i, j - 1) + bs%h(min(nx, i + 1), j - 1))
         h_uf_NW = 0.5_wp*(bs%h(max(1, i - 1), j) + bs%h(i, j))
         h_uf_NE = 0.5_wp*(bs%h(i, j) + bs%h(min(nx, i + 1), j))
         uh_sum = (bs%u_face_x(i, j - 1)*h_uf_SW + bs%u_face_x(i + 1, j - 1)*h_uf_SE) + &
                  (bs%u_face_x(i, j)*h_uf_NW + bs%u_face_x(i + 1, j)*h_uf_NE)
         h_eff_sum = (h_uf_SW + h_uf_SE) + (h_uf_NW + h_uf_NE)
         if (h_eff_sum > 0.0_wp) then
            u_at_v = uh_sum/h_eff_sum
         else
            u_at_v = 0.0_wp
         end if
         zeta_at_v = 0.5_wp*(this%q_corner%data(i, j, 1) + &
                             this%q_corner%data(i + 1, j, 1))
         f_at_v = 0.5_wp*(this%f_corner(i, j) + this%f_corner(i + 1, j))
         ke_grad_y = (this%ke_centre%data(i, j, 1) - &
                      this%ke_centre%data(i, j - 1, 1))*metrics%idyCv(i, j)
         this%pv_flux_y%data(i, j, 1) = -(zeta_at_v + f_at_v)*u_at_v - ke_grad_y
      end do
      do concurrent(i=1:nx)
         this%pv_flux_y%data(i, 1, 1) = 0.0_wp
         this%pv_flux_y%data(i, ny + 1, 1) = 0.0_wp
      end do
   end subroutine coriolis_adv_compute_tendencies_barotropic