ice_cat_flux_x_impl Subroutine

   pure subroutine ice_cat_flux_x_impl(wet_T, dy_cu, idxT, u_ice, mca, htot_work, hl_x_work, &
                                       hr_x_work, uhtot_work, uh_out, dt_adv, nghost, nx_phys, &
                                       ncat, nx, ny, wall_w, wall_e)
      !! SIS2 `zonal_mass_flux` (`SIS_continuity.F90:1064`): PPM
      !! reconstruction of the category-SUMMED mass `htot`, ONE total
      !! face transport `uhtot` from the swept-volume parabola integral
      !! (`volcfl_face`, bit-for-bit the SIS2 face expression), then the
      !! PROPORTIONATE split `uh(c) = uhtot*mca(donor,c)*I_htot(donor)`
      !! (`SIS_continuity.F90:1199-1205`, Adcroft reciprocal — `I_htot=0`
      !! when `htot(donor)<=0`).  Stencil + edge STORAGE CONVENTION +
      !! swept-face orientation copied verbatim from `continuity_zonal_flux`
      !! (`rdb_continuity.F90:871-948`): `hl_x_work(i)` == `h_face_left_x(i)`
      !! is the value AT east face `i` from the LEFT cell `i-1` (that
      !! cell's OWN downwind edge); `hr_x_work(i)` == `h_face_right_x(i)`
      !! is from the RIGHT cell `i` (its OWN left edge).  H3-style limited
      !! edges + CW84 + `ppm_limit_pos` (SIS2 runs `PPM_limit_pos`
      !! UNCONDITIONALLY on the PD scheme, so it is not optional here).
      !! CFL metric: SIS2's shipped default `vol_CFL=.false.` uses
      !! `CFL = |u|*dt*IdxT(donor)` (`SIS_continuity.F90:1165`) — the
      !! T-cell inverse spacing, NOT the `dy_cu*iareaT` swept-area ratio
      !! (== SIS2's `vol_CFL=.true.` variant).  Identical on uniform
      !! Cartesian; correct on spherical / anisotropic.
      integer, intent(in) :: nghost, nx_phys, ncat, nx, ny
      real(wp), intent(in) :: wet_T(nx, ny)
      real(wp), intent(in) :: dy_cu(nx + 1, ny)
      real(wp), intent(in) :: idxT(nx, ny)
      real(wp), intent(in) :: u_ice(nx + 1, ny)
      real(wp), intent(in) :: mca(nx, ny, ncat)
      real(wp), intent(inout) :: htot_work(nx, ny)
      real(wp), intent(inout) :: hl_x_work(nx + 1, ny)
      real(wp), intent(inout) :: hr_x_work(nx + 1, ny)
      real(wp), intent(inout) :: uhtot_work(nx + 1, ny)
      real(wp), intent(out) :: uh_out(nx + 1, ny, ncat)
      real(wp), intent(in) :: dt_adv
      logical, intent(in) :: wall_w, wall_e
         !! Zero face `nghost+1` / `nghost+nx_phys+1`: only on a physical,
         !! non-periodic tile edge.  At an MPI or periodic seam the face is
         !! an ordinary face whose donor ghosts X4 filled (F1).

      integer :: i, j, c
      real(wp) :: dh_m1, dh_0, dh_p1, h_left, h_right
      real(wp) :: hm2, hm1, h0, hp1, hp2
      real(wp) :: u, cfl, dh_d, curv3_d, h_face, i_htot

      ! ---- Category-summed mass ----
      do concurrent(j=1:ny, i=1:nx)
         htot_work(i, j) = sum(mca(i, j, :))
      end do

      ! ---- PPM reconstruction (interior 5-point stencil) ----
      ! Cell i's LEFT edge -> hr_x_work(i) (right-cell state at face i);
      ! cell i's RIGHT edge -> hl_x_work(i+1) (left-cell state at face
      ! i+1).  Exactly `continuity_compute_fluxes`'s write pattern.
      do concurrent(j=1:ny, i=3:nx - 2) &
         local(dh_m1, dh_0, dh_p1, h_left, h_right, hm2, hm1, h0, hp1, hp2)
         h0 = htot_work(i, j)
         hm1 = ppm_mirror_h(htot_work(i - 1, j), h0, wet_T(i - 1, j))
         hp1 = ppm_mirror_h(htot_work(i + 1, j), h0, wet_T(i + 1, j))
         hm2 = ppm_mirror_h(htot_work(i - 2, j), hm1, wet_T(i - 2, j))
         hp2 = ppm_mirror_h(htot_work(i + 2, j), hp1, 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*wet_T(i - 1, j)*wet_T(i, j)*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)
         call ppm_limit_pos(h0, h_left, h_right, 0.0_wp)
         hr_x_work(i, j) = h_left
         hl_x_work(i + 1, j) = h_right
      end do
      ! First-order fallback for the edges the 5-point loop above cannot
      ! reach (cells 1, 2, nx-1, nx; face-shaped writes, in bounds:
      ! `nx+1` is the last valid face).  NOT `hr_x_work(3)`: cell 3 is
      ! the loop's first cell, so its left edge is already a full PPM
      ! value — and with `nghost = 3` it is the donor-side edge the
      ! `u > 0` flux at the west tile-edge face `nghost+1` reads, so
      ! overwriting it made an MPI seam face first-order where the serial
      ! run's interior face is PPM (F2; the y twin is the same).
      do concurrent(j=1:ny)
         hl_x_work(1, j) = htot_work(1, j)
         hr_x_work(1, j) = htot_work(1, j)
         hl_x_work(2, j) = htot_work(1, j)
         hr_x_work(2, j) = htot_work(2, j)
         hl_x_work(3, j) = htot_work(2, j)
         hr_x_work(nx - 1, j) = htot_work(nx - 1, j)
         hl_x_work(nx, j) = htot_work(nx - 1, j)
         hr_x_work(nx, j) = htot_work(nx, j)
         hl_x_work(nx + 1, j) = htot_work(nx, j)
         hr_x_work(nx + 1, j) = htot_work(nx, j)
      end do

      ! ---- Total-mass face transport (swept-volume PPM integral) ----
      ! Donor-edge orientation copied verbatim from
      ! `continuity_zonal_flux` (`rdb_continuity.F90:923-948`):
      !   u>0 (donor i-1): edge = hl_x_work(i) [cell i-1's downwind
      !     edge]; dh = hr_x_work(i-1) - hl_x_work(i); curv3 =
      !     hr_x_work(i-1) + hl_x_work(i) - 2*htot(i-1).
      !   u<0 (donor i):   edge = hr_x_work(i) [cell i's downwind edge];
      !     dh = hl_x_work(i+1) - hr_x_work(i); curv3 = hr_x_work(i) +
      !     hl_x_work(i+1) - 2*htot(i).
      ! No `max(...,0)` clamp: a `ppm_limit_pos`'d parabola has a
      ! non-negative swept mean for CFL <= 1, so `volcfl_face` stays >= 0
      ! here — the earlier negative values were the SIGNATURE of the
      ! wrong (now-fixed) donor-edge stencil feeding mixed-cell inputs,
      ! not a property of the scheme (neither SIS2 nor
      ! continuity_zonal_flux clamps).
      do concurrent(j=1:ny, i=2:nx) local(u, cfl, dh_d, curv3_d, h_face)
         u = u_ice(i, j)
         if (u > 0.0_wp) then
            cfl = u*dt_adv*idxT(i - 1, j)
            h_face = hl_x_work(i, j)
            dh_d = hr_x_work(i - 1, j) - h_face
            curv3_d = hr_x_work(i - 1, j) + h_face - 2.0_wp*htot_work(i - 1, j)
            uhtot_work(i, j) = dy_cu(i, j)*u*volcfl_face(h_face, dh_d, curv3_d, cfl)
         else if (u < 0.0_wp) then
            cfl = (-u)*dt_adv*idxT(i, j)
            h_face = hr_x_work(i, j)
            dh_d = hl_x_work(i + 1, j) - h_face
            curv3_d = h_face + hl_x_work(i + 1, j) - 2.0_wp*htot_work(i, j)
            uhtot_work(i, j) = dy_cu(i, j)*u*volcfl_face(h_face, dh_d, curv3_d, cfl)
         else
            uhtot_work(i, j) = 0.0_wp
         end if
      end do
      do concurrent(j=1:ny)
         uhtot_work(1, j) = 0.0_wp
         uhtot_work(nx + 1, j) = 0.0_wp
      end do
      ! Physical wall faces (not just array edges) — mirrors
      ! `continuity_zonal_flux`'s wall-zeroing rationale: a moving ocean
      ! surface layer can leave a nonzero sampled velocity at the
      ! interior physical wall even on a closed-boundary configuration.
      ! Only a PHYSICAL, non-periodic tile edge is a wall (F1).
      if (wall_w) then
         do concurrent(j=1:ny)
            uhtot_work(nghost + 1, j) = 0.0_wp
         end do
      end if
      if (wall_e) then
         do concurrent(j=1:ny)
            uhtot_work(nghost + nx_phys + 1, j) = 0.0_wp
         end do
      end if

      ! ---- Proportionate category split (Adcroft reciprocal) ----
      do concurrent(j=1:ny, i=1:nx + 1, c=1:ncat) local(i_htot)
         if (uhtot_work(i, j) == 0.0_wp) then
            uh_out(i, j, c) = 0.0_wp
         else if (u_ice(i, j) >= 0.0_wp) then
            if (htot_work(i - 1, j) > 0.0_wp) then
               i_htot = 1.0_wp/htot_work(i - 1, j)
            else
               i_htot = 0.0_wp
            end if
            uh_out(i, j, c) = uhtot_work(i, j)*mca(i - 1, j, c)*i_htot
         else
            if (htot_work(i, j) > 0.0_wp) then
               i_htot = 1.0_wp/htot_work(i, j)
            else
               i_htot = 0.0_wp
            end if
            uh_out(i, j, c) = uhtot_work(i, j)*mca(i, j, c)*i_htot
         end if
      end do
   end subroutine ice_cat_flux_x_impl