continuity_meridional_flux Subroutine

   pure subroutine continuity_meridional_flux(grid, metrics, this, ms, dt, vhbt, bc, &
                                              visc_rem, v_cor)
      !! Meridional (y-only) PPM reconstruction + per-face mass
      !! flux.  Mirror of `continuity_zonal_flux`, with the same
      !! optional `vhbt` transport-constraint renormalisation.
      !! Writes `ms%mass_flux_y_layer`.  Walls at j=1 and j=ny+1
      !! zeroed.  In the Lie split this runs *after* the zonal
      !! apply, so it reconstructs against the already-updated
      !! `ms%h_layer`.
      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
      real(wp), intent(in) :: dt
         !! Time increment (s).  Used only for the swept-volume CFL when
         !! `this%vol_cfl = .true.`; ignored (bit-identical) otherwise.
      real(wp), intent(in), optional :: vhbt(:, :)
         !! Time-mean north-face transport from barotropic substep (m³/s),
         !! shape `(nx, ny+1)`.  Width-weighted (carries `dx_cv`).
      type(ocean_bc_state_t), intent(in), optional :: bc
      ! assumed-shape-ok: pure passthrough to the renormaliser.
      real(wp), intent(in), optional :: visc_rem(:, :, :)
         !! Per-layer viscous remnant gamma_k. Absent => 1, bit-identical.
      ! assumed-shape-ok: pure passthrough to the renormaliser.
      real(wp), intent(inout), optional :: v_cor(:, :, :)
         !! Forwarded MOM6 `v_cor` destination — separate time-mean field.

      integer :: i, j, k, nx, ny, nz
      integer :: bc_s_tag, bc_n_tag
      real(wp) :: dh_m1, dh_0, dh_p1, h_left, h_right, v, h_face
      real(wp) :: hm2, hm1, h0, hp1, hp2
      real(wp) :: cfl_d, curv3_d, dh_d
      logical :: do_pos, do_volcfl, do_pd
      logical :: has_s_flux, has_n_flux, renorm_skip_walls
      real(wp) :: h_min_pos, two_h_lim

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      do_pos = this%use_ppm_limit_pos
      do_volcfl = this%vol_cfl
      h_min_pos = this%h_min
      ! P1 positive-definite reconstruction floor (MOM6's positive-definite
      ! PPM).  Scalars copied to locals so the DC loop never walks the
      ! continuity_t descriptor.  Off ⇒ untaken branch ⇒ bit-identical.
      do_pd = this%positive_definite
      two_h_lim = 2.0_wp*this%h_lim

      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)
         if (do_pd) then
            h_left = max(h_left, two_h_lim)
            h_right = max(h_right, two_h_lim)
         end if
         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, cfl_d, curv3_d, dh_d)
         v = ms%v_face_y_layer(i, j, k)
         if (v >= 0.0_wp) then
            ! Donor = cell j-1; downwind (north) edge = h_face_left_y(j).
            h_face = this%h_face_left_y%data(i, j, k)
            if (do_volcfl) then
               ! Swept-oriented donor edge diff dh = h_S - h_N (south-north).
               dh_d = this%h_face_right_y%data(i, j - 1, k) - h_face
               curv3_d = this%h_face_right_y%data(i, j - 1, k) + h_face &
                         - 2.0_wp*ms%h_layer(i, j - 1, k)
               cfl_d = v*dt*metrics%dx_cv(i, j)*metrics%iareaT(i, j - 1)
               h_face = volcfl_face(h_face, dh_d, curv3_d, cfl_d)
            end if
         else
            ! Donor = cell j; downwind (south) edge = h_face_right_y(j).
            h_face = this%h_face_right_y%data(i, j, k)
            if (do_volcfl) then
               ! Swept-oriented donor edge diff dh = h_N - h_S (north-south).
               dh_d = this%h_face_left_y%data(i, j + 1, k) - h_face
               curv3_d = h_face + this%h_face_left_y%data(i, j + 1, k) &
                         - 2.0_wp*ms%h_layer(i, j, k)
               cfl_d = (-v)*dt*metrics%dx_cv(i, j)*metrics%iareaT(i, j)
               h_face = volcfl_face(h_face, dh_d, curv3_d, cfl_d)
            end if
         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 — see `continuity_zonal_flux` for the
      ! detailed rationale.  Without this, v_face_y at the physical
      ! south/north walls picks up Coriolis / wind contributions and
      ! the PPM transport leaks tracer through to ghost cells.
      ! OBC dispatch matches the zonal helper — see header.
      ! Under MPI decomposition bc%has_south / bc%has_north gates the
      ! zeroing so a subdomain seam face is left for the halo exchange.
      bc_s_tag = OBC_WALL
      bc_n_tag = OBC_WALL
      ! Physical-edge flags: default .true. => single-rank bit-identical
      ! behaviour; .false. at an MPI seam skips the hard zero (O0, D0).
      has_s_flux = .true.
      has_n_flux = .true.
      if (present(bc)) then
         bc_s_tag = ocean_bc_outer_face_tag(bc%south%bc_type)
         bc_n_tag = ocean_bc_outer_face_tag(bc%north%bc_type)
         has_s_flux = bc%has_south
         has_n_flux = bc%has_north
      end if
      do concurrent(k=1:nz, i=1:nx)
         if (bc_s_tag == OBC_WALL .and. has_s_flux) ms%mass_flux_y_layer(i, grid%nghost + 1, k) = 0.0_wp
         if (bc_n_tag == OBC_WALL .and. has_n_flux) 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

      ! MOM6-style transport constraint — see continuity_zonal_flux
      ! for the rationale.
      if (present(vhbt)) then
         ! Any non-WALL meridional edge (periodic OR open-class) carries
         ! genuine transport at the physical-wall face and must be
         ! renormalised — see continuity_zonal_flux for the full rationale.
         ! Genuine closed WALL is skipped (default); all-WALL stays
         ! bit-identical to pre-OBC behaviour.
         ! One call site per `por` actual — see the zonal twin for why
         ! folding the two BC branches together is exact and why it
         ! matters for the default-path cost.
         renorm_skip_walls = (bc_s_tag == OBC_WALL .and. bc_n_tag == OBC_WALL)
         ! Four branches — see the zonal twin for why.
         if (metrics%use_porous) then
            if (metrics%use_closed_faces) then
               call renormalise_meridional_flux_to_vhbt(grid, metrics, this, ms, vhbt, dt, &
                                                        skip_walls=renorm_skip_walls, &
                                                        has_south=has_s_flux, has_north=has_n_flux, &
                                                        visc_rem=visc_rem, v_cor=v_cor, &
                                                        use_por=.true., por=metrics%por_face_area_v, &
                                                        use_open=.true., open_f=metrics%open_v)
            else
               call renormalise_meridional_flux_to_vhbt(grid, metrics, this, ms, vhbt, dt, &
                                                        skip_walls=renorm_skip_walls, &
                                                        has_south=has_s_flux, has_north=has_n_flux, &
                                                        visc_rem=visc_rem, v_cor=v_cor, &
                                                        use_por=.true., por=metrics%por_face_area_v, &
                                                        use_open=.false., open_f=this%h_face_right_y%data)
            end if
         else if (metrics%use_closed_faces) then
            call renormalise_meridional_flux_to_vhbt(grid, metrics, this, ms, vhbt, dt, &
                                                     skip_walls=renorm_skip_walls, &
                                                     has_south=has_s_flux, has_north=has_n_flux, &
                                                     visc_rem=visc_rem, v_cor=v_cor, &
                                                     use_por=.false., por=this%h_face_left_y%data, &
                                                     use_open=.true., open_f=metrics%open_v)
         else
            call renormalise_meridional_flux_to_vhbt(grid, metrics, this, ms, vhbt, dt, &
                                                     skip_walls=renorm_skip_walls, &
                                                     has_south=has_s_flux, has_north=has_n_flux, &
                                                     visc_rem=visc_rem, v_cor=v_cor, &
                                                     use_por=.false., por=this%h_face_left_y%data, &
                                                     use_open=.false., open_f=this%h_face_right_y%data)
         end if
      end if
   end subroutine continuity_meridional_flux