!> @file boundary_2d.f90
!> @brief 2D boundary halo filling: doubly-periodic wrap plus per-edge physical
!! BCs (reflecting wall / outflow / supersonic inlet), normal-aware on
!! curvilinear grids, with an MPI dual-path for inter-rank edges.
module boundary_2d
  use precision, only: wp
  use solver_state_2d, only: solver_state_2d_t
  use option_registry, only: bc_periodic, bc_outflow, bc_reflecting, bc_supersonic_inlet, &
                             bc_farfield, method_fdm
  use euler_physics_2d, only: primitives_2d
  implicit none
  private
  public :: apply_periodic_halos_2d, apply_halos_2d, wall_edges_2d

contains

  !> Report which of the four block edges are PHYSICAL solid (reflecting) walls,
  !! for the curvilinear residual's weak pressure-wall flux.  Mirrors
  !! apply_halos_2d's serial/MPI dual-path edge selection exactly (single source
  !! of truth): serial (np=1 or no Cart comm) => all four edges are physical;
  !! MPI => an edge is physical only if its Cartesian neighbour is MPI_PROC_NULL.
  !! An edge is a wall iff it is physical AND its cfg % bc_* is 'reflecting'.
  subroutine wall_edges_2d(state, w_xlo, w_xhi, w_ylo, w_yhi)
    use mpi_cart_2d, only: cart_is_setup
    use mpi_runtime, only: mpi_proc_null_value, n_ranks
    type(solver_state_2d_t), intent(in) :: state
    logical, intent(out) :: w_xlo, w_xhi, w_ylo, w_yhi
    logical :: phys_xlo, phys_xhi, phys_ylo, phys_yhi
    integer :: null_rank

    if (n_ranks() == 1 .or. .not. cart_is_setup()) then
      phys_xlo = .true.; phys_xhi = .true.; phys_ylo = .true.; phys_yhi = .true.
    else
      null_rank = mpi_proc_null_value()
      phys_xlo = (state % decomp_2d % x_lo_neighbour == null_rank)
      phys_xhi = (state % decomp_2d % x_hi_neighbour == null_rank)
      phys_ylo = (state % decomp_2d % y_lo_neighbour == null_rank)
      phys_yhi = (state % decomp_2d % y_hi_neighbour == null_rank)
    end if
    w_xlo = phys_xlo .and. trim(state % cfg % bc_left) == bc_reflecting
    w_xhi = phys_xhi .and. trim(state % cfg % bc_right) == bc_reflecting
    w_ylo = phys_ylo .and. trim(state % cfg % bc_bottom) == bc_reflecting
    w_yhi = phys_yhi .and. trim(state % cfg % bc_top) == bc_reflecting
  end subroutine wall_edges_2d

  !> Fill all four halo bands by periodic wrap of the opposite interior edge.
  subroutine apply_periodic_halos_2d(state)
    type(solver_state_2d_t), intent(inout) :: state
    integer :: i, j, k, h, nx, ny
    h = state % halo_width; nx = state % nx_local; ny = state % ny_local

    ! x halos (wrap columns), full y interior
    do j = 1, ny
      do k = 1, h
        state % ub(:, 1 - k, j) = state % ub(:, nx + 1 - k, j)   ! left  <- right interior
        state % ub(:, nx + k, j) = state % ub(:, k, j)           ! right <- left interior
      end do
    end do
    ! y halos (wrap rows), full x range INCLUDING the x-halos just filled (corner consistency)
    do k = 1, h
      do i = 1 - h, nx + h
        state % ub(:, i, 1 - k) = state % ub(:, i, ny + 1 - k)   ! bottom <- top interior
        state % ub(:, i, ny + k) = state % ub(:, i, k)           ! top    <- bottom interior
      end do
    end do
  end subroutine apply_periodic_halos_2d

  !> Fill all four halo bands per cfg % bc_* edge type.
  !!
  !! Dual-path: if no MPI Cartesian communicator has been set up (serial unit
  !! tests), falls back to the Phase-2C pure per-edge fills for all four edges.
  !! When the cart comm exists (MPI app / multi-rank tests), calls
  !! exchange_halos_2d to fill from real neighbours and then applies the
  !! physical BC fill only on edges whose Cart neighbour is MPI_PROC_NULL.
  !!
  !! Diagonal corner ghosts are intentionally left unfilled (the
  !! dimension-by-dimension residual never reads them).
  subroutine apply_halos_2d(state)
    use halo_exchange_2d, only: exchange_halos_2d
    use mpi_cart_2d, only: cart_is_setup
    use mpi_runtime, only: mpi_proc_null_value, n_ranks
    type(solver_state_2d_t), intent(inout) :: state
    integer :: null_rank

    ! NODAL finite-difference path: route to the node-anchored BC kernel.  The
    ! FDM grid reflects/wraps about the boundary NODE (not the cell face), so the
    ! cell-centred FVM fills below would be off by half a cell.  Keep the FVM body
    ! byte-identical for the default (FVM) path.
    if (trim(state % blocks(1) % method) == method_fdm) then
      call apply_bcs_fdm_2d(state)
      return
    end if

    if (n_ranks() == 1 .or. .not. cart_is_setup()) then
      ! Serial (single rank, or no MPI cart): pure per-edge physical/periodic fills.
      call fill_x_edge(state, .true., trim(state % cfg % bc_left))
      call fill_x_edge(state, .false., trim(state % cfg % bc_right))
      call fill_y_edge(state, .true., trim(state % cfg % bc_bottom))
      call fill_y_edge(state, .false., trim(state % cfg % bc_top))
      return
    end if

    ! MPI mode: exchange with real neighbours, then physical fills at PROC_NULL edges.
    null_rank = mpi_proc_null_value()
    call exchange_halos_2d(state, state % decomp_2d)
    if (state % decomp_2d % x_lo_neighbour == null_rank) &
      call fill_x_edge(state, .true., trim(state % cfg % bc_left))
    if (state % decomp_2d % x_hi_neighbour == null_rank) &
      call fill_x_edge(state, .false., trim(state % cfg % bc_right))
    if (state % decomp_2d % y_lo_neighbour == null_rank) &
      call fill_y_edge(state, .true., trim(state % cfg % bc_bottom))
    if (state % decomp_2d % y_hi_neighbour == null_rank) &
      call fill_y_edge(state, .false., trim(state % cfg % bc_top))
  end subroutine apply_halos_2d

  !> Fill the four halo bands for the NODAL finite-difference (FDM) path.
  !!
  !! Node anchoring: the physical boundary lies ON the boundary node (column 1 /
  !! nx_local, row 1 / ny_local), not on a cell face, so the ghost fills differ
  !! from the cell-centred FVM kernel:
  !!   * reflecting: mirror about the boundary NODE — ghost(1-k)=node(1+k) with the
  !!     wall-normal momentum negated (vs. the FVM face mirror ghost(1-k)=cell(k)).
  !!   * periodic: duplicate-endpoint-aware nodal wrap.  The nodal mesh has nodes
  !!     ON both boundaries with node 1 == node nx_local the SAME physical point
  !!     (period = (nx_local-1)*dx), so the wrap SKIPS the duplicate endpoint:
  !!     ub(1-k)=ub(nx_local-k), ub(nx_local+k)=ub(1+k).  This is the 2D analogue
  !!     of the verified-correct 1D convention (boundary_conditions::apply_bcs_fdm,
  !!     n_pt = n_cell+1 with ub(:,1)==ub(:,n_pt)).  It differs from the FVM
  !!     cell-centred wrap at the boundary (FVM cell averages have no duplicate
  !!     endpoint), so the FDM operator equals FVM only in the INTERIOR away from
  !!     the boundary ghosts — which is what the equivalence oracle in
  !!     test_spatial_discretization_fdm_2d now pins.
  !!   * outflow: zeroth-order extrapolation (copy the boundary node).
  !!   * supersonic_inlet: free-stream conserved state.
  !!
  !! Dual-path identical to apply_halos_2d: serial (np=1 / no cart comm) fills all
  !! four physical edges; under MPI the interior partition edges are filled by
  !! exchange_halos_2d and only the global-boundary edges (Cart neighbour ==
  !! MPI_PROC_NULL) get a physical-BC fill, so no interior edge is double-filled.
  !!
  !! NOTE (np>1 periodic): the MPI cartesian wrap in exchange_halos_2d performs
  !! the neighbour copy across the periodic boundary; verifying it stays
  !! consistent with the nodal grid under decomposition is deferred to the
  !! Phase 5 np>1 full-solve task.
  subroutine apply_bcs_fdm_2d(state)
    use halo_exchange_2d, only: exchange_halos_2d
    use mpi_cart_2d, only: cart_is_setup
    use mpi_runtime, only: mpi_proc_null_value, n_ranks
    type(solver_state_2d_t), intent(inout) :: state
    integer :: null_rank

    if (n_ranks() == 1 .or. .not. cart_is_setup()) then
      call fill_x_edge_fdm(state, .true., trim(state % cfg % bc_left))
      call fill_x_edge_fdm(state, .false., trim(state % cfg % bc_right))
      call fill_y_edge_fdm(state, .true., trim(state % cfg % bc_bottom))
      call fill_y_edge_fdm(state, .false., trim(state % cfg % bc_top))
      return
    end if

    null_rank = mpi_proc_null_value()
    call exchange_halos_2d(state, state % decomp_2d)
    if (state % decomp_2d % x_lo_neighbour == null_rank) &
      call fill_x_edge_fdm(state, .true., trim(state % cfg % bc_left))
    if (state % decomp_2d % x_hi_neighbour == null_rank) &
      call fill_x_edge_fdm(state, .false., trim(state % cfg % bc_right))
    if (state % decomp_2d % y_lo_neighbour == null_rank) &
      call fill_y_edge_fdm(state, .true., trim(state % cfg % bc_bottom))
    if (state % decomp_2d % y_hi_neighbour == null_rank) &
      call fill_y_edge_fdm(state, .false., trim(state % cfg % bc_top))
  end subroutine apply_bcs_fdm_2d

  !> Fill a left (is_min=.true.) or right x-edge band over rows j=1..ny, NODAL.
  subroutine fill_x_edge_fdm(state, is_min, bc)
    type(solver_state_2d_t), intent(inout) :: state
    logical, intent(in) :: is_min
    character(len=*), intent(in) :: bc
    integer :: j, k, h, nx, gc, src
    real(wp) :: q_inf(size(state % ub, 1)), q_ff_ref(size(state % ub, 1))
    h = state % halo_width; nx = state % nx_local
    if (bc == bc_supersonic_inlet) q_inf = ref_state(state)
    if (bc == bc_farfield) q_ff_ref = ref_state(state)
    do j = 1, state % ny_local
      do k = 1, h
        if (is_min) then; gc = 1 - k; else; gc = nx + k; end if
        select case (bc)
        case (bc_periodic)
          ! Duplicate-endpoint-aware nodal wrap (period (nx-1)*dx; node 1 == node
          ! nx are the SAME physical point).  Skip the duplicate endpoint so the
          ! periodic stencil has no repeated point: ub(1-k)=ub(nx-k),
          ! ub(nx+k)=ub(1+k).  This is the 2D analogue of the verified-correct 1D
          ! convention (boundary_conditions::apply_bcs_fdm).
          if (is_min) then; src = nx - k; else; src = 1 + k; end if
          state % ub(:, gc, j) = state % ub(:, src, j)
        case (bc_outflow)
          if (is_min) then; src = 1; else; src = nx; end if
          state % ub(:, gc, j) = state % ub(:, src, j)
        case (bc_reflecting)
          ! Node mirror about the boundary node; wall-normal momentum negated.
          if (is_min) then; src = 1 + k; else; src = nx - k; end if
          if (state % mesh % fdm_curvilinear) then
            ! Curvilinear nodal slip-wall: reflect the interior partner's
            ! velocity about the boundary-NODE metric normal sx_xi(:,fcol,j)
            ! (S_xi = (y_eta,-x_eta)).  Every ghost layer uses the single
            ! boundary-node normal, matching the FVM curvilinear branch in
            ! fill_x_edge (node-anchored here rather than face-anchored).
            block
              real(wp) :: nrm(2), smag
              integer :: fcol
              if (is_min) then; fcol = 1; else; fcol = nx + 1; end if
              smag = sqrt(state % mesh % sx_xi(1, fcol, j)**2 + state % mesh % sx_xi(2, fcol, j)**2)
              if (smag <= 0.0_wp) error stop 'boundary_2d: degenerate (zero-length) x-boundary node metric'
              nrm = state % mesh % sx_xi(:, fcol, j) / smag
              state % ub(:, gc, j) = reflect_velocity_2d(state % ub(:, src, j), nrm)
            end block
          else
            state % ub(:, gc, j) = state % ub(:, src, j)
            state % ub(2, gc, j) = -state % ub(2, src, j)
          end if
        case (bc_supersonic_inlet)
          state % ub(:, gc, j) = q_inf
        case (bc_farfield)
          ! Characteristic far-field about the OUTWARD boundary-NODE metric
          ! normal (sx_xi points +xi; left-edge outward normal is -sx_xi).
          ! Uniform nodal FDM has no nodal metric -> use the axis normal.
          block
            real(wp) :: nrm(2), smag
            integer :: fcol, isrc
            if (is_min) then; fcol = 1; isrc = 1; else; fcol = nx + 1; isrc = nx; end if
            if (state % mesh % fdm_curvilinear) then
              smag = sqrt(state % mesh % sx_xi(1, fcol, j)**2 + state % mesh % sx_xi(2, fcol, j)**2)
              if (smag <= 0.0_wp) error stop 'boundary_2d: degenerate x-boundary node metric (farfield)'
              nrm = state % mesh % sx_xi(:, fcol, j) / smag
              if (is_min) nrm = -nrm
            else
              nrm = 0.0_wp
              if (is_min) then; nrm(1) = -1.0_wp; else; nrm(1) = 1.0_wp; end if
            end if
            state % ub(:, gc, j) = set_farfield_ghost(state % ub(:, isrc, j), q_ff_ref, nrm, state % cfg % gam)
          end block
        case default
          error stop 'boundary_2d: unknown x-edge FDM BC "'//bc//'"'
        end select
      end do
    end do
  end subroutine fill_x_edge_fdm

  !> Fill a bottom (is_min=.true.) or top y-edge band over columns i=1..nx, NODAL.
  subroutine fill_y_edge_fdm(state, is_min, bc)
    type(solver_state_2d_t), intent(inout) :: state
    logical, intent(in) :: is_min
    character(len=*), intent(in) :: bc
    integer :: i, k, h, ny, gc, src
    real(wp) :: q_inf(size(state % ub, 1)), q_ff_ref(size(state % ub, 1))
    h = state % halo_width; ny = state % ny_local
    if (bc == bc_supersonic_inlet) q_inf = ref_state(state)
    if (bc == bc_farfield) q_ff_ref = ref_state(state)
    do k = 1, h
      do i = 1, state % nx_local
        if (is_min) then; gc = 1 - k; else; gc = ny + k; end if
        select case (bc)
        case (bc_periodic)
          ! Duplicate-endpoint-aware nodal wrap (period (ny-1)*dy; row 1 == row ny
          ! are the SAME physical point).  Skip the duplicate endpoint:
          ! ub(1-k)=ub(ny-k), ub(ny+k)=ub(1+k).  Matches the 1D convention.
          if (is_min) then; src = ny - k; else; src = 1 + k; end if
          state % ub(:, i, gc) = state % ub(:, i, src)
        case (bc_outflow)
          if (is_min) then; src = 1; else; src = ny; end if
          state % ub(:, i, gc) = state % ub(:, i, src)
        case (bc_reflecting)
          ! Node mirror about the boundary node; wall-normal momentum negated.
          if (is_min) then; src = 1 + k; else; src = ny - k; end if
          if (state % mesh % fdm_curvilinear) then
            ! Curvilinear nodal slip-wall: reflect about the boundary-NODE metric
            ! normal sx_eta(:,i,frow) (S_eta = (-y_xi,x_xi)); see fill_x_edge_fdm.
            block
              real(wp) :: nrm(2), smag
              integer :: frow
              if (is_min) then; frow = 1; else; frow = ny + 1; end if
              smag = sqrt(state % mesh % sx_eta(1, i, frow)**2 + state % mesh % sx_eta(2, i, frow)**2)
              if (smag <= 0.0_wp) error stop 'boundary_2d: degenerate (zero-length) y-boundary node metric'
              nrm = state % mesh % sx_eta(:, i, frow) / smag
              state % ub(:, i, gc) = reflect_velocity_2d(state % ub(:, i, src), nrm)
            end block
          else
            state % ub(:, i, gc) = state % ub(:, i, src)
            state % ub(3, i, gc) = -state % ub(3, i, src)
          end if
        case (bc_supersonic_inlet)
          state % ub(:, i, gc) = q_inf
        case (bc_farfield)
          ! Characteristic far-field about the OUTWARD boundary-NODE metric
          ! normal (sx_eta points +eta; bottom-edge outward normal is -sx_eta).
          block
            real(wp) :: nrm(2), smag
            integer :: frow, jsrc
            if (is_min) then; frow = 1; jsrc = 1; else; frow = ny + 1; jsrc = ny; end if
            if (state % mesh % fdm_curvilinear) then
              smag = sqrt(state % mesh % sx_eta(1, i, frow)**2 + state % mesh % sx_eta(2, i, frow)**2)
              if (smag <= 0.0_wp) error stop 'boundary_2d: degenerate y-boundary node metric (farfield)'
              nrm = state % mesh % sx_eta(:, i, frow) / smag
              if (is_min) nrm = -nrm
            else
              nrm = 0.0_wp
              if (is_min) then; nrm(2) = -1.0_wp; else; nrm(2) = 1.0_wp; end if
            end if
            state % ub(:, i, gc) = set_farfield_ghost(state % ub(:, i, jsrc), q_ff_ref, nrm, state % cfg % gam)
          end block
        case default
          error stop 'boundary_2d: unknown y-edge FDM BC "'//bc//'"'
        end select
      end do
    end do
  end subroutine fill_y_edge_fdm

  !> Reflect the velocity of a conserved state about unit normal nrm, preserving
  !! density and total energy (|v| is unchanged). Inviscid no-penetration wall.
  !! The density is floored to eps (matching limit_pos_2d): the halo fill runs
  !! BEFORE any residual-side positivity guard can fire, so a transient
  !! non-positive boundary-adjacent state (near-vacuum off a strong expansion)
  !! must not divide-by-zero / SIGFPE here (audit 2026-07-06 N2); the residual
  !! guards still abort cleanly on the interior state itself.
  pure function reflect_velocity_2d(q, nrm) result(qr)
    real(wp), intent(in) :: q(:), nrm(2)
    real(wp) :: qr(size(q))
    real(wp), parameter :: eps = 1.0e-13_wp
    real(wp) :: rho, u, v, vn
    rho = max(q(1), eps)
    u = q(2) / rho
    v = q(3) / rho
    vn = u * nrm(1) + v * nrm(2)
    qr(1) = rho
    qr(2) = rho * (u - 2.0_wp * vn * nrm(1))
    qr(3) = rho * (v - 2.0_wp * vn * nrm(2))
    qr(4) = q(4)
  end function reflect_velocity_2d

  !> Conserved free-stream reference state from config.
  pure function ref_state(state) result(q)
    type(solver_state_2d_t), intent(in) :: state
    real(wp) :: q(size(state % ub, 1))
    associate (c => state % cfg)
      q(1) = c % rho_inf
      q(2) = c % rho_inf * c % u_inf
      q(3) = c % rho_inf * c % v_inf
      q(4) = c % p_inf / (c % gam - 1.0_wp) + 0.5_wp * c % rho_inf * (c % u_inf**2 + c % v_inf**2)
    end associate
  end function ref_state

  !> Characteristic (Riemann-invariant) far-field ghost state for an OUTWARD unit
  !! normal nrm.  q_int is the boundary-adjacent interior conserved state, q_inf
  !! the free-stream.  With interior/free-stream normal velocities Un_i, Un_inf
  !! and sound speeds c_i, c_inf:
  !!   * supersonic in/out (|Un_i| >= c_i) short-circuit to pure free-stream /
  !!     pure interior extrapolation;
  !!   * otherwise combine the outgoing (interior) invariant
  !!     R+ = Un_i + 2 c_i/(gam-1) with the incoming (free-stream)
  !!     R- = Un_inf - 2 c_inf/(gam-1) to get Un_b = (R+ + R-)/2,
  !!     c_b = (gam-1)/4 (R+ - R-).  Entropy and tangential velocity are taken
  !!     from the free-stream on inflow (Un_b <= 0) and from the interior on
  !!     outflow; rho_b, p_b follow from entropy + c_b.
  !! gam is passed explicitly (no module-level thermodynamic state here).
  pure function set_farfield_ghost(q_int, q_inf, nrm, gam) result(q_ghost)
    real(wp), intent(in) :: q_int(:), q_inf(:), nrm(2), gam
    real(wp) :: q_ghost(size(q_int))
    real(wp), parameter :: eps = 1.0e-13_wp
    real(wp) :: qi(size(q_int))
    real(wp) :: gm1, rho_i, u_i, v_i, p_i, c_i, un_i
    real(wp) :: rho_o, u_o, v_o, p_o, c_o, un_o
    real(wp) :: rp, rm, un_b, c_b, s_b, ut, vt, rho_b, p_b, ux, uy
    gm1 = gam - 1.0_wp
    ! Positivity guard (audit 2026-07-06 N2): the halo fill runs before the
    ! residual guards, so a transient non-positive interior density/pressure
    ! must not reach the divide / sqrt below (NaN ghosts, or SIGFPE under
    ! -ffpe-trap, with no clean abort). Floor to eps like limit_pos_2d; with
    ! rho_i, p_i > 0 every downstream expression (c_i, s_b, rho_b) stays
    ! finite. The free-stream reference is config-validated positive.
    qi = q_int
    qi(1) = max(qi(1), eps)
    call primitives_2d(qi, gam, rho_i, u_i, v_i, p_i)
    p_i = max(p_i, eps)
    call primitives_2d(q_inf, gam, rho_o, u_o, v_o, p_o)
    c_i = sqrt(gam * p_i / rho_i)
    c_o = sqrt(gam * p_o / rho_o)
    un_i = u_i * nrm(1) + v_i * nrm(2)
    un_o = u_o * nrm(1) + v_o * nrm(2)
    if (un_i <= -c_i) then
      q_ghost = q_inf; return          ! supersonic inflow -> pure free-stream
    else if (un_i >= c_i) then
      q_ghost = q_int; return          ! supersonic outflow -> pure extrapolation
    end if
    rp = un_i + 2.0_wp * c_i / gm1      ! outgoing invariant (from interior)
    rm = un_o - 2.0_wp * c_o / gm1      ! incoming invariant (from free-stream)
    un_b = 0.5_wp * (rp + rm)
    c_b = 0.25_wp * gm1 * (rp - rm)
    if (un_b <= 0.0_wp) then            ! inflow: entropy + tangential from free-stream
      s_b = p_o / rho_o**gam
      ut = u_o - un_o * nrm(1)
      vt = v_o - un_o * nrm(2)
    else                               ! outflow: entropy + tangential from interior
      s_b = p_i / rho_i**gam
      ut = u_i - un_i * nrm(1)
      vt = v_i - un_i * nrm(2)
    end if
    rho_b = (c_b * c_b / (gam * s_b))**(1.0_wp / gm1)
    p_b = rho_b * c_b * c_b / gam
    ux = ut + un_b * nrm(1)
    uy = vt + un_b * nrm(2)
    q_ghost(1) = rho_b
    q_ghost(2) = rho_b * ux
    q_ghost(3) = rho_b * uy
    q_ghost(4) = p_b / gm1 + 0.5_wp * rho_b * (ux * ux + uy * uy)
  end function set_farfield_ghost

  !> Fill a left (is_min=.true.) or right x-edge band over rows j=1..ny.
  subroutine fill_x_edge(state, is_min, bc)
    type(solver_state_2d_t), intent(inout) :: state
    logical, intent(in) :: is_min
    character(len=*), intent(in) :: bc
    integer :: j, k, h, nx, gc, src
    real(wp) :: q_inf(size(state % ub, 1)), q_ff_ref(size(state % ub, 1))
    h = state % halo_width; nx = state % nx_local
    if (bc == bc_supersonic_inlet) q_inf = ref_state(state)
    if (bc == bc_farfield) q_ff_ref = ref_state(state)
    do j = 1, state % ny_local
      do k = 1, h
        if (is_min) then; gc = 1 - k; else; gc = nx + k; end if
        select case (bc)
        case (bc_periodic)
          if (is_min) then; src = nx + 1 - k; else; src = k; end if
          state % ub(:, gc, j) = state % ub(:, src, j)
        case (bc_outflow)
          if (is_min) then; src = 1; else; src = nx; end if
          state % ub(:, gc, j) = state % ub(:, src, j)
        case (bc_reflecting)
          if (is_min) then; src = k; else; src = nx + 1 - k; end if
          if (state % mesh % uniform) then
            state % ub(:, gc, j) = state % ub(:, src, j)
            state % ub(2, gc, j) = -state % ub(2, src, j)
          else
            ! Curvilinear reflecting wall: all halo layers (gc = 1-k) reflect
            ! their interior partner about the SINGLE boundary-face normal
            ! s_xi(:,fcol,j). The ghost layers march along the wall-normal
            ! direction (into the domain), not along the wall, so there is no
            ! distinct per-layer wall normal to use; the wall face is the
            ! correct reflection plane for every layer. No-penetration is
            ! enforced exactly at the boundary face. This single-normal
            ! treatment is standard practice for curvilinear FV walls. See
            ! doc/theory-guide.org (curvilinear boundary conditions).
            block
              real(wp) :: nrm(2), smag
              integer :: fcol
              if (is_min) then; fcol = 1; else; fcol = nx + 1; end if
              smag = sqrt(state % mesh % s_xi(1, fcol, j)**2 + state % mesh % s_xi(2, fcol, j)**2)
              if (smag <= 0.0_wp) error stop 'boundary_2d: degenerate (zero-length) x-boundary face'
              nrm = state % mesh % s_xi(:, fcol, j) / smag
              state % ub(:, gc, j) = reflect_velocity_2d(state % ub(:, src, j), nrm)
            end block
          end if
        case (bc_supersonic_inlet)
          state % ub(:, gc, j) = q_inf
        case (bc_farfield)
          ! Characteristic far-field about the OUTWARD face normal.  s_xi points
          ! +xi (into the domain at the left face, out at the right), so the
          ! left-edge outward normal is -s_xi.  Uniform FVM has no s_xi -> use
          ! the axis normal.
          block
            real(wp) :: nrm(2), smag
            integer :: fcol, isrc
            if (is_min) then; fcol = 1; isrc = 1; else; fcol = nx + 1; isrc = nx; end if
            if (state % mesh % uniform) then
              nrm = 0.0_wp
              if (is_min) then; nrm(1) = -1.0_wp; else; nrm(1) = 1.0_wp; end if
            else
              smag = sqrt(state % mesh % s_xi(1, fcol, j)**2 + state % mesh % s_xi(2, fcol, j)**2)
              if (smag <= 0.0_wp) error stop 'boundary_2d: degenerate x-boundary face (farfield)'
              nrm = state % mesh % s_xi(:, fcol, j) / smag
              if (is_min) nrm = -nrm
            end if
            state % ub(:, gc, j) = set_farfield_ghost(state % ub(:, isrc, j), q_ff_ref, nrm, state % cfg % gam)
          end block
        case default
          error stop 'boundary_2d: unknown x-edge BC "'//bc//'"'
        end select
      end do
    end do
  end subroutine fill_x_edge

  !> Fill a bottom (is_min=.true.) or top y-edge band over columns i=1..nx.
  subroutine fill_y_edge(state, is_min, bc)
    type(solver_state_2d_t), intent(inout) :: state
    logical, intent(in) :: is_min
    character(len=*), intent(in) :: bc
    integer :: i, k, h, ny, gc, src
    real(wp) :: q_inf(size(state % ub, 1)), q_ff_ref(size(state % ub, 1))
    h = state % halo_width; ny = state % ny_local
    if (bc == bc_supersonic_inlet) q_inf = ref_state(state)
    if (bc == bc_farfield) q_ff_ref = ref_state(state)
    do k = 1, h
      do i = 1, state % nx_local
        if (is_min) then; gc = 1 - k; else; gc = ny + k; end if
        select case (bc)
        case (bc_periodic)
          if (is_min) then; src = ny + 1 - k; else; src = k; end if
          state % ub(:, i, gc) = state % ub(:, i, src)
        case (bc_outflow)
          if (is_min) then; src = 1; else; src = ny; end if
          state % ub(:, i, gc) = state % ub(:, i, src)
        case (bc_reflecting)
          if (is_min) then; src = k; else; src = ny + 1 - k; end if
          if (state % mesh % uniform) then
            state % ub(:, i, gc) = state % ub(:, i, src)
            state % ub(3, i, gc) = -state % ub(3, i, src)
          else
            ! Curvilinear reflecting wall: see fill_x_edge — all halo layers
            ! reflect about the single boundary-face normal s_eta(:,i,frow);
            ! the wall face is the correct reflection plane for every layer.
            block
              real(wp) :: nrm(2), smag
              integer :: frow
              if (is_min) then; frow = 1; else; frow = ny + 1; end if
              smag = sqrt(state % mesh % s_eta(1, i, frow)**2 + state % mesh % s_eta(2, i, frow)**2)
              if (smag <= 0.0_wp) error stop 'boundary_2d: degenerate (zero-length) y-boundary face'
              nrm = state % mesh % s_eta(:, i, frow) / smag
              state % ub(:, i, gc) = reflect_velocity_2d(state % ub(:, i, src), nrm)
            end block
          end if
        case (bc_supersonic_inlet)
          state % ub(:, i, gc) = q_inf
        case (bc_farfield)
          ! Characteristic far-field about the OUTWARD face normal (s_eta points
          ! +eta; bottom-edge outward normal is -s_eta).  Uniform FVM -> axis.
          block
            real(wp) :: nrm(2), smag
            integer :: frow, jsrc
            if (is_min) then; frow = 1; jsrc = 1; else; frow = ny + 1; jsrc = ny; end if
            if (state % mesh % uniform) then
              nrm = 0.0_wp
              if (is_min) then; nrm(2) = -1.0_wp; else; nrm(2) = 1.0_wp; end if
            else
              smag = sqrt(state % mesh % s_eta(1, i, frow)**2 + state % mesh % s_eta(2, i, frow)**2)
              if (smag <= 0.0_wp) error stop 'boundary_2d: degenerate y-boundary face (farfield)'
              nrm = state % mesh % s_eta(:, i, frow) / smag
              if (is_min) nrm = -nrm
            end if
            state % ub(:, i, gc) = set_farfield_ghost(state % ub(:, i, jsrc), q_ff_ref, nrm, state % cfg % gam)
          end block
        case default
          error stop 'boundary_2d: unknown y-edge BC "'//bc//'"'
        end select
      end do
    end do
  end subroutine fill_y_edge

end module boundary_2d
