ice_cat_flux_y_impl Subroutine

   pure subroutine ice_cat_flux_y_impl(wet_T, dx_cv, idyT, v_ice, mca, htot_work, hl_y_work, &
                                       hr_y_work, vhtot_work, vh_out, dt_adv, nghost, ny_phys, &
                                       ncat, nx, ny, wall_s, wall_n)
      !! Meridional twin of `ice_cat_flux_x_impl`.  Same face-indexed edge
      !! STORAGE + swept orientation as `continuity_meridional_flux`
      !! (`rdb_continuity.F90:1159-1202`): `hl_y_work(i,j)` == north-face
      !! `h_face_left_y` (from the SOUTH cell `j-1`), `hr_y_work(i,j)` ==
      !! `h_face_right_y` (from the NORTH cell `j`).  CFL uses
      !! `IdyT(donor)` (SIS2 `vol_CFL=.false.` default).
      integer, intent(in) :: nghost, ny_phys, ncat, nx, ny
      real(wp), intent(in) :: wet_T(nx, ny)
      real(wp), intent(in) :: dx_cv(nx, ny + 1)
      real(wp), intent(in) :: idyT(nx, ny)
      real(wp), intent(in) :: v_ice(nx, ny + 1)
      real(wp), intent(in) :: mca(nx, ny, ncat)
      real(wp), intent(inout) :: htot_work(nx, ny)
      real(wp), intent(inout) :: hl_y_work(nx, ny + 1)
      real(wp), intent(inout) :: hr_y_work(nx, ny + 1)
      real(wp), intent(inout) :: vhtot_work(nx, ny + 1)
      real(wp), intent(out) :: vh_out(nx, ny + 1, ncat)
      real(wp), intent(in) :: dt_adv
      logical, intent(in) :: wall_s, wall_n
         !! Zero face `nghost+1` / `nghost+ny_phys+1`: only on a physical,
         !! non-periodic tile edge (F1, see the x twin).

      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) :: v, cfl, dh_d, curv3_d, h_face, i_htot

      do concurrent(j=1:ny, i=1:nx)
         htot_work(i, j) = sum(mca(i, j, :))
      end do

      do concurrent(j=3:ny - 2, i=1:nx) &
         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, j - 1), h0, wet_T(i, j - 1))
         hp1 = ppm_mirror_h(htot_work(i, j + 1), h0, wet_T(i, j + 1))
         hm2 = ppm_mirror_h(htot_work(i, j - 2), hm1, wet_T(i, j - 2))
         hp2 = ppm_mirror_h(htot_work(i, j + 2), hp1, 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*wet_T(i, j - 1)*wet_T(i, j)*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)
         call ppm_limit_pos(h0, h_left, h_right, 0.0_wp)
         hr_y_work(i, j) = h_left
         hl_y_work(i, j + 1) = h_right
      end do
      ! Edge fallback — NOT `hr_y_work(:, 3)` (F2, see the x twin).
      do concurrent(i=1:nx)
         hl_y_work(i, 1) = htot_work(i, 1)
         hr_y_work(i, 1) = htot_work(i, 1)
         hl_y_work(i, 2) = htot_work(i, 1)
         hr_y_work(i, 2) = htot_work(i, 2)
         hl_y_work(i, 3) = htot_work(i, 2)
         hr_y_work(i, ny - 1) = htot_work(i, ny - 1)
         hl_y_work(i, ny) = htot_work(i, ny - 1)
         hr_y_work(i, ny) = htot_work(i, ny)
         hl_y_work(i, ny + 1) = htot_work(i, ny)
         hr_y_work(i, ny + 1) = htot_work(i, ny)
      end do

      ! Donor-edge orientation from `continuity_meridional_flux`:
      !   v>0 (donor j-1): edge = hl_y_work(i,j); dh = hr_y_work(i,j-1)
      !     - hl_y_work(i,j); curv3 = hr_y_work(i,j-1) + hl_y_work(i,j)
      !     - 2*htot(i,j-1).
      !   v<0 (donor j):   edge = hr_y_work(i,j); dh = hl_y_work(i,j+1)
      !     - hr_y_work(i,j); curv3 = hr_y_work(i,j) + hl_y_work(i,j+1)
      !     - 2*htot(i,j).
      ! No `max(...,0)` clamp (see the x-pass twin).
      do concurrent(j=2:ny, i=1:nx) local(v, cfl, dh_d, curv3_d, h_face)
         v = v_ice(i, j)
         if (v > 0.0_wp) then
            cfl = v*dt_adv*idyT(i, j - 1)
            h_face = hl_y_work(i, j)
            dh_d = hr_y_work(i, j - 1) - h_face
            curv3_d = hr_y_work(i, j - 1) + h_face - 2.0_wp*htot_work(i, j - 1)
            vhtot_work(i, j) = dx_cv(i, j)*v*volcfl_face(h_face, dh_d, curv3_d, cfl)
         else if (v < 0.0_wp) then
            cfl = (-v)*dt_adv*idyT(i, j)
            h_face = hr_y_work(i, j)
            dh_d = hl_y_work(i, j + 1) - h_face
            curv3_d = h_face + hl_y_work(i, j + 1) - 2.0_wp*htot_work(i, j)
            vhtot_work(i, j) = dx_cv(i, j)*v*volcfl_face(h_face, dh_d, curv3_d, cfl)
         else
            vhtot_work(i, j) = 0.0_wp
         end if
      end do
      do concurrent(i=1:nx)
         vhtot_work(i, 1) = 0.0_wp
         vhtot_work(i, ny + 1) = 0.0_wp
      end do
      if (wall_s) then
         do concurrent(i=1:nx)
            vhtot_work(i, nghost + 1) = 0.0_wp
         end do
      end if
      if (wall_n) then
         do concurrent(i=1:nx)
            vhtot_work(i, nghost + ny_phys + 1) = 0.0_wp
         end do
      end if

      do concurrent(j=1:ny + 1, i=1:nx, c=1:ncat) local(i_htot)
         if (vhtot_work(i, j) == 0.0_wp) then
            vh_out(i, j, c) = 0.0_wp
         else if (v_ice(i, j) >= 0.0_wp) then
            if (htot_work(i, j - 1) > 0.0_wp) then
               i_htot = 1.0_wp/htot_work(i, j - 1)
            else
               i_htot = 0.0_wp
            end if
            vh_out(i, j, c) = vhtot_work(i, j)*mca(i, j - 1, 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
            vh_out(i, j, c) = vhtot_work(i, j)*mca(i, j, c)*i_htot
         end if
      end do
   end subroutine ice_cat_flux_y_impl