continuity_compute_fluxes Subroutine

   pure subroutine continuity_compute_fluxes(grid, metrics, this, ms)
      !! **Test-only** (no production caller): the unsplit reference path,
      !! kept as the oracle the split production path is checked against.
      !! Multilayer counterpart to
      !! `continuity_compute_fluxes_barotropic`: identical PPM
      !! reconstruction + upwind face pick + flux divergence, lifted
      !! per-layer.  Each k-slice is independent (the PPM stencil
      !! reads only the same k), so the do-concurrent kernels
      !! parallelize over (k, j, i) simultaneously for GPU
      !! occupancy.
      !!
      !! Workspaces (`this%h_face_*_x/y`) must have been initialised
      !! with `nz_ml` matching `ms%nz_ml` — handled by passing the
      !! optional `nz_ml` to `continuity_init`.
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
      type(continuity_t), intent(inout) :: this
      type(multilayer_state_t), intent(inout) :: ms

      integer :: i, j, k, nx, ny, nz
      real(wp) :: dh_m1, dh_0, dh_p1, h_left, h_right, u, v, h_face
      real(wp) :: hm2, hm1, h0, hp1, hp2
      logical :: do_pos
      real(wp) :: h_min_pos

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      do_pos = this%use_ppm_limit_pos
      h_min_pos = this%h_min

      ! ============================================================
      ! X-DIRECTION reconstruction (5-point stencil per cell)
      ! ============================================================
      do concurrent(k=1:nz, j=1:ny, i=3:nx - 2) &
         local(dh_m1, dh_0, dh_p1, h_left, h_right, hm2, hm1, h0, hp1, hp2)
         ! Mirror-h at land neighbours (C2); bit-identical for all-wet.
         h0 = ms%h_layer(i, j, k)
         hm1 = ppm_mirror_h(ms%h_layer(i - 1, j, k), h0, metrics%wet_T(i - 1, j))
         hp1 = ppm_mirror_h(ms%h_layer(i + 1, j, k), h0, metrics%wet_T(i + 1, j))
         hm2 = ppm_mirror_h(ms%h_layer(i - 2, j, k), hm1, metrics%wet_T(i - 2, j))
         hp2 = ppm_mirror_h(ms%h_layer(i + 2, j, k), hp1, metrics%wet_T(i + 2, j))
         call ppm_limited_slope(hm2, hm1, h0, dh_m1)
         call ppm_limited_slope(hm1, h0, hp1, dh_0)
         call ppm_limited_slope(h0, hp1, hp2, dh_p1)
         dh_0 = dh_0*metrics%wet_T(i - 1, j)*metrics%wet_T(i, j)*metrics%wet_T(i + 1, j)
         h_left = 0.5_wp*(hm1 + h0) - (dh_0 - dh_m1)/6.0_wp
         h_right = 0.5_wp*(h0 + hp1) - (dh_p1 - dh_0)/6.0_wp
         call ppm_cell_limiter(h0, h_left, h_right)
         if (do_pos) call ppm_limit_pos(h0, h_left, h_right, h_min_pos)
         this%h_face_right_x%data(i, j, k) = h_left
         this%h_face_left_x%data(i + 1, j, k) = h_right
      end do
      ! Boundary cells: 1st-order fallback
      do concurrent(k=1:nz, j=1:ny)
         this%h_face_left_x%data(1, j, k) = ms%h_layer(1, j, k)
         this%h_face_right_x%data(1, j, k) = ms%h_layer(1, j, k)
         this%h_face_left_x%data(2, j, k) = ms%h_layer(1, j, k)
         this%h_face_right_x%data(2, j, k) = ms%h_layer(2, j, k)
         ! Face 3: cell 2's right edge falls back to 1st order (its
         ! 5-point stencil needs cell 0); cell 3's left edge is set
         ! by the interior loop above.
         this%h_face_left_x%data(3, j, k) = ms%h_layer(2, j, k)
         this%h_face_right_x%data(3, j, k) = ms%h_layer(2, j, k)
         ! Face nx-1: mirror of face 3.  Cell nx-1's left edge falls
         ! back to 1st order; cell nx-2's right edge came from the
         ! interior loop.
         this%h_face_right_x%data(nx - 1, j, k) = ms%h_layer(nx - 1, j, k)
         this%h_face_left_x%data(nx, j, k) = ms%h_layer(nx - 1, j, k)
         this%h_face_right_x%data(nx, j, k) = ms%h_layer(nx, j, k)
         this%h_face_left_x%data(nx + 1, j, k) = ms%h_layer(nx, j, k)
         this%h_face_right_x%data(nx + 1, j, k) = ms%h_layer(nx, j, k)
      end do

      ! X-direction face transport (upwind pick) — width-weighted dy_cu
      do concurrent(k=1:nz, j=1:ny, i=2:nx) local(u, h_face)
         u = ms%u_face_x_layer(i, j, k)
         if (u >= 0.0_wp) then
            h_face = this%h_face_left_x%data(i, j, k)
         else
            h_face = this%h_face_right_x%data(i, j, k)
         end if
         ms%mass_flux_x_layer(i, j, k) = u*h_face*metrics%dy_cu(i, j)
      end do
      do concurrent(k=1:nz, j=1:ny)
         ms%mass_flux_x_layer(1, j, k) = 0.0_wp
         ms%mass_flux_x_layer(nx + 1, j, k) = 0.0_wp
      end do
      ! Physical wall zeroing — see `continuity_zonal_flux` for the
      ! full rationale.  Must mirror the split form's wall closure or
      ! `test_split_zonal_only_matches_unsplit` breaks.
      do concurrent(k=1:nz, j=1:ny)
         ms%mass_flux_x_layer(grid%nghost + 1, j, k) = 0.0_wp
         ms%mass_flux_x_layer(grid%nghost + grid%nx_phys + 1, j, k) = 0.0_wp
      end do

      ! ---- Porous barriers (Adcroft 2013) ----
      ! Narrow the layer transport by the OPEN-AREA fraction of the face.
      ! Host-side gate: with the knob off there is no kernel launch and
      ! the loops above are textually unchanged, so the whole path is
      ! byte-identical.  Applied AFTER the wall zeroing (0 stays 0) and
      ! BEFORE the barotropic renormalisation, which must see the narrowed
      ! transports it is constraining.
      !
      ! WRITTEN INLINE, not as a call to `porous_narrow_3d`.  Handing
      ! `ms%mass_flux_*_layer` to an external subroutine as an
      ! `intent(inout)` actual makes nvfortran treat the array as ESCAPING,
      ! which pessimises every `do concurrent` in this routine — even
      ! though the branch never runs with the knob off.  Measured on a
      ! 600x600x50 default-path (porous OFF) double-gyre, V100: the call
      ! form costs `ocean_continuity` 10.58 s vs 9.11 s inline (+4.8% on
      ! total solver time vs origin/main; the inline form is +0.5%, i.e.
      ! noise).  `rdb_coriolis_adv` keeps the shared `porous_narrow_3d`
      ! helper — measured there at no cost.
      if (metrics%use_porous) then
         do concurrent(k=1:nz, j=1:ny, i=1:nx + 1)
            ms%mass_flux_x_layer(i, j, k) = ms%mass_flux_x_layer(i, j, k)* &
                                            metrics%por_face_area_u(i, j, k)
         end do
      end if

      ! ---- z-level closed faces: see the composition rule on
      ! `ocean_metrics_t%open_v`.  A SEPARATE pass, not composed into
      ! `por_face_area_u`: the two gates are independent and the porous
      ! fraction is refreshed per outer step while this mask is static.
      ! Inline for the same escaping-array reason as the porous pass.
      if (metrics%use_closed_faces) then
         do concurrent(k=1:nz, j=1:ny, i=1:nx + 1)
            ms%mass_flux_x_layer(i, j, k) = ms%mass_flux_x_layer(i, j, k)* &
                                            metrics%open_u(i, j, k)
         end do
      end if

      ! ============================================================
      ! Y-DIRECTION reconstruction
      ! ============================================================
      do concurrent(k=1:nz, j=3:ny - 2, i=1:nx) &
         local(dh_m1, dh_0, dh_p1, h_left, h_right, hm2, hm1, h0, hp1, hp2)
         h0 = ms%h_layer(i, j, k)
         hm1 = ppm_mirror_h(ms%h_layer(i, j - 1, k), h0, metrics%wet_T(i, j - 1))
         hp1 = ppm_mirror_h(ms%h_layer(i, j + 1, k), h0, metrics%wet_T(i, j + 1))
         hm2 = ppm_mirror_h(ms%h_layer(i, j - 2, k), hm1, metrics%wet_T(i, j - 2))
         hp2 = ppm_mirror_h(ms%h_layer(i, j + 2, k), hp1, metrics%wet_T(i, j + 2))
         call ppm_limited_slope(hm2, hm1, h0, dh_m1)
         call ppm_limited_slope(hm1, h0, hp1, dh_0)
         call ppm_limited_slope(h0, hp1, hp2, dh_p1)
         dh_0 = dh_0*metrics%wet_T(i, j - 1)*metrics%wet_T(i, j)*metrics%wet_T(i, j + 1)
         h_left = 0.5_wp*(hm1 + h0) - (dh_0 - dh_m1)/6.0_wp
         h_right = 0.5_wp*(h0 + hp1) - (dh_p1 - dh_0)/6.0_wp
         call ppm_cell_limiter(h0, h_left, h_right)
         if (do_pos) call ppm_limit_pos(h0, h_left, h_right, h_min_pos)
         this%h_face_right_y%data(i, j, k) = h_left
         this%h_face_left_y%data(i, j + 1, k) = h_right
      end do
      do concurrent(k=1:nz, i=1:nx)
         this%h_face_left_y%data(i, 1, k) = ms%h_layer(i, 1, k)
         this%h_face_right_y%data(i, 1, k) = ms%h_layer(i, 1, k)
         this%h_face_left_y%data(i, 2, k) = ms%h_layer(i, 1, k)
         this%h_face_right_y%data(i, 2, k) = ms%h_layer(i, 2, k)
         this%h_face_left_y%data(i, 3, k) = ms%h_layer(i, 2, k)
         this%h_face_right_y%data(i, 3, k) = ms%h_layer(i, 2, k)
         this%h_face_right_y%data(i, ny - 1, k) = ms%h_layer(i, ny - 1, k)
         this%h_face_left_y%data(i, ny, k) = ms%h_layer(i, ny - 1, k)
         this%h_face_right_y%data(i, ny, k) = ms%h_layer(i, ny, k)
         this%h_face_left_y%data(i, ny + 1, k) = ms%h_layer(i, ny, k)
         this%h_face_right_y%data(i, ny + 1, k) = ms%h_layer(i, ny, k)
      end do

      do concurrent(k=1:nz, j=2:ny, i=1:nx) local(v, h_face)
         v = ms%v_face_y_layer(i, j, k)
         if (v >= 0.0_wp) then
            h_face = this%h_face_left_y%data(i, j, k)
         else
            h_face = this%h_face_right_y%data(i, j, k)
         end if
         ms%mass_flux_y_layer(i, j, k) = v*h_face*metrics%dx_cv(i, j)
      end do
      do concurrent(k=1:nz, i=1:nx)
         ms%mass_flux_y_layer(i, 1, k) = 0.0_wp
         ms%mass_flux_y_layer(i, ny + 1, k) = 0.0_wp
      end do
      ! Physical wall zeroing — mirrors split form's
      ! `continuity_meridional_flux`.
      do concurrent(k=1:nz, i=1:nx)
         ms%mass_flux_y_layer(i, grid%nghost + 1, k) = 0.0_wp
         ms%mass_flux_y_layer(i, grid%nghost + grid%ny_phys + 1, k) = 0.0_wp
      end do

      ! ---- Porous barriers: see the zonal twin (incl. why it is inline) ----
      if (metrics%use_porous) then
         do concurrent(k=1:nz, j=1:ny + 1, i=1:nx)
            ms%mass_flux_y_layer(i, j, k) = ms%mass_flux_y_layer(i, j, k)* &
                                            metrics%por_face_area_v(i, j, k)
         end do
      end if

      ! ---- z-level closed faces: see the zonal twin ----
      if (metrics%use_closed_faces) then
         do concurrent(k=1:nz, j=1:ny + 1, i=1:nx)
            ms%mass_flux_y_layer(i, j, k) = ms%mass_flux_y_layer(i, j, k)* &
                                            metrics%open_v(i, j, k)
         end do
      end if

      ! ============================================================
      ! Per-layer flux divergence (transport divergence · iareaT)
      ! ============================================================
      do concurrent(k=1:nz, j=1:ny, i=1:nx)
         ms%flux_h_layer(i, j, k) = &
            ((ms%mass_flux_x_layer(i + 1, j, k) - &
              ms%mass_flux_x_layer(i, j, k)) + &
             (ms%mass_flux_y_layer(i, j + 1, k) - &
              ms%mass_flux_y_layer(i, j, k)))*metrics%iareaT(i, j)
      end do
   end subroutine continuity_compute_fluxes