tracer_advect_one_impl Subroutine

   pure subroutine tracer_advect_one_impl(nx, ny, nz, dt, iareaT, wet_T, h, hTr, &
                                          mass_flux_x, mass_flux_y, &
                                          Tr_face_left_x, Tr_face_right_x, &
                                          Tr_face_left_y, Tr_face_right_y)
      !! One-tracer PPM advection.  Flat-impl: takes bare 3D arrays
      !! (no derived-type deref inside do-concurrent), so NVHPC
      !! stdpar handles it cleanly even for tracers stored in an
      !! array-of-derived-types registry.
      !!
      !! Three passes per direction:
      !!   1. PPM reconstruction on Tr = hTr/h (Colella-Woodward
      !!      monotonic limiter) — same helpers as continuity-PPM.
      !!      Stores left/right face states in the scratch buffers.
      !!   2. Upwind pick using mass_flux sign — overwrites the same
      !!      scratch buffer with mass_flux * Tr_face (the tracer
      !!      mass flux at each face).
      !!   3. Flux divergence + forward-Euler update of hTr.
      !!
      !! Wall faces are gated by mass_flux (= 0 at walls from
      !! continuity), so tracer flux through walls is identically
      !! zero — no special wall handling needed in this kernel.
      integer, intent(in) :: nx, ny, nz
      real(wp), intent(in) :: dt
      real(wp), intent(in) :: iareaT(nx, ny)
      real(wp), intent(in) :: wet_T(nx, ny)
      real(wp), intent(in) :: h(nx, ny, nz)
      real(wp), intent(inout) :: hTr(nx, ny, nz)
      real(wp), intent(in) :: mass_flux_x(nx + 1, ny, nz)
      real(wp), intent(in) :: mass_flux_y(nx, ny + 1, nz)
      real(wp), intent(inout) :: Tr_face_left_x(nx + 1, ny, nz)
      real(wp), intent(inout) :: Tr_face_right_x(nx + 1, ny, nz)
      real(wp), intent(inout) :: Tr_face_left_y(nx, ny + 1, nz)
      real(wp), intent(inout) :: Tr_face_right_y(nx, ny + 1, nz)

      integer :: i, j, k
      real(wp) :: Tr_m2, Tr_m1, Tr_0, Tr_p1, Tr_p2
      real(wp) :: dh_m1, dh_0, dh_p1, Tr_left, Tr_right
      real(wp) :: mass_x, mass_y, Tr_face

      ! ===== X-direction PPM reconstruction =====
      ! Interior cells with full 5-point stencil.  Mirror-T at land
      ! neighbours (C2); bit-identical for all-wet.
      do concurrent(k=1:nz, j=1:ny, i=3:nx - 2) &
         local(Tr_m2, Tr_m1, Tr_0, Tr_p1, Tr_p2, &
               dh_m1, dh_0, dh_p1, Tr_left, Tr_right)
         Tr_0 = hTr(i, j, k)/h(i, j, k)
         Tr_m1 = ppm_mirror_h(hTr(i - 1, j, k)/h(i - 1, j, k), Tr_0, wet_T(i - 1, j))
         Tr_p1 = ppm_mirror_h(hTr(i + 1, j, k)/h(i + 1, j, k), Tr_0, wet_T(i + 1, j))
         Tr_m2 = ppm_mirror_h(hTr(i - 2, j, k)/h(i - 2, j, k), Tr_m1, wet_T(i - 2, j))
         Tr_p2 = ppm_mirror_h(hTr(i + 2, j, k)/h(i + 2, j, k), Tr_p1, wet_T(i + 2, j))
         call ppm_limited_slope(Tr_m2, Tr_m1, Tr_0, dh_m1)
         call ppm_limited_slope(Tr_m1, Tr_0, Tr_p1, dh_0)
         call ppm_limited_slope(Tr_0, Tr_p1, Tr_p2, dh_p1)
         dh_0 = dh_0*wet_T(i - 1, j)*wet_T(i, j)*wet_T(i + 1, j)
         Tr_left = 0.5_wp*(Tr_m1 + Tr_0) - (dh_0 - dh_m1)/6.0_wp
         Tr_right = 0.5_wp*(Tr_0 + Tr_p1) - (dh_p1 - dh_0)/6.0_wp
         call ppm_cell_limiter(Tr_0, Tr_left, Tr_right)
         Tr_face_right_x(i, j, k) = Tr_left
         Tr_face_left_x(i + 1, j, k) = Tr_right
      end do
      ! Boundary cells (i=1, 2, nx-1, nx): 1st-order fallback (face
      ! value = abutting cell's concentration).  Mass flux is zero
      ! at i=1 and i=nx+1 walls anyway, so the wall faces' Tr value
      ! doesn't enter the divergence.
      do concurrent(k=1:nz, j=1:ny)
         Tr_face_left_x(1, j, k) = hTr(1, j, k)/h(1, j, k)
         Tr_face_right_x(1, j, k) = hTr(1, j, k)/h(1, j, k)
         Tr_face_left_x(2, j, k) = hTr(1, j, k)/h(1, j, k)
         Tr_face_right_x(2, j, k) = hTr(2, j, k)/h(2, j, k)
         ! Face 3: cell 2's 5-point stencil is short, so its right
         ! edge falls back to 1st order — match continuity's fix.
         Tr_face_left_x(3, j, k) = hTr(2, j, k)/h(2, j, k)
         Tr_face_right_x(3, j, k) = hTr(2, j, k)/h(2, j, k)
         ! Face nx-1: mirror of face 3.
         Tr_face_right_x(nx - 1, j, k) = hTr(nx - 1, j, k)/h(nx - 1, j, k)
         Tr_face_left_x(nx, j, k) = hTr(nx - 1, j, k)/h(nx - 1, j, k)
         Tr_face_right_x(nx, j, k) = hTr(nx, j, k)/h(nx, j, k)
         Tr_face_left_x(nx + 1, j, k) = hTr(nx, j, k)/h(nx, j, k)
         Tr_face_right_x(nx + 1, j, k) = hTr(nx, j, k)/h(nx, j, k)
      end do

      ! Upwind pick at east faces -> overwrite Tr_face_left_x with
      ! the per-face tracer mass flux (mass_flux_x * Tr_face_upwind).
      do concurrent(k=1:nz, j=1:ny, i=1:nx + 1) local(mass_x, Tr_face)
         mass_x = mass_flux_x(i, j, k)
         if (mass_x >= 0.0_wp) then
            Tr_face = Tr_face_left_x(i, j, k)
         else
            Tr_face = Tr_face_right_x(i, j, k)
         end if
         Tr_face_left_x(i, j, k) = mass_x*Tr_face
      end do

      ! ===== Y-direction PPM reconstruction =====
      do concurrent(k=1:nz, j=3:ny - 2, i=1:nx) &
         local(Tr_m2, Tr_m1, Tr_0, Tr_p1, Tr_p2, &
               dh_m1, dh_0, dh_p1, Tr_left, Tr_right)
         Tr_0 = hTr(i, j, k)/h(i, j, k)
         Tr_m1 = ppm_mirror_h(hTr(i, j - 1, k)/h(i, j - 1, k), Tr_0, wet_T(i, j - 1))
         Tr_p1 = ppm_mirror_h(hTr(i, j + 1, k)/h(i, j + 1, k), Tr_0, wet_T(i, j + 1))
         Tr_m2 = ppm_mirror_h(hTr(i, j - 2, k)/h(i, j - 2, k), Tr_m1, wet_T(i, j - 2))
         Tr_p2 = ppm_mirror_h(hTr(i, j + 2, k)/h(i, j + 2, k), Tr_p1, wet_T(i, j + 2))
         call ppm_limited_slope(Tr_m2, Tr_m1, Tr_0, dh_m1)
         call ppm_limited_slope(Tr_m1, Tr_0, Tr_p1, dh_0)
         call ppm_limited_slope(Tr_0, Tr_p1, Tr_p2, dh_p1)
         dh_0 = dh_0*wet_T(i, j - 1)*wet_T(i, j)*wet_T(i, j + 1)
         Tr_left = 0.5_wp*(Tr_m1 + Tr_0) - (dh_0 - dh_m1)/6.0_wp
         Tr_right = 0.5_wp*(Tr_0 + Tr_p1) - (dh_p1 - dh_0)/6.0_wp
         call ppm_cell_limiter(Tr_0, Tr_left, Tr_right)
         Tr_face_right_y(i, j, k) = Tr_left
         Tr_face_left_y(i, j + 1, k) = Tr_right
      end do
      do concurrent(k=1:nz, i=1:nx)
         Tr_face_left_y(i, 1, k) = hTr(i, 1, k)/h(i, 1, k)
         Tr_face_right_y(i, 1, k) = hTr(i, 1, k)/h(i, 1, k)
         Tr_face_left_y(i, 2, k) = hTr(i, 1, k)/h(i, 1, k)
         Tr_face_right_y(i, 2, k) = hTr(i, 2, k)/h(i, 2, k)
         Tr_face_left_y(i, 3, k) = hTr(i, 2, k)/h(i, 2, k)
         Tr_face_right_y(i, 3, k) = hTr(i, 2, k)/h(i, 2, k)
         Tr_face_right_y(i, ny - 1, k) = hTr(i, ny - 1, k)/h(i, ny - 1, k)
         Tr_face_left_y(i, ny, k) = hTr(i, ny - 1, k)/h(i, ny - 1, k)
         Tr_face_right_y(i, ny, k) = hTr(i, ny, k)/h(i, ny, k)
         Tr_face_left_y(i, ny + 1, k) = hTr(i, ny, k)/h(i, ny, k)
         Tr_face_right_y(i, ny + 1, k) = hTr(i, ny, k)/h(i, ny, k)
      end do

      ! Upwind pick at north faces
      do concurrent(k=1:nz, j=1:ny + 1, i=1:nx) local(mass_y, Tr_face)
         mass_y = mass_flux_y(i, j, k)
         if (mass_y >= 0.0_wp) then
            Tr_face = Tr_face_left_y(i, j, k)
         else
            Tr_face = Tr_face_right_y(i, j, k)
         end if
         Tr_face_left_y(i, j, k) = mass_y*Tr_face
      end do

      ! ===== Apply: forward-Euler hTr update (transport div · iareaT) =====
      do concurrent(k=1:nz, j=1:ny, i=1:nx)
         hTr(i, j, k) = hTr(i, j, k) - dt*( &
                        (Tr_face_left_x(i + 1, j, k) - Tr_face_left_x(i, j, k)) + &
                        (Tr_face_left_y(i, j + 1, k) - Tr_face_left_y(i, j, k)))*iareaT(i, j)
      end do
   end subroutine tracer_advect_one_impl