!> @file mesh_2d.f90
!> @brief 2D single-block structured (curvilinear) grid metric layer.
!!
!! Owns node coordinates and the derived finite-volume metrics for the 2D Euler
!! solver: cell areas (vol), face area-normal vectors (s_xi, s_eta), and cell
!! centres (xc, yc). Geometry is read from a single-block Plot3D grid file
!! (ASCII or binary). On a uniform Cartesian grid the metrics are unused — the
!! residual keeps its scalar dx/dy fast path (bit-for-bit). Phase A builds the
!! whole-domain (serial) metrics; local slicing for MPI is added in Phase D.
!!
!! Dependency rule: this module must NOT `use config_2d` (config_2d uses it).
module mesh_2d
  use precision, only: wp
  use, intrinsic :: iso_fortran_env, only: int64
  use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
  implicit none
  private

  public :: mesh_2d_t
  public :: read_plot3d_2d
  public :: check_plot3d_dims
  public :: build_mesh_2d_global
  public :: build_mesh_2d_uniform_nodal
  public :: compute_metrics_2d
  public :: compute_metrics_fdm_2d
  public :: build_mesh_2d_fdm_curvilinear
  public :: build_mesh_2d_local
  public :: MESH2D_MAX_DIM

  !> Grid metrics for the 2D structured scheme.
  type :: mesh_2d_t
    logical :: uniform = .true.            !< true => scalar dx/dy fast path
    integer :: nx = 0, ny = 0             !< cell count per axis (FVM) or node count (FDM nodal mesh)
    real(wp) :: dx_uniform = 0.0_wp, dy_uniform = 0.0_wp
    real(wp), allocatable :: x_node(:, :)  !< node coords; shape (nx+1,ny+1) for FVM; (mesh%nx,mesh%ny) for the FDM nodal mesh
    real(wp), allocatable :: y_node(:, :)
    real(wp), allocatable :: vol(:, :)     !< cell area V_{i,j} (nx, ny)
    real(wp), allocatable :: s_xi(:, :, :) !< xi-face vectors (2, nx+1, ny)
    real(wp), allocatable :: s_eta(:, :, :)!< eta-face vectors (2, nx, ny+1)
    real(wp), allocatable :: xc(:, :)      !< cell-centre x (nx, ny)
    real(wp), allocatable :: yc(:, :)      !< cell-centre y (nx, ny)
    !> Nodal FDM transformation metrics (allocated only for curvilinear FDM).
    real(wp), allocatable :: sx_xi(:, :, :)   !< S_xi=(y_eta,-x_eta) at xi-faces (2, nx+1, ny)
    real(wp), allocatable :: sx_eta(:, :, :)  !< S_eta=(-y_xi,x_xi) at eta-faces (2, nx, ny+1)
    real(wp), allocatable :: jac(:, :)        !< Jacobian J>0 at nodes (nx, ny)
    logical :: fdm_curvilinear = .false.      !< true when nodal metrics are present
  end type mesh_2d_t

  !> Sanity caps on Plot3D header dimensions. Per-axis node cap mirrors the
  !! spirit of the 1D n_cell schema bound; the total-node cap is checked via an
  !! int64 product so a pathological header cannot overflow the default-integer
  !! ni*nj before the allocation's stat= is consulted.
  integer, parameter :: MESH2D_MAX_DIM = 1000000
  integer(int64), parameter :: MESH2D_MAX_POINTS = 100000000_int64

contains

  !> Read a single-block 2D Plot3D grid. Header `ni nj`, then all x(i,j)
  !! (i fastest), then all y(i,j). ASCII (binary=.false.) or Fortran
  !! unformatted sequential (binary=.true.).
  subroutine read_plot3d_2d(path, binary, xg, yg, ni, nj, is_ok, message)
    character(len=*), intent(in) :: path
    logical, intent(in) :: binary
    real(wp), allocatable, intent(out) :: xg(:, :), yg(:, :)
    integer, intent(out) :: ni, nj
    logical, intent(out) :: is_ok
    character(len=*), intent(out) :: message

    integer :: u, info, astat, i, j

    is_ok = .false.
    message = ''
    ni = 0; nj = 0

    if (binary) then
      open (newunit=u, file=trim(path), status='old', action='read', &
            form='unformatted', access='sequential', iostat=info)
      if (info /= 0) then
        message = 'mesh_2d: cannot open binary grid_file "'//trim(path)//'"'
        return
      end if
      read (u, iostat=info) ni, nj
      if (info /= 0) then
        close (u, iostat=info)
        message = 'mesh_2d: malformed Plot3D binary header'
        return
      end if
      block
        logical :: dims_ok
        call check_plot3d_dims(ni, nj, dims_ok, message)
        if (.not. dims_ok) then
          close (u, iostat=info)
          return
        end if
      end block
      allocate (xg(ni, nj), yg(ni, nj), stat=astat)
      if (astat /= 0) then
        close (u, iostat=info)
        message = 'mesh_2d: allocation failed for grid arrays'
        return
      end if
      read (u, iostat=info) ((xg(i, j), i=1, ni), j=1, nj)
      if (info == 0) read (u, iostat=info) ((yg(i, j), i=1, ni), j=1, nj)
      close (u, iostat=astat)
      if (info /= 0) then
        deallocate (xg, yg, stat=astat)
        message = 'mesh_2d: malformed/short Plot3D binary coordinate block'// &
                  ' (also check host endianness for binary grids)'
        return
      end if
      if (.not. (all_finite_2d(xg) .and. all_finite_2d(yg))) then
        deallocate (xg, yg, stat=astat)
        message = 'mesh_2d: non-finite (NaN/Inf) coordinate in binary Plot3D grid'// &
                  ' (truncated record or wrong endianness?)'
        return
      end if
      is_ok = .true.
      return
    end if

    ! ASCII path.
    open (newunit=u, file=trim(path), status='old', action='read', iostat=info)
    if (info /= 0) then
      message = 'mesh_2d: cannot open grid_file "'//trim(path)//'"'
      return
    end if
    read (u, *, iostat=info) ni, nj
    if (info /= 0) then
      close (u, iostat=info)
      message = 'mesh_2d: malformed Plot3D header (expected "ni nj")'
      return
    end if
    block
      logical :: dims_ok
      call check_plot3d_dims(ni, nj, dims_ok, message)
      if (.not. dims_ok) then
        close (u, iostat=info)
        return
      end if
    end block
    allocate (xg(ni, nj), yg(ni, nj), stat=astat)
    if (astat /= 0) then
      close (u, iostat=info)
      message = 'mesh_2d: allocation failed for grid arrays'
      return
    end if
    read (u, *, iostat=info) ((xg(i, j), i=1, ni), j=1, nj)
    if (info == 0) read (u, *, iostat=info) ((yg(i, j), i=1, ni), j=1, nj)
    close (u, iostat=astat)
    if (info /= 0) then
      deallocate (xg, yg, stat=astat)
      message = 'mesh_2d: malformed/short Plot3D coordinate block'
      return
    end if
    if (.not. (all_finite_2d(xg) .and. all_finite_2d(yg))) then
      deallocate (xg, yg, stat=astat)
      message = 'mesh_2d: non-finite (NaN/Inf) coordinate in Plot3D grid'
      return
    end if
    is_ok = .true.
  end subroutine read_plot3d_2d

  !> True iff every element of a 2D coordinate array is finite (no NaN/Inf).
  pure logical function all_finite_2d(a) result(ok)
    real(wp), intent(in) :: a(:, :)
    ok = all(ieee_is_finite(a))
  end function all_finite_2d

  !> Validate a Plot3D header: lower bound (>= 2 nodes per axis) and upper
  !! sanity caps (per-axis and total node count). Returns is_ok=.false. with a
  !! descriptive message on violation.
  subroutine check_plot3d_dims(ni, nj, is_ok, message)
    integer, intent(in) :: ni, nj
    logical, intent(out) :: is_ok
    character(len=*), intent(out) :: message
    is_ok = .false.
    message = ''
    ! Per-axis and total-node caps are checked in separate branches so the
    ! diagnostic names the cap that was actually exceeded. (The int64 product
    ! cannot overflow: both operands are promoted to int64 before multiplying,
    ! and any int32 ni*nj is <= ~4.6e18 < huge(int64).)
    if (ni < 2 .or. nj < 2) then
      write (message, '(A,I0,A,I0)') 'mesh_2d: grid too small ni=', ni, ' nj=', nj
    else if (ni > MESH2D_MAX_DIM .or. nj > MESH2D_MAX_DIM) then
      write (message, '(A,I0,A,I0,A,I0,A)') 'mesh_2d: grid dimension exceeds per-axis cap ni=', &
        ni, ' nj=', nj, ' (max ', MESH2D_MAX_DIM, ' per axis)'
    else if (int(ni, int64) * int(nj, int64) > MESH2D_MAX_POINTS) then
      write (message, '(A,I0,A,I0,A,I0,A)') 'mesh_2d: grid total node count exceeds cap ni=', &
        ni, ' nj=', nj, ' (max ', MESH2D_MAX_POINTS, ' total)'
    else
      is_ok = .true.
    end if
  end subroutine check_plot3d_dims

  !> Build the global (whole-domain) mesh from either a uniform specification or a
  !! Plot3D grid file.  For 'uniform', nx/ny are read from the caller.  For
  !! 'plot3d', ni/nj are read from the file and nx/ny are written back
  !! (intent(inout)) so the caller knows the derived sizes.  Calls
  !! compute_metrics_2d at the end and propagates its status — a metric failure
  !! (e.g. a folded or zero-area cell) fails the build.
  subroutine build_mesh_2d_global(mesh, grid_type, grid_file, grid_binary, &
                                  nx, ny, x_left, x_right, y_left, y_right, &
                                  is_ok, message)
    type(mesh_2d_t), intent(inout) :: mesh
    character(len=*), intent(in) :: grid_type, grid_file
    logical, intent(in) :: grid_binary
    integer, intent(inout) :: nx, ny
    real(wp), intent(in) :: x_left, x_right, y_left, y_right
    logical, intent(out) :: is_ok
    character(len=*), intent(out) :: message

    integer :: ni, nj, i, j, astat
    real(wp), allocatable :: xg(:, :), yg(:, :)

    is_ok = .true.
    message = ''

    if (trim(grid_type) == 'plot3d') then
      call read_plot3d_2d(grid_file, grid_binary, xg, yg, ni, nj, is_ok, message)
      if (.not. is_ok) return
      mesh % uniform = .false.
      mesh % nx = ni - 1
      mesh % ny = nj - 1
      mesh % dx_uniform = 0.0_wp
      mesh % dy_uniform = 0.0_wp
      call move_alloc(xg, mesh % x_node)
      call move_alloc(yg, mesh % y_node)
      nx = mesh % nx     ! file is authoritative
      ny = mesh % ny
    else if (trim(grid_type) == 'uniform') then
      if (nx <= 0 .or. ny <= 0) then
        is_ok = .false.
        message = 'mesh_2d: uniform grid requires nx > 0 and ny > 0'
        return
      end if
      mesh % uniform = .true.
      mesh % nx = nx
      mesh % ny = ny
      mesh % dx_uniform = (x_right - x_left) / real(nx, wp)
      mesh % dy_uniform = (y_right - y_left) / real(ny, wp)
      ni = nx + 1; nj = ny + 1
      if (allocated(mesh % x_node)) then
        deallocate (mesh % x_node, stat=astat)
        if (astat /= 0) error stop 'mesh_2d: x_node deallocate failed'
      end if
      if (allocated(mesh % y_node)) then
        deallocate (mesh % y_node, stat=astat)
        if (astat /= 0) error stop 'mesh_2d: y_node deallocate failed'
      end if
      allocate (mesh % x_node(ni, nj), mesh % y_node(ni, nj), stat=astat)
      if (astat /= 0) then
        is_ok = .false.
        message = 'mesh_2d: allocation failed for uniform node arrays'
        return
      end if
      do j = 1, nj
        do i = 1, ni
          mesh % x_node(i, j) = x_left + mesh % dx_uniform * real(i - 1, wp)
          mesh % y_node(i, j) = y_left + mesh % dy_uniform * real(j - 1, wp)
        end do
      end do
    else
      is_ok = .false.
      message = "mesh_2d: unknown grid_type '"//trim(grid_type)// &
                "' (expected 'uniform' or 'plot3d')"
      return
    end if

    call compute_metrics_2d(mesh, is_ok, message)
  end subroutine build_mesh_2d_global

  !> Compute finite-volume metrics from node coordinates already in mesh%x_node /
  !! mesh%y_node.  Allocates (with stat=) and fills vol(nx,ny), s_xi(2,nx+1,ny),
  !! s_eta(2,nx,ny+1), xc(nx,ny), yc(nx,ny).  Returns is_ok=.false. and a
  !! diagnostic message naming the first (i,j) with non-positive cell area.
  !!
  !! Metric definitions
  !! ------------------
  !! * xi-face vector at column iface, row j: edge (iface,j)->(iface,j+1).
  !!   s_xi = (dy, -dx), i.e. (+dy, -dx) where dx/dy are the edge increments.
  !! * eta-face vector at row jface, column i: edge (i,jface)->(i+1,jface).
  !!   s_eta = (-dy, dx).
  !! * Cell area via diagonal cross-product (exact for arbitrary quads):
  !!   vol = 0.5 * [(d1x*d2y) - (d2x*d1y)] where d1,d2 are the two diagonals.
  !! * Cell centre = arithmetic mean of four corners.
  subroutine compute_metrics_2d(mesh, is_ok, message)
    type(mesh_2d_t), intent(inout) :: mesh
    logical, intent(out) :: is_ok
    character(len=*), intent(out) :: message

    integer :: nx, ny, i, j, astat
    real(wp) :: dx, dy

    is_ok = .true.
    message = ''
    nx = mesh % nx; ny = mesh % ny

    if (allocated(mesh % vol)) then
      deallocate (mesh % vol, stat=astat)
      if (astat /= 0) error stop 'mesh_2d: vol deallocate failed'
    end if
    if (allocated(mesh % s_xi)) then
      deallocate (mesh % s_xi, stat=astat)
      if (astat /= 0) error stop 'mesh_2d: s_xi deallocate failed'
    end if
    if (allocated(mesh % s_eta)) then
      deallocate (mesh % s_eta, stat=astat)
      if (astat /= 0) error stop 'mesh_2d: s_eta deallocate failed'
    end if
    if (allocated(mesh % xc)) then
      deallocate (mesh % xc, stat=astat)
      if (astat /= 0) error stop 'mesh_2d: xc deallocate failed'
    end if
    if (allocated(mesh % yc)) then
      deallocate (mesh % yc, stat=astat)
      if (astat /= 0) error stop 'mesh_2d: yc deallocate failed'
    end if
    allocate (mesh % vol(nx, ny), mesh % xc(nx, ny), mesh % yc(nx, ny), &
              mesh % s_xi(2, nx + 1, ny), mesh % s_eta(2, nx, ny + 1), stat=astat)
    if (astat /= 0) then
      is_ok = .false.
      message = 'mesh_2d: allocation failed for metric arrays'
      return
    end if

    ! xi-face vectors at node columns iface=1..nx+1, edge (iface,j)->(iface,j+1).
    do j = 1, ny
      do i = 1, nx + 1
        dx = mesh % x_node(i, j + 1) - mesh % x_node(i, j)
        dy = mesh % y_node(i, j + 1) - mesh % y_node(i, j)
        mesh % s_xi(1, i, j) = dy
        mesh % s_xi(2, i, j) = -dx
        if (mesh % s_xi(1, i, j)**2 + mesh % s_xi(2, i, j)**2 <= 0.0_wp) then
          is_ok = .false.
          write (message, '(A,I0,A,I0,A)') 'mesh_2d: zero-length xi-face at (', i, ',', j, ')'
          return
        end if
      end do
    end do

    ! eta-face vectors at node rows jface=1..ny+1, edge (i,jface)->(i+1,jface).
    do j = 1, ny + 1
      do i = 1, nx
        dx = mesh % x_node(i + 1, j) - mesh % x_node(i, j)
        dy = mesh % y_node(i + 1, j) - mesh % y_node(i, j)
        mesh % s_eta(1, i, j) = -dy
        mesh % s_eta(2, i, j) = dx
        if (mesh % s_eta(1, i, j)**2 + mesh % s_eta(2, i, j)**2 <= 0.0_wp) then
          is_ok = .false.
          write (message, '(A,I0,A,I0,A)') 'mesh_2d: zero-length eta-face at (', i, ',', j, ')'
          return
        end if
      end do
    end do

    ! Cell areas (diagonal cross-product) and centres (corner mean).
    do j = 1, ny
      do i = 1, nx
        mesh % vol(i, j) = 0.5_wp * ( &
                           (mesh % x_node(i + 1, j + 1) - mesh % x_node(i, j)) * &
                           (mesh % y_node(i, j + 1) - mesh % y_node(i + 1, j)) - &
                           (mesh % x_node(i, j + 1) - mesh % x_node(i + 1, j)) * &
                           (mesh % y_node(i + 1, j + 1) - mesh % y_node(i, j)))
        if (mesh % vol(i, j) <= 0.0_wp) then
          is_ok = .false.
          write (message, '(A,I0,A,I0,A)') 'mesh_2d: non-positive cell area at (', &
            i, ',', j, ') — check node ordering / folded grid'
          return
        end if
        mesh % xc(i, j) = 0.25_wp * (mesh % x_node(i, j) + mesh % x_node(i + 1, j) + &
                                     mesh % x_node(i + 1, j + 1) + mesh % x_node(i, j + 1))
        mesh % yc(i, j) = 0.25_wp * (mesh % y_node(i, j) + mesh % y_node(i + 1, j) + &
                                     mesh % y_node(i + 1, j + 1) + mesh % y_node(i, j + 1))
      end do
    end do
  end subroutine compute_metrics_2d

  !> Stencil width for the high-order conservative-flux metric operator on n
  !! nodes: 6 (6th-order) when n>=6, 4 (4th-order) when n>=4, else 2 (the
  !! classical central scheme).  Even widths keep the interior operator a
  !! symmetric central difference, whose conservative flux telescopes cleanly.
  pure integer function cons_width(n) result(w)
    integer, intent(in) :: n
    if (n >= 6) then
      w = 6
    else if (n >= 4) then
      w = 4
    else
      w = 2
    end if
  end function cons_width

  !> Lagrange weights c(1:w) for interpolating/extrapolating a unit-spacing line
  !! sampled at local positions 0,1,...,w-1 to the (possibly out-of-range) target
  !! offset t.  Used to build the ghost values that close the symmetric stencil at
  !! domain boundaries.
  pure subroutine lagrange_unit_weights(w, t, c)
    integer, intent(in) :: w
    real(wp), intent(in) :: t
    real(wp), intent(out) :: c(:)
    integer :: m, k
    real(wp) :: p
    do m = 0, w - 1
      p = 1.0_wp
      do k = 0, w - 1
        if (k /= m) p = p * (t - real(k, wp)) / real(m - k, wp)
      end do
      c(m + 1) = p
    end do
  end subroutine lagrange_unit_weights

  !> High-order freestream-consistent CONSERVATIVE-FLUX face values of a
  !! unit-spacing nodal line g(1..n): hf(1..n+1), where face iface lies between
  !! node iface-1 and node iface.  The construction guarantees the telescoping
  !! identity
  !!
  !!     hf(i+1) - hf(i) = D[g](i)   (the matching high-order central derivative)
  !!
  !! EXACTLY (to round-off) at every node, interior AND boundary.  Because the
  !! nodal derivatives AND the transverse face metrics are BOTH built from this
  !! single operator, and central operators on the two computational axes commute
  !! discretely, the discrete metric identity Delta_xi(S_xi) + Delta_eta(S_eta)=0
  !! holds to round-off on ANY map — the high-order freestream-consistent metric
  !! of Visbal-Gaitonde (2002) / Nonomura-Iizuka-Fujii (2010) (design spec §5.3).
  !!
  !! Interior faces use the symmetric width-w conservative flux (w = cons_width):
  !!   * w=6: [1,-8,37,37,-8,1]/60, whose undivided difference is the 6th-order
  !!     central first derivative;
  !!   * w=4: [-1,7,7,-1]/12 (4th-order central);
  !!   * w=2: the classical [1,1]/2 (2nd-order central — the legacy operator).
  !! Near a boundary the symmetric stencil is closed by high-order polynomial
  !! extrapolation of the w nearest interior nodes into the w/2 ghost values it
  !! needs.  Extrapolation keeps the operator a FIXED LINEAR stencil (so the
  !! cross-axis operators still commute, hence freestream is exact) while keeping
  !! the closure high order; for the smooth body-fitted maps of interest this
  !! boundary closure does not degrade the global convergence order.
  pure subroutine cons_flux_1d_unit(g, hf)
    real(wp), intent(in) :: g(:)
    real(wp), intent(out) :: hf(:)
    integer :: n, w, half, iface, ll, s, m
    real(wp) :: cf(6), cl(6), acc
    real(wp) :: gp(1 - 3:size(g) + 3)
    n = size(g)
    w = cons_width(n)
    half = w / 2
    select case (w)
    case (6)
      cf(1:6) = [1.0_wp, -8.0_wp, 37.0_wp, 37.0_wp, -8.0_wp, 1.0_wp] / 60.0_wp
    case (4)
      cf(1:4) = [-1.0_wp, 7.0_wp, 7.0_wp, -1.0_wp] / 12.0_wp
    case default
      cf(1:2) = [0.5_wp, 0.5_wp]
    end select

    ! Interior nodes verbatim; ghost nodes by high-order Lagrange extrapolation of
    ! the w nearest boundary nodes, so every face can use the centred stencil.
    do s = 1, n
      gp(s) = g(s)
    end do
    do m = 1, half
      ! Left ghost at position 1-m, extrapolated from nodes 1..w (local 0..w-1).
      call lagrange_unit_weights(w, real(-m, wp), cl)
      acc = 0.0_wp
      do s = 1, w
        acc = acc + cl(s) * g(s)
      end do
      gp(1 - m) = acc
      ! Right ghost at position n+m, extrapolated from nodes n-w+1..n.
      call lagrange_unit_weights(w, real(m + w - 1, wp), cl)
      acc = 0.0_wp
      do s = 1, w
        acc = acc + cl(s) * g(n - w + s)
      end do
      gp(n + m) = acc
    end do

    ! Centred width-w conservative flux at every face.  Face iface sits between
    ! node iface-1 and node iface; the stencil spans nodes (iface-1)-(half-1) ..
    ! (iface-1)+half, i.e. local index ll..ll+w-1.
    do iface = 1, n + 1
      ll = iface - 1 - (half - 1)
      acc = 0.0_wp
      do s = 1, w
        acc = acc + cf(s) * gp(ll + s - 1)
      end do
      hf(iface) = acc
    end do
  end subroutine cons_flux_1d_unit

  !> Matching high-order first derivative of a unit-spacing nodal line: the
  !! undivided difference of the conservative-flux face values
  !! (cons_flux_1d_unit), so dg(i) reproduces the high-order central derivative
  !! AND stays bit-consistent with the face metrics built from the same operator
  !! (the property that makes freestream telescope to round-off, design spec §5.3).
  pure subroutine deriv_1d_unit(g, dg)
    real(wp), intent(in) :: g(:)
    real(wp), intent(out) :: dg(:)
    real(wp) :: hf(size(g) + 1)
    integer :: i, n
    n = size(g)
    call cons_flux_1d_unit(g, hf)
    do i = 1, n
      dg(i) = hf(i + 1) - hf(i)
    end do
  end subroutine deriv_1d_unit

  !> Compute NODAL transformation metrics for curvilinear FDM from the node
  !! coordinates already in mesh%x_node / mesh%y_node (shape (nx, ny), the Plot3D
  !! points treated as nodes; unit computational spacing dxi=deta=1).
  !!
  !! Allocates (with stat=) and fills:
  !!   * jac(nx, ny)         = x_xi*y_eta - x_eta*y_xi  at NODES (must be > 0)
  !!   * sx_xi(2, nx+1, ny)  = (y_eta, -x_eta)          at XI-faces
  !!   * sx_eta(2, nx, ny+1) = (-y_xi,  x_xi)           at ETA-faces
  !! and sets mesh%fdm_curvilinear = .true.
  !!
  !! Geometry derivatives AND the transverse face metrics are built from the
  !! SAME high-order conservative-flux operator (cons_flux_1d_unit: a symmetric
  !! width-w central conservative flux — 6th order for n>=6 — closed by high-order
  !! extrapolation at boundaries).  The nodal derivatives are the undivided
  !! difference of that operator's face values (deriv_1d_unit), and the face
  !! metric at face iface IS the operator's conservative-flux value of the
  !! transverse nodal derivative, so the telescoped first difference of the face
  !! metric reproduces the high-order nodal derivative EXACTLY at EVERY node
  !! (interior and boundary).  Because deriv_xi and deriv_eta are tensor-product
  !! operators on separate axes they commute discretely, so with this consistent
  !! HIGH-ORDER evaluation the discrete identities
  !!   Delta_xi(y_eta) - Delta_eta(y_xi) = 0,  Delta_xi(x_eta) - Delta_eta(x_xi) = 0
  !! hold to machine precision at every node even on a fully non-separable map —
  !! the constant flow telescopes to round-off (the freestream/GCL gate) WITHOUT
  !! capping the curved-grid order at 2 (the high-order metric of Visbal-Gaitonde
  !! 2002 / Nonomura-Iizuka-Fujii 2010, design spec §5.3).
  !!
  !! Face index convention: the xi-face between node i and node i+1 is stored at
  !! index i+1 in 1..nx+1, so index 1 is the left boundary face (node 1) and
  !! index nx+1 is the right boundary face (node nx).  Likewise for eta-faces.
  subroutine compute_metrics_fdm_2d(mesh, is_ok, message)
    type(mesh_2d_t), intent(inout) :: mesh
    logical, intent(out) :: is_ok
    character(len=*), intent(out) :: message

    integer :: nx, ny, i, j, astat
    real(wp), allocatable :: x_xi(:, :), x_eta(:, :), y_xi(:, :), y_eta(:, :)
    real(wp), allocatable :: line(:), dline(:), hfx(:), hfy(:)

    is_ok = .true.
    message = ''
    nx = mesh % nx; ny = mesh % ny

    if (nx < 3 .or. ny < 3) then
      is_ok = .false.
      message = 'mesh_2d: curvilinear FDM requires nx >= 3 and ny >= 3 nodes per axis'
      return
    end if
    if (.not. (allocated(mesh % x_node) .and. allocated(mesh % y_node))) then
      is_ok = .false.
      message = 'mesh_2d: compute_metrics_fdm_2d requires x_node/y_node'
      return
    end if

    allocate (x_xi(nx, ny), x_eta(nx, ny), y_xi(nx, ny), y_eta(nx, ny), stat=astat)
    if (astat /= 0) then
      is_ok = .false.; message = 'mesh_2d: allocation failed for nodal derivative scratch'; return
    end if

    ! xi-derivatives (along i) row by row.
    allocate (line(nx), dline(nx), stat=astat)
    if (astat /= 0) then
      is_ok = .false.; message = 'mesh_2d: allocation failed for xi-derivative line'; return
    end if
    do j = 1, ny
      line = mesh % x_node(:, j); call deriv_1d_unit(line, dline); x_xi(:, j) = dline
      line = mesh % y_node(:, j); call deriv_1d_unit(line, dline); y_xi(:, j) = dline
    end do
    deallocate (line, dline, stat=astat)
    if (astat /= 0) error stop 'mesh_2d: xi-derivative line deallocate failed'

    ! eta-derivatives (along j) column by column.
    allocate (line(ny), dline(ny), stat=astat)
    if (astat /= 0) then
      is_ok = .false.; message = 'mesh_2d: allocation failed for eta-derivative line'; return
    end if
    do i = 1, nx
      line = mesh % x_node(i, :); call deriv_1d_unit(line, dline); x_eta(i, :) = dline
      line = mesh % y_node(i, :); call deriv_1d_unit(line, dline); y_eta(i, :) = dline
    end do
    deallocate (line, dline, stat=astat)
    if (astat /= 0) error stop 'mesh_2d: eta-derivative line deallocate failed'

    ! Allocate the metric arrays (re-allocate defensively on a dimension change).
    ! Deallocate any stale array with checked stat= (do not silently overwrite the
    ! shared stat between the defensive deallocates and the re-allocate).
    if (allocated(mesh % jac)) then
      deallocate (mesh % jac, stat=astat)
      if (astat /= 0) then
        is_ok = .false.; message = 'mesh_2d: deallocation of stale jac failed'; return
      end if
    end if
    if (allocated(mesh % sx_xi)) then
      deallocate (mesh % sx_xi, stat=astat)
      if (astat /= 0) then
        is_ok = .false.; message = 'mesh_2d: deallocation of stale sx_xi failed'; return
      end if
    end if
    if (allocated(mesh % sx_eta)) then
      deallocate (mesh % sx_eta, stat=astat)
      if (astat /= 0) then
        is_ok = .false.; message = 'mesh_2d: deallocation of stale sx_eta failed'; return
      end if
    end if
    allocate (mesh % jac(nx, ny), mesh % sx_xi(2, nx + 1, ny), mesh % sx_eta(2, nx, ny + 1), stat=astat)
    if (astat /= 0) then
      is_ok = .false.; message = 'mesh_2d: allocation failed for nodal metric arrays'; return
    end if

    ! Jacobian at nodes (guard J > 0).
    do j = 1, ny
      do i = 1, nx
        mesh % jac(i, j) = x_xi(i, j) * y_eta(i, j) - x_eta(i, j) * y_xi(i, j)
        if (mesh % jac(i, j) <= 0.0_wp) then
          is_ok = .false.
          write (message, '(A,I0,A,I0,A)') 'mesh_2d: non-positive Jacobian at node (', &
            i, ',', j, ') — check node ordering / folded grid'
          return
        end if
      end do
    end do

    allocate (hfx(nx + 1), hfy(ny + 1), stat=astat)
    if (astat /= 0) then
      is_ok = .false.; message = 'mesh_2d: allocation failed for conservative-flux face scratch'; return
    end if

    ! xi-face metrics S_xi = (y_eta, -x_eta): the high-order conservative-flux
    ! interpolation (cons_flux_1d_unit) of the transverse nodal derivative ALONG i,
    ! evaluated at faces 1..nx+1.  Its telescoped difference is the high-order
    ! nodal xi-derivative, so freestream telescopes to round-off (see header).
    do j = 1, ny
      call cons_flux_1d_unit(y_eta(:, j), hfx); mesh % sx_xi(1, :, j) = hfx
      call cons_flux_1d_unit(x_eta(:, j), hfx); mesh % sx_xi(2, :, j) = -hfx
    end do

    ! eta-face metrics S_eta = (-y_xi, x_xi): the same operator ALONG j, evaluated
    ! at faces 1..ny+1.
    do i = 1, nx
      call cons_flux_1d_unit(y_xi(i, :), hfy); mesh % sx_eta(1, i, :) = -hfy
      call cons_flux_1d_unit(x_xi(i, :), hfy); mesh % sx_eta(2, i, :) = hfy
    end do

    mesh % fdm_curvilinear = .true.

    deallocate (x_xi, x_eta, y_xi, y_eta, hfx, hfy, stat=astat)
    if (astat /= 0) error stop 'mesh_2d: nodal derivative scratch deallocate failed'
  end subroutine compute_metrics_fdm_2d

  !> Build a curvilinear-FDM mesh from a single-block Plot3D grid: the Plot3D
  !! points are treated as NODES (mesh % nx = ni, mesh % ny = nj, nodes lying ON
  !! the boundaries), mesh % uniform = .false., dx/dy = 0 (unused on the
  !! curvilinear path), then compute_metrics_fdm_2d fills the nodal metrics.
  !! Distinct from build_mesh_2d_global's 'plot3d' branch, which treats the points
  !! as CELL CORNERS (nx = ni-1 cells) and builds finite-volume metrics.
  subroutine build_mesh_2d_fdm_curvilinear(mesh, grid_file, grid_binary, is_ok, message)
    type(mesh_2d_t), intent(inout) :: mesh
    character(len=*), intent(in) :: grid_file
    logical, intent(in) :: grid_binary
    logical, intent(out) :: is_ok
    character(len=*), intent(out) :: message

    integer :: ni, nj, astat
    real(wp), allocatable :: xg(:, :), yg(:, :)

    call read_plot3d_2d(grid_file, grid_binary, xg, yg, ni, nj, is_ok, message)
    if (.not. is_ok) return

    mesh % uniform = .false.
    mesh % nx = ni
    mesh % ny = nj
    mesh % dx_uniform = 0.0_wp
    mesh % dy_uniform = 0.0_wp
    call move_alloc(xg, mesh % x_node)
    call move_alloc(yg, mesh % y_node)

    ! Drop any stale finite-volume metric arrays so a re-used mesh_2d_t cannot
    ! carry FVM metrics into the curvilinear-FDM residual.
    if (allocated(mesh % vol)) deallocate (mesh % vol, stat=astat)
    if (allocated(mesh % s_xi)) deallocate (mesh % s_xi, stat=astat)
    if (allocated(mesh % s_eta)) deallocate (mesh % s_eta, stat=astat)
    if (allocated(mesh % xc)) deallocate (mesh % xc, stat=astat)
    if (allocated(mesh % yc)) deallocate (mesh % yc, stat=astat)

    call compute_metrics_fdm_2d(mesh, is_ok, message)
  end subroutine build_mesh_2d_fdm_curvilinear

  !> Build a uniform Cartesian NODAL mesh for the finite-difference (FDM) path.
  !! Unlike build_mesh_2d_global's 'uniform' branch (which stores CELL counts and
  !! derives FV metrics), this stores NODE counts directly: mesh % nx = npx,
  !! mesh % ny = npy, with grid points lying ON both boundaries. The spacing is
  !! dx = (x_right - x_left) / (npx - 1) (equal to the FV cell size for the same
  !! number of cells). The FV metric arrays (vol, s_xi, s_eta, xc, yc) are left
  !! UNALLOCATED on purpose — the uniform FDM residual divides by the scalar
  !! dx/dy, not by vol — so this builder does NOT call compute_metrics_2d.
  subroutine build_mesh_2d_uniform_nodal(mesh, npx, npy, x_left, x_right, y_left, y_right)
    type(mesh_2d_t), intent(inout) :: mesh
    integer, intent(in) :: npx, npy
    real(wp), intent(in) :: x_left, x_right, y_left, y_right
    integer :: i, j, astat

    if (npx < 2 .or. npy < 2) &
      error stop 'mesh_2d: build_mesh_2d_uniform_nodal requires npx >= 2 and npy >= 2'

    mesh % uniform = .true.
    mesh % nx = npx
    mesh % ny = npy
    mesh % dx_uniform = (x_right - x_left) / real(npx - 1, wp)
    mesh % dy_uniform = (y_right - y_left) / real(npy - 1, wp)

    if (allocated(mesh % x_node)) then
      deallocate (mesh % x_node, stat=astat)
      if (astat /= 0) error stop 'mesh_2d: x_node deallocate failed'
    end if
    if (allocated(mesh % y_node)) then
      deallocate (mesh % y_node, stat=astat)
      if (astat /= 0) error stop 'mesh_2d: y_node deallocate failed'
    end if
    ! FV metrics are left UNALLOCATED on purpose for the nodal mesh.  If this
    ! mesh_2d_t was previously initialised via build_mesh_2d_global +
    ! compute_metrics_2d, stale vol/s_xi/s_eta/xc/yc would remain allocated with
    ! the wrong dimensions.  Deallocate them so the nodal mesh is truly metric-free.
    if (allocated(mesh % vol)) then
      deallocate (mesh % vol, stat=astat)
      if (astat /= 0) error stop 'mesh_2d: vol deallocate failed'
    end if
    if (allocated(mesh % s_xi)) then
      deallocate (mesh % s_xi, stat=astat)
      if (astat /= 0) error stop 'mesh_2d: s_xi deallocate failed'
    end if
    if (allocated(mesh % s_eta)) then
      deallocate (mesh % s_eta, stat=astat)
      if (astat /= 0) error stop 'mesh_2d: s_eta deallocate failed'
    end if
    if (allocated(mesh % xc)) then
      deallocate (mesh % xc, stat=astat)
      if (astat /= 0) error stop 'mesh_2d: xc deallocate failed'
    end if
    if (allocated(mesh % yc)) then
      deallocate (mesh % yc, stat=astat)
      if (astat /= 0) error stop 'mesh_2d: yc deallocate failed'
    end if
    allocate (mesh % x_node(npx, npy), mesh % y_node(npx, npy), stat=astat)
    if (astat /= 0) error stop 'mesh_2d: allocation failed for nodal node arrays'

    do j = 1, npy
      do i = 1, npx
        mesh % x_node(i, j) = x_left + mesh % dx_uniform * real(i - 1, wp)
        mesh % y_node(i, j) = y_left + mesh % dy_uniform * real(j - 1, wp)
      end do
    end do
  end subroutine build_mesh_2d_uniform_nodal

  !> Slice the (global) metric arrays of `mesh` into this rank's local sub-block.
  !! Cells ix_first..ix_first+nx_local-1, iy_first..iy_first+ny_local-1 (1-based
  !! global). xi-faces span nx_local+1 columns starting at global column
  !! ix_first; eta-faces span ny_local+1 rows starting at iy_first. Shared
  !! inter-rank faces are identical on both ranks because both slice the same
  !! global metrics. No-op for uniform grids (they use scalar dx/dy).
  subroutine build_mesh_2d_local(mesh, ix_first, iy_first, nx_local, ny_local, is_ok, message)
    type(mesh_2d_t), intent(inout) :: mesh
    integer, intent(in) :: ix_first, iy_first, nx_local, ny_local
    logical, intent(out) :: is_ok
    character(len=*), intent(out) :: message
    real(wp), allocatable :: vol_l(:, :), xc_l(:, :), yc_l(:, :)
    real(wp), allocatable :: sxi_l(:, :, :), seta_l(:, :, :)
    real(wp), allocatable :: xn_l(:, :), yn_l(:, :), jac_l(:, :)
    integer :: i, j, astat, i0, j0

    is_ok = .true.
    message = ''
    if (mesh % uniform) return            ! uniform path uses scalar dx/dy; metrics unused

    if (ix_first < 1 .or. iy_first < 1 .or. &
        ix_first + nx_local - 1 > mesh % nx .or. iy_first + ny_local - 1 > mesh % ny) then
      is_ok = .false.
      message = 'mesh_2d: build_mesh_2d_local index range out of global bounds'
      return
    end if
    i0 = ix_first - 1; j0 = iy_first - 1

    if (mesh % fdm_curvilinear) then
      ! Curvilinear-FDM nodal-metric slice.  The Plot3D points are NODES, so
      ! mesh % nx/ny are NODE counts and the local block owns nx_local x ny_local
      ! NODES at global (ix_first, iy_first).  xi-faces number nx+1 (sx_xi), so the
      ! local block has nx_local+1 xi-faces starting at global face ix_first;
      ! eta-faces number ny+1.  Both ranks slice the SAME global metric, so a
      ! shared inter-rank face is bit-identical on the two ranks (the FVM
      ! curvilinear MPI invariant, applied to the nodal metrics).  Distinct from
      ! the FVM branch below, whose x_node array is nx+1 (cell corners); here
      ! x_node/y_node are node-count arrays sliced to (nx_local, ny_local).
      allocate (jac_l(nx_local, ny_local), &
                sxi_l(2, nx_local + 1, ny_local), seta_l(2, nx_local, ny_local + 1), &
                xn_l(nx_local, ny_local), yn_l(nx_local, ny_local), stat=astat)
      if (astat /= 0) then
        is_ok = .false.; message = 'mesh_2d: build_mesh_2d_local (fdm) allocation failed'; return
      end if
      do j = 1, ny_local
        do i = 1, nx_local
          jac_l(i, j) = mesh % jac(i0 + i, j0 + j)
          xn_l(i, j) = mesh % x_node(i0 + i, j0 + j)
          yn_l(i, j) = mesh % y_node(i0 + i, j0 + j)
        end do
      end do
      do j = 1, ny_local
        do i = 1, nx_local + 1
          sxi_l(:, i, j) = mesh % sx_xi(:, i0 + i, j0 + j)
        end do
      end do
      do j = 1, ny_local + 1
        do i = 1, nx_local
          seta_l(:, i, j) = mesh % sx_eta(:, i0 + i, j0 + j)
        end do
      end do
      call move_alloc(jac_l, mesh % jac)
      call move_alloc(sxi_l, mesh % sx_xi)
      call move_alloc(seta_l, mesh % sx_eta)
      call move_alloc(xn_l, mesh % x_node)
      call move_alloc(yn_l, mesh % y_node)
      mesh % nx = nx_local
      mesh % ny = ny_local
      return
    end if

    allocate (vol_l(nx_local, ny_local), xc_l(nx_local, ny_local), yc_l(nx_local, ny_local), &
              sxi_l(2, nx_local + 1, ny_local), seta_l(2, nx_local, ny_local + 1), &
              xn_l(nx_local + 1, ny_local + 1), yn_l(nx_local + 1, ny_local + 1), stat=astat)
    if (astat /= 0) then
      is_ok = .false.; message = 'mesh_2d: build_mesh_2d_local allocation failed'; return
    end if

    do j = 1, ny_local
      do i = 1, nx_local
        vol_l(i, j) = mesh % vol(i0 + i, j0 + j)
        xc_l(i, j) = mesh % xc(i0 + i, j0 + j)
        yc_l(i, j) = mesh % yc(i0 + i, j0 + j)
      end do
    end do
    ! Local node block: nodes ix_first..ix_first+nx_local (nx_local+1 of them).
    do j = 1, ny_local + 1
      do i = 1, nx_local + 1
        xn_l(i, j) = mesh % x_node(i0 + i, j0 + j)
        yn_l(i, j) = mesh % y_node(i0 + i, j0 + j)
      end do
    end do
    do j = 1, ny_local
      do i = 1, nx_local + 1
        sxi_l(:, i, j) = mesh % s_xi(:, i0 + i, j0 + j)
      end do
    end do
    do j = 1, ny_local + 1
      do i = 1, nx_local
        seta_l(:, i, j) = mesh % s_eta(:, i0 + i, j0 + j)
      end do
    end do

    call move_alloc(vol_l, mesh % vol)
    call move_alloc(xc_l, mesh % xc)
    call move_alloc(yc_l, mesh % yc)
    call move_alloc(sxi_l, mesh % s_xi)
    call move_alloc(seta_l, mesh % s_eta)
    call move_alloc(xn_l, mesh % x_node)
    call move_alloc(yn_l, mesh % y_node)
    mesh % nx = nx_local
    mesh % ny = ny_local
  end subroutine build_mesh_2d_local

end module mesh_2d