!! WENO tracer face-reconstruction ladder for the layered tracer transport.
!!
!! Provides a drop-in replacement for PPM tracer advection, offering a
!! rung-adaptive reconstruction that automatically degrades to lower-order
!! schemes near boundaries where the full stencil is unavailable.
!!
!! Schemes (all with Z-weights, Borges et al. 2008):
!!   WENO5  — Jiang & Shu (1996),  3 quadratic candidates, nghost >= 3
!!   WENO7  — Balsara & Shu (2000), 4 cubic     candidates, nghost >= 4
!!   WENO9  — Balsara & Shu (2000), 5 quartic   candidates, nghost >= 5
!!
!! Caveat on WENO9 order: this rung uses simplified 3-term smoothness
!! indicators rather than the full quartic betas, so its effective
!! convergence is ~order 6 (not the nominal 9), and at a sharp front it
!! overshoots slightly MORE than WENO7 — WENO7 gives the best front L1.
!! Do not assume weno9 > weno7. See the function docstring and
!! docs/REFERENCE.md for the measured comparison.
!!
!! At interior cells far from boundaries the requested scheme is used.
!! As a face approaches the boundary the stencil is automatically degraded:
!!   WENO9 → WENO7 → WENO5 → PLM → donor-cell
!!
!! The swept-average form is used throughout: for upwind cell c with
!! CFL sigma, the face value integrates the reconstruction over
!! [c + 1/2 - sigma, c + 1/2], giving a time-averaged flux consistent
!! with the subcycled PPM approach in the rest of the solver.
module rdb_recon_weno
#ifdef LFORTRAN_PASSING
   use rdb_constants, only: wp, DRY_TOLERANCE
#else
   use rdb_constants, only: NZ_STACK_MAX, wp, DRY_TOLERANCE
#endif
   use rdb_grid, only: hgrid_t
   implicit none
   private

#ifdef LFORTRAN_PASSING
   integer, parameter :: NZ_STACK_MAX = 64
      !! LFortran 0.64 workaround: module-local copy of the rdb_constants value
      !! (an imported parameter used as an explicit-shape dummy bound inside a
      !! PURE call becomes an impure getter under LFortran). Keep in sync (=64).
#endif

   ! ---- Public integer parameters (rung codes) ----
   integer, parameter, public :: TRACER_RECON_PPM = 0
      !! Use PPM (the default path — does not call this module).
   integer, parameter, public :: TRACER_RECON_WENO5 = 1
      !! WENO5-Z (Jiang-Shu 1996 + Borges 2008 Z-weights), nghost >= 3.
   integer, parameter, public :: TRACER_RECON_WENO7 = 2
      !! WENO7-Z (Balsara-Shu 2000), nghost >= 4.
   integer, parameter, public :: TRACER_RECON_WENO9 = 3
      !! WENO9-Z (Balsara-Shu 2000), nghost >= 5.

   ! ---- Support-status codes (tracer_recon_support_status) ----
   integer, parameter, public :: TRACER_RECON_OK = 0
      !! Config honours the requested tracer_recon code.
   integer, parameter, public :: TRACER_RECON_REJECT_PQM = 1
      !! "pqm" requested — deferred, not implemented on any config.
   integer, parameter, public :: TRACER_RECON_REJECT_UNKNOWN = 2
      !! Unrecognised tracer_recon string.
   integer, parameter, public :: TRACER_RECON_REJECT_OCEAN = 3
      !! weno* requested on sim_type='ocean' (coastal-only in wave 1).
   integer, parameter, public :: TRACER_RECON_REJECT_NO_ML = 4
      !! weno* requested on a coastal run without multilayer tracers.

   ! ---- Public procedures ----
   public :: parse_tracer_recon
   public :: tracer_recon_required_nghost
   public :: tracer_recon_support_status
   public :: recon_rung_for_face
   ! Face helpers exported for the polynomial-exactness unit tests
   ! (test_tracer_recon); production code calls them only from the kernel.
   public :: plm_face_swept
   public :: weno5_face_swept
   public :: weno7_face_swept
   public :: weno9_face_swept

   ! The per-call concentration snapshot is caller-owned scratch
   ! (`tracer_scratch_t%weno_tr_snap`), passed in as an `intent(out)` dummy
   ! so `ml_advect_tracer_weno` stays `pure` — no module state, no per-step
   ! allocation (mirrors ml_advect_tracer_ppm; see the #315 refactor).

contains

   ! ==========================================================================
   !  Utility: parse / query
   ! ==========================================================================

   pure function parse_tracer_recon(str) result(code)
      !! Map a namelist string onto a TRACER_RECON_* code.
      !!
      !! Returns:
      !!   TRACER_RECON_PPM  (0) for "ppm"
      !!   TRACER_RECON_WENO5(1) for "weno5"
      !!   TRACER_RECON_WENO7(2) for "weno7"
      !!   TRACER_RECON_WENO9(3) for "weno9"
      !!   -2                    for "pqm" (deferred)
      !!   -1                    for anything else (unrecognised)
      character(len=*), intent(in) :: str
      integer :: code
      select case (trim(str))
      case ("ppm")
         code = TRACER_RECON_PPM
      case ("weno5")
         code = TRACER_RECON_WENO5
      case ("weno7")
         code = TRACER_RECON_WENO7
      case ("weno9")
         code = TRACER_RECON_WENO9
      case ("pqm")
         code = -2
      case default
         code = -1
      end select
   end function parse_tracer_recon

   pure function tracer_recon_required_nghost(code) result(ng)
      !! Minimum nghost required for a given TRACER_RECON_* code.
      integer, intent(in) :: code
         !! TRACER_RECON_* constant.
      integer :: ng
      select case (code)
      case (TRACER_RECON_WENO5)
         ng = 3
      case (TRACER_RECON_WENO7)
         ng = 4
      case (TRACER_RECON_WENO9)
         ng = 5
      case default
         ng = 2  ! PPM or unrecognised — PPM requires nghost >= 2
      end select
   end function tracer_recon_required_nghost

   pure function tracer_recon_support_status(recon_code, is_ocean, &
                                             has_multilayer_tracers) result(status)
      !! Decide whether a parsed tracer_recon code is honoured on this
      !! config, independent of the multilayer nghost/use_ppm_tracer
      !! quality guards (those apply only once this returns OK for a weno
      !! rung).  Factored out as a pure predicate so the fail-loud
      !! decision is unit-testable without triggering error stop.
      !!
      !! `recon_code` is the parse_tracer_recon result (may be < 0).
      !! The coastal multilayer solver is the ONLY consumer, so a weno*
      !! rung is rejected fail-loud on every other config:
      !!   * "pqm"        → REJECT_PQM      (deferred on all configs)
      !!   * unrecognised → REJECT_UNKNOWN
      !!   * sim_type='ocean' → REJECT_OCEAN  (coastal-only in wave 1)
      !!   * coastal, not multilayer → REJECT_NO_ML (2D barotropic / NH)
      !! "ppm" is always OK (the default path; bit-identical).
      integer, intent(in) :: recon_code
         !! parse_tracer_recon(cfg%tracer_recon) result.
      logical, intent(in) :: is_ocean
         !! .true. when sim_type == 'ocean'.
      logical, intent(in) :: has_multilayer_tracers
         !! cfg%use_multilayer — the coastal multilayer tracer registry.
      integer :: status
      if (recon_code == -2) then
         status = TRACER_RECON_REJECT_PQM
      else if (recon_code < 0) then
         status = TRACER_RECON_REJECT_UNKNOWN
      else if (recon_code == TRACER_RECON_PPM) then
         status = TRACER_RECON_OK
      else if (is_ocean) then
         status = TRACER_RECON_REJECT_OCEAN
      else if (.not. has_multilayer_tracers) then
         status = TRACER_RECON_REJECT_NO_ML
      else
         status = TRACER_RECON_OK
      end if
   end function tracer_recon_support_status

   ! ==========================================================================
   !  Rung-degradation helper
   ! ==========================================================================

   pure function recon_rung_for_face(avail_up, avail_down, rung_max) result(rung)
      !! Return the highest feasible reconstruction rung for one face.
      !!
      !! Internal rung scale (higher = more accurate):
      !!   4 = WENO9 (needs avail_up >= 5, avail_down >= 4)
      !!   3 = WENO7 (needs avail_up >= 4, avail_down >= 3)
      !!   2 = WENO5 (needs avail_up >= 3, avail_down >= 2)
      !!   1 = PLM   (needs avail_up >= 2, avail_down >= 1)
      !!   0 = donor (needs avail_up >= 1)
      !!  -1 = no data (should never happen at interior cells)
      !!
      !! rung_max is the PUBLIC TRACER_RECON_* code plus one:
      !!   WENO5(1) → rung_max_internal = 2
      !!   WENO7(2) → rung_max_internal = 3
      !!   WENO9(3) → rung_max_internal = 4
      !! The kernel maps the public code to rung_max with `recon + 1`.
      !$acc routine seq
      integer, intent(in) :: avail_up
         !! Interior cells from domain-west to face (including upwind cell c).
      integer, intent(in) :: avail_down
         !! Interior cells from face downwind side to domain-east.
      integer, intent(in) :: rung_max
         !! Maximum internal rung allowed (4 = WENO9, 3 = WENO7, 2 = WENO5).
      integer :: rung
      if (rung_max >= 4 .and. avail_up >= 5 .and. avail_down >= 4) then
         rung = 4
      else if (rung_max >= 3 .and. avail_up >= 4 .and. avail_down >= 3) then
         rung = 3
      else if (rung_max >= 2 .and. avail_up >= 3 .and. avail_down >= 2) then
         rung = 2
      else if (avail_up >= 2 .and. avail_down >= 1) then
         rung = 1
      else if (avail_up >= 1) then
         rung = 0
      else
         rung = -1
      end if
   end function recon_rung_for_face

   ! ==========================================================================
   !  Low-level face helpers (pure, !$acc routine seq)
   ! ==========================================================================

   pure function plm_face_swept(qm1, q0, qp1, sigma) result(face)
      !! PLM (piecewise-linear) swept-average face value.
      !!
      !! For u > 0 the upwind cell is q0; stencil: qm1 = i-1, q0 = i, qp1 = i+1.
      !! slope = minmod(q0 - qm1, qp1 - q0)
      !! face  = q0 + (1 - sigma)/2 * slope
      !!
      !! Reference: van Leer (1977) slope limiter.
      !$acc routine seq
      real(wp), intent(in) :: qm1, q0, qp1, sigma
      real(wp) :: face
      real(wp) :: d0, d1, slope
      d0 = q0 - qm1
      d1 = qp1 - q0
      if (d0*d1 > 0.0_wp) then
         if (abs(d0) < abs(d1)) then
            slope = d0
         else
            slope = d1
         end if
      else
         slope = 0.0_wp
      end if
      face = q0 + 0.5_wp*(1.0_wp - sigma)*slope
   end function plm_face_swept

   pure function weno5_face_swept(qm2, qm1, q0, qp1, qp2, sigma) result(face)
      !! WENO5-Z swept-average face value (u > 0, upwind cell = q0).
      !!
      !! Three quadratic candidates reconstructing the right edge of cell q0.
      !! Smoothness indicators: Jiang & Shu (1996), eq. 3.1.
      !! Z-weights: Borges et al. (2008).
      !!
      !! Stencil offsets relative to upwind cell i=q0:
      !!   r=0: {i-2, i-1, i}     d=1/10
      !!   r=1: {i-1, i,   i+1}   d=6/10
      !!   r=2: {i,   i+1, i+2}   d=3/10
      !$acc routine seq
      real(wp), intent(in) :: qm2, qm1, q0, qp1, qp2
         !! Cell averages at i-2, i-1, i, i+1, i+2.
      real(wp), intent(in) :: sigma
         !! CFL of upwind cell: |u_eff|*dt/dx.
      real(wp) :: face

      real(wp), parameter :: d0 = 0.1_wp, d1 = 0.6_wp, d2 = 0.3_wp
      real(wp), parameter :: eps = 1.0e-36_wp
      real(wp), parameter :: c13_12 = 13.0_wp/12.0_wp
      real(wp), parameter :: c1_4 = 0.25_wp

      ! Curvature/slope/intercept for each candidate (quadratic: A + B*xi + C*xi^2)
      ! Cell-centre xi=0, right edge xi=1/2.
      real(wp) :: C0, B0, A0
      real(wp) :: C1, B1, A1
      real(wp) :: C2, B2, A2

      ! Smoothness indicators
      real(wp) :: beta0, beta1, beta2, tau5

      ! Z-weights
      real(wp) :: alpha0, alpha1, alpha2, alpha_sum
      real(wp) :: omega0, omega1, omega2

      ! Swept-average face values per candidate
      real(wp) :: face0, face1, face2

      ! Candidate r=0: stencil {i-2, i-1, i}, offsets (-2,-1,0)
      C0 = 0.5_wp*(qm2 - 2.0_wp*qm1 + q0)
      B0 = C0 - (qm1 - q0)
      A0 = q0 - C0/12.0_wp

      ! Candidate r=1: stencil {i-1, i, i+1}, offsets (-1,0,+1)
      C1 = 0.5_wp*(qm1 - 2.0_wp*q0 + qp1)
      B1 = 0.5_wp*(qp1 - qm1)
      A1 = q0 - C1/12.0_wp

      ! Candidate r=2: stencil {i, i+1, i+2}, offsets (0,+1,+2)
      C2 = 0.5_wp*(q0 - 2.0_wp*qp1 + qp2)
      B2 = qp1 - q0 - C2
      A2 = q0 - C2/12.0_wp

      ! JS96 smoothness indicators
      beta0 = c13_12*(qm2 - 2.0_wp*qm1 + q0)**2 + c1_4*(qm2 - 4.0_wp*qm1 + 3.0_wp*q0)**2
      beta1 = c13_12*(qm1 - 2.0_wp*q0 + qp1)**2 + c1_4*(qm1 - qp1)**2
      beta2 = c13_12*(q0 - 2.0_wp*qp1 + qp2)**2 + c1_4*(3.0_wp*q0 - 4.0_wp*qp1 + qp2)**2

      ! Z-weights (Borges 2008): tau5 = |beta0 - beta2|
      tau5 = abs(beta0 - beta2)
      alpha0 = d0*(1.0_wp + (tau5/(eps + beta0))**2)
      alpha1 = d1*(1.0_wp + (tau5/(eps + beta1))**2)
      alpha2 = d2*(1.0_wp + (tau5/(eps + beta2))**2)
      alpha_sum = alpha0 + alpha1 + alpha2
      omega0 = alpha0/alpha_sum
      omega1 = alpha1/alpha_sum
      omega2 = alpha2/alpha_sum

      ! Swept-average of quadratic A + B*xi + C*xi^2 over [1/2-sigma, 1/2]:
      !   = (A + B/2 + C/4) - sigma/2*(B + C) + sigma^2/6*(2*C)
      face0 = (A0 + 0.5_wp*B0 + 0.25_wp*C0) - 0.5_wp*sigma*(B0 + C0) &
              + sigma**2/3.0_wp*C0
      face1 = (A1 + 0.5_wp*B1 + 0.25_wp*C1) - 0.5_wp*sigma*(B1 + C1) &
              + sigma**2/3.0_wp*C1
      face2 = (A2 + 0.5_wp*B2 + 0.25_wp*C2) - 0.5_wp*sigma*(B2 + C2) &
              + sigma**2/3.0_wp*C2

      face = omega0*face0 + omega1*face1 + omega2*face2
   end function weno5_face_swept

   pure function weno7_face_swept(qm3, qm2, qm1, q0, qp1, qp2, qp3, sigma) result(face)
      !! WENO7-Z swept-average face value (u > 0, upwind cell = q0).
      !!
      !! Four cubic candidates with optimal weights d=(4,18,12,1)/35.
      !! Z-weights: tau7 = |beta0 - beta3|.
      !!
      !! Coefficient matrices from Balsara & Shu (2000), Table 1.
      !! Smoothness indicators from Balsara & Shu (2000), eq. 2.17.
      !$acc routine seq
      real(wp), intent(in) :: qm3, qm2, qm1, q0, qp1, qp2, qp3
         !! Cell averages at i-3 through i+3.
      real(wp), intent(in) :: sigma
         !! CFL of upwind cell: |u_eff|*dt/dx.
      real(wp) :: face

      ! Optimal weights for WENO7 (BS00)
      real(wp), parameter :: d0_w = 4.0_wp/35.0_wp
      real(wp), parameter :: d1_w = 18.0_wp/35.0_wp
      real(wp), parameter :: d2_w = 12.0_wp/35.0_wp
      real(wp), parameter :: d3_w = 1.0_wp/35.0_wp
      real(wp), parameter :: eps = 1.0e-36_wp

      ! Cubic polynomial coefficients A + B*xi + C*xi^2 + D*xi^3
      ! for each of the 4 stencils (via exact rational matrices from BS00)
      real(wp) :: A0, B0, C0, D0
      real(wp) :: A1, B1, C1, D1
      real(wp) :: A2, B2, C2, D2
      real(wp) :: A3, B3, C3, D3

      ! Smoothness indicators (BS00 eq 2.17 exact)
      real(wp) :: beta0, beta1, beta2, beta3, tau7

      ! Z-weights
      real(wp) :: alpha0, alpha1, alpha2, alpha3, alpha_sum
      real(wp) :: omega0, omega1, omega2, omega3

      ! Swept-average face values per candidate
      real(wp) :: face0, face1, face2, face3

      ! Cubic swept-average helpers: p(1/2), p'(1/2), p''(1/2), p'''(1/2)
      real(wp) :: p_half, pp_half, ppp_half, pppp_half

      ! Stencil 0: offsets (-3,-2,-1,0) — q values: qm3,qm2,qm1,q0
      A0 = (1.0_wp/24.0_wp)*qm3 + (-1.0_wp/6.0_wp)*qm2 &
           + (5.0_wp/24.0_wp)*qm1 + (11.0_wp/12.0_wp)*q0
      B0 = (-7.0_wp/24.0_wp)*qm3 + (11.0_wp/8.0_wp)*qm2 &
           + (-23.0_wp/8.0_wp)*qm1 + (43.0_wp/24.0_wp)*q0
      C0 = (-1.0_wp/2.0_wp)*qm3 + 2.0_wp*qm2 &
           + (-5.0_wp/2.0_wp)*qm1 + 1.0_wp*q0
      D0 = (-1.0_wp/6.0_wp)*qm3 + (1.0_wp/2.0_wp)*qm2 &
           + (-1.0_wp/2.0_wp)*qm1 + (1.0_wp/6.0_wp)*q0

      ! Stencil 1: offsets (-2,-1,0,+1) — q values: qm2,qm1,q0,qp1
      A1 = 0.0_wp*qm2 + (-1.0_wp/24.0_wp)*qm1 &
           + (13.0_wp/12.0_wp)*q0 + (-1.0_wp/24.0_wp)*qp1
      B1 = (5.0_wp/24.0_wp)*qm2 + (-9.0_wp/8.0_wp)*qm1 &
           + (5.0_wp/8.0_wp)*q0 + (7.0_wp/24.0_wp)*qp1
      C1 = 0.0_wp*qm2 + (1.0_wp/2.0_wp)*qm1 &
           + (-1.0_wp)*q0 + (1.0_wp/2.0_wp)*qp1
      D1 = (-1.0_wp/6.0_wp)*qm2 + (1.0_wp/2.0_wp)*qm1 &
           + (-1.0_wp/2.0_wp)*q0 + (1.0_wp/6.0_wp)*qp1

      ! Stencil 2: offsets (-1,0,+1,+2) — q values: qm1,q0,qp1,qp2
      A2 = (-1.0_wp/24.0_wp)*qm1 + (13.0_wp/12.0_wp)*q0 &
           + (-1.0_wp/24.0_wp)*qp1 + 0.0_wp*qp2
      B2 = (-7.0_wp/24.0_wp)*qm1 + (-5.0_wp/8.0_wp)*q0 &
           + (9.0_wp/8.0_wp)*qp1 + (-5.0_wp/24.0_wp)*qp2
      C2 = (1.0_wp/2.0_wp)*qm1 + (-1.0_wp)*q0 &
           + (1.0_wp/2.0_wp)*qp1 + 0.0_wp*qp2
      D2 = (-1.0_wp/6.0_wp)*qm1 + (1.0_wp/2.0_wp)*q0 &
           + (-1.0_wp/2.0_wp)*qp1 + (1.0_wp/6.0_wp)*qp2

      ! Stencil 3: offsets (0,+1,+2,+3) — q values: q0,qp1,qp2,qp3
      A3 = (11.0_wp/12.0_wp)*q0 + (5.0_wp/24.0_wp)*qp1 &
           + (-1.0_wp/6.0_wp)*qp2 + (1.0_wp/24.0_wp)*qp3
      B3 = (-43.0_wp/24.0_wp)*q0 + (23.0_wp/8.0_wp)*qp1 &
           + (-11.0_wp/8.0_wp)*qp2 + (7.0_wp/24.0_wp)*qp3
      C3 = 1.0_wp*q0 + (-5.0_wp/2.0_wp)*qp1 &
           + 2.0_wp*qp2 + (-1.0_wp/2.0_wp)*qp3
      D3 = (-1.0_wp/6.0_wp)*q0 + (1.0_wp/2.0_wp)*qp1 &
           + (-1.0_wp/2.0_wp)*qp2 + (1.0_wp/6.0_wp)*qp3

      ! Smoothness indicators: the published integer-coefficient quadratic
      ! forms (Balsara & Shu 2000, eq. 2.17), one per candidate stencil.
      ! These are the forms validated in the numpy prototype
      ! (local_archive/prototypes/tracer_weno/) — do not substitute re-derivations.
      beta0 = qm3*(547.0_wp*qm3 - 3882.0_wp*qm2 + 4642.0_wp*qm1 - 1854.0_wp*q0) &
              + qm2*(7043.0_wp*qm2 - 17246.0_wp*qm1 + 7042.0_wp*q0) &
              + qm1*(11003.0_wp*qm1 - 9402.0_wp*q0) &
              + 2107.0_wp*q0**2
      beta1 = qm2*(267.0_wp*qm2 - 1642.0_wp*qm1 + 1602.0_wp*q0 - 494.0_wp*qp1) &
              + qm1*(2843.0_wp*qm1 - 5966.0_wp*q0 + 1922.0_wp*qp1) &
              + q0*(3443.0_wp*q0 - 2522.0_wp*qp1) &
              + 547.0_wp*qp1**2
      beta2 = qm1*(547.0_wp*qm1 - 2522.0_wp*q0 + 1922.0_wp*qp1 - 494.0_wp*qp2) &
              + q0*(3443.0_wp*q0 - 5966.0_wp*qp1 + 1602.0_wp*qp2) &
              + qp1*(2843.0_wp*qp1 - 1642.0_wp*qp2) &
              + 267.0_wp*qp2**2
      beta3 = q0*(2107.0_wp*q0 - 9402.0_wp*qp1 + 7042.0_wp*qp2 - 1854.0_wp*qp3) &
              + qp1*(11003.0_wp*qp1 - 17246.0_wp*qp2 + 4642.0_wp*qp3) &
              + qp2*(7043.0_wp*qp2 - 3882.0_wp*qp3) &
              + 547.0_wp*qp3**2

      tau7 = abs(beta0 - beta3)
      alpha0 = d0_w*(1.0_wp + (tau7/(eps + beta0))**2)
      alpha1 = d1_w*(1.0_wp + (tau7/(eps + beta1))**2)
      alpha2 = d2_w*(1.0_wp + (tau7/(eps + beta2))**2)
      alpha3 = d3_w*(1.0_wp + (tau7/(eps + beta3))**2)
      alpha_sum = alpha0 + alpha1 + alpha2 + alpha3
      omega0 = alpha0/alpha_sum
      omega1 = alpha1/alpha_sum
      omega2 = alpha2/alpha_sum
      omega3 = alpha3/alpha_sum

      ! Cubic swept-average over [1/2 - sigma, 1/2]:
      !   q_face = p(1/2) - sigma/2 * p'(1/2) + sigma^2/6 * p''(1/2) - sigma^3/24 * p'''(1/2)
      ! where:
      !   p(1/2)   = A + B/2 + C/4 + D/8
      !   p'(1/2)  = B + C + 3D/4
      !   p''(1/2) = 2C + 3D
      !   p'''(1/2)= 6D

      ! Stencil 0
      p_half = A0 + 0.5_wp*B0 + 0.25_wp*C0 + 0.125_wp*D0
      pp_half = B0 + C0 + 0.75_wp*D0
      ppp_half = 2.0_wp*C0 + 3.0_wp*D0
      pppp_half = 6.0_wp*D0
      face0 = p_half - 0.5_wp*sigma*pp_half &
              + sigma**2/6.0_wp*ppp_half - sigma**3/24.0_wp*pppp_half

      ! Stencil 1
      p_half = A1 + 0.5_wp*B1 + 0.25_wp*C1 + 0.125_wp*D1
      pp_half = B1 + C1 + 0.75_wp*D1
      ppp_half = 2.0_wp*C1 + 3.0_wp*D1
      pppp_half = 6.0_wp*D1
      face1 = p_half - 0.5_wp*sigma*pp_half &
              + sigma**2/6.0_wp*ppp_half - sigma**3/24.0_wp*pppp_half

      ! Stencil 2
      p_half = A2 + 0.5_wp*B2 + 0.25_wp*C2 + 0.125_wp*D2
      pp_half = B2 + C2 + 0.75_wp*D2
      ppp_half = 2.0_wp*C2 + 3.0_wp*D2
      pppp_half = 6.0_wp*D2
      face2 = p_half - 0.5_wp*sigma*pp_half &
              + sigma**2/6.0_wp*ppp_half - sigma**3/24.0_wp*pppp_half

      ! Stencil 3
      p_half = A3 + 0.5_wp*B3 + 0.25_wp*C3 + 0.125_wp*D3
      pp_half = B3 + C3 + 0.75_wp*D3
      ppp_half = 2.0_wp*C3 + 3.0_wp*D3
      pppp_half = 6.0_wp*D3
      face3 = p_half - 0.5_wp*sigma*pp_half &
              + sigma**2/6.0_wp*ppp_half - sigma**3/24.0_wp*pppp_half

      face = omega0*face0 + omega1*face1 + omega2*face2 + omega3*face3
   end function weno7_face_swept

   pure function weno9_face_swept(qm4, qm3, qm2, qm1, q0, qp1, qp2, qp3, qp4, sigma) &
      result(face)
      !! WENO9-Z swept-average face value (u > 0, upwind cell = q0).
      !!
      !! Five quartic candidates with optimal weights
      !! d = (1/126, 10/63, 10/21, 20/63, 5/126).
      !! Z-weights: tau9 = |beta0 - beta4|.
      !! Smoothness indicators use the simplified 3-term form.
      !!
      !! Reference: Balsara & Shu (2000).
      !$acc routine seq
      real(wp), intent(in) :: qm4, qm3, qm2, qm1, q0, qp1, qp2, qp3, qp4
         !! Cell averages at i-4 through i+4.
      real(wp), intent(in) :: sigma
         !! CFL of upwind cell: |u_eff|*dt/dx.
      real(wp) :: face

      ! Optimal weights (BS00): d = (1/126, 10/63, 10/21, 20/63, 5/126)
      real(wp), parameter :: d0_w = 1.0_wp/126.0_wp
      real(wp), parameter :: d1_w = 10.0_wp/63.0_wp
      real(wp), parameter :: d2_w = 10.0_wp/21.0_wp
      real(wp), parameter :: d3_w = 20.0_wp/63.0_wp
      real(wp), parameter :: d4_w = 5.0_wp/126.0_wp
      real(wp), parameter :: eps = 1.0e-36_wp
      real(wp), parameter :: c13_12 = 13.0_wp/12.0_wp
      real(wp), parameter :: c1_4 = 0.25_wp
      real(wp), parameter :: c1_80 = 1.0_wp/80.0_wp

      ! Quartic polynomial A + B*xi + C*xi^2 + D*xi^3 + E*xi^4 per stencil
      real(wp) :: A0, B0, C0, D0, E0
      real(wp) :: A1, B1, C1, D1, E1
      real(wp) :: A2, B2, C2, D2, E2
      real(wp) :: A3, B3, C3, D3, E3
      real(wp) :: A4, B4, C4, D4, E4

      ! Simplified 3-term beta (centred on stencil centre c):
      !   beta = (13/12)*(q_{c-1}-2*q_c+q_{c+1})^2 + (1/4)*(q_{c-1}-q_{c+1})^2
      !          + (1/80)*(q_{c-2}-4*q_{c-1}+6*q_c-4*q_{c+1}+q_{c+2})^2
      real(wp) :: beta0, beta1, beta2, beta3, beta4, tau9

      ! Z-weights
      real(wp) :: alpha0, alpha1, alpha2, alpha3, alpha4, alpha_sum
      real(wp) :: omega0, omega1, omega2, omega3, omega4

      ! Swept-average face values per candidate
      real(wp) :: face0, face1, face2, face3, face4

      ! Quartic swept-average helpers
      real(wp) :: p_half, pp_half, ppp_half, pppp_half, ppppp_half

      ! -----------------------------------------------------------------------
      ! Quartic coefficients (exact rational from BS00)
      ! -----------------------------------------------------------------------

      ! Stencil 0: offsets (-4,-3,-2,-1,0) → qm4,qm3,qm2,qm1,q0
      A0 = (-71.0_wp/1920.0_wp)*qm4 + (91.0_wp/480.0_wp)*qm3 &
           + (-373.0_wp/960.0_wp)*qm2 + (57.0_wp/160.0_wp)*qm1 &
           + (563.0_wp/640.0_wp)*q0
      B0 = (3.0_wp/16.0_wp)*qm4 + (-25.0_wp/24.0_wp)*qm3 &
           + (5.0_wp/2.0_wp)*qm2 + (-29.0_wp/8.0_wp)*qm1 &
           + (95.0_wp/48.0_wp)*q0
      C0 = (7.0_wp/16.0_wp)*qm4 + (-9.0_wp/4.0_wp)*qm3 &
           + (37.0_wp/8.0_wp)*qm2 + (-17.0_wp/4.0_wp)*qm1 &
           + (23.0_wp/16.0_wp)*q0
      D0 = (1.0_wp/4.0_wp)*qm4 + (-7.0_wp/6.0_wp)*qm3 &
           + 2.0_wp*qm2 + (-3.0_wp/2.0_wp)*qm1 &
           + (5.0_wp/12.0_wp)*q0
      E0 = (1.0_wp/24.0_wp)*qm4 + (-1.0_wp/6.0_wp)*qm3 &
           + (1.0_wp/4.0_wp)*qm2 + (-1.0_wp/6.0_wp)*qm1 &
           + (1.0_wp/24.0_wp)*q0

      ! Stencil 1: offsets (-3,-2,-1,0,+1) → qm3,qm2,qm1,q0,qp1
      A1 = (3.0_wp/640.0_wp)*qm3 + (-3.0_wp/160.0_wp)*qm2 &
           + (-13.0_wp/960.0_wp)*qm1 + (511.0_wp/480.0_wp)*q0 &
           + (-71.0_wp/1920.0_wp)*qp1
      B1 = (-5.0_wp/48.0_wp)*qm3 + (5.0_wp/8.0_wp)*qm2 &
           + (-7.0_wp/4.0_wp)*qm1 + (25.0_wp/24.0_wp)*q0 &
           + (3.0_wp/16.0_wp)*qp1
      C1 = (-1.0_wp/16.0_wp)*qm3 + (1.0_wp/4.0_wp)*qm2 &
           + (1.0_wp/8.0_wp)*qm1 + (-3.0_wp/4.0_wp)*q0 &
           + (7.0_wp/16.0_wp)*qp1
      D1 = (1.0_wp/12.0_wp)*qm3 + (-1.0_wp/2.0_wp)*qm2 &
           + 1.0_wp*qm1 + (-5.0_wp/6.0_wp)*q0 &
           + (1.0_wp/4.0_wp)*qp1
      E1 = (1.0_wp/24.0_wp)*qm3 + (-1.0_wp/6.0_wp)*qm2 &
           + (1.0_wp/4.0_wp)*qm1 + (-1.0_wp/6.0_wp)*q0 &
           + (1.0_wp/24.0_wp)*qp1

      ! Stencil 2: offsets (-2,-1,0,+1,+2) → qm2,qm1,q0,qp1,qp2
      A2 = (3.0_wp/640.0_wp)*qm2 + (-29.0_wp/480.0_wp)*qm1 &
           + (1067.0_wp/960.0_wp)*q0 + (-29.0_wp/480.0_wp)*qp1 &
           + (3.0_wp/640.0_wp)*qp2
      B2 = (5.0_wp/48.0_wp)*qm2 + (-17.0_wp/24.0_wp)*qm1 &
           + 0.0_wp*q0 + (17.0_wp/24.0_wp)*qp1 &
           + (-5.0_wp/48.0_wp)*qp2
      C2 = (-1.0_wp/16.0_wp)*qm2 + (3.0_wp/4.0_wp)*qm1 &
           + (-11.0_wp/8.0_wp)*q0 + (3.0_wp/4.0_wp)*qp1 &
           + (-1.0_wp/16.0_wp)*qp2
      D2 = (-1.0_wp/12.0_wp)*qm2 + (1.0_wp/6.0_wp)*qm1 &
           + 0.0_wp*q0 + (-1.0_wp/6.0_wp)*qp1 &
           + (1.0_wp/12.0_wp)*qp2
      E2 = (1.0_wp/24.0_wp)*qm2 + (-1.0_wp/6.0_wp)*qm1 &
           + (1.0_wp/4.0_wp)*q0 + (-1.0_wp/6.0_wp)*qp1 &
           + (1.0_wp/24.0_wp)*qp2

      ! Stencil 3: offsets (-1,0,+1,+2,+3) → qm1,q0,qp1,qp2,qp3
      A3 = (-71.0_wp/1920.0_wp)*qm1 + (511.0_wp/480.0_wp)*q0 &
           + (-13.0_wp/960.0_wp)*qp1 + (-3.0_wp/160.0_wp)*qp2 &
           + (3.0_wp/640.0_wp)*qp3
      B3 = (-3.0_wp/16.0_wp)*qm1 + (-25.0_wp/24.0_wp)*q0 &
           + (7.0_wp/4.0_wp)*qp1 + (-5.0_wp/8.0_wp)*qp2 &
           + (5.0_wp/48.0_wp)*qp3
      C3 = (7.0_wp/16.0_wp)*qm1 + (-3.0_wp/4.0_wp)*q0 &
           + (1.0_wp/8.0_wp)*qp1 + (1.0_wp/4.0_wp)*qp2 &
           + (-1.0_wp/16.0_wp)*qp3
      D3 = (-1.0_wp/4.0_wp)*qm1 + (5.0_wp/6.0_wp)*q0 &
           + (-1.0_wp)*qp1 + (1.0_wp/2.0_wp)*qp2 &
           + (-1.0_wp/12.0_wp)*qp3
      E3 = (1.0_wp/24.0_wp)*qm1 + (-1.0_wp/6.0_wp)*q0 &
           + (1.0_wp/4.0_wp)*qp1 + (-1.0_wp/6.0_wp)*qp2 &
           + (1.0_wp/24.0_wp)*qp3

      ! Stencil 4: offsets (0,+1,+2,+3,+4) → q0,qp1,qp2,qp3,qp4
      A4 = (563.0_wp/640.0_wp)*q0 + (57.0_wp/160.0_wp)*qp1 &
           + (-373.0_wp/960.0_wp)*qp2 + (91.0_wp/480.0_wp)*qp3 &
           + (-71.0_wp/1920.0_wp)*qp4
      B4 = (-95.0_wp/48.0_wp)*q0 + (29.0_wp/8.0_wp)*qp1 &
           + (-5.0_wp/2.0_wp)*qp2 + (25.0_wp/24.0_wp)*qp3 &
           + (-3.0_wp/16.0_wp)*qp4
      C4 = (23.0_wp/16.0_wp)*q0 + (-17.0_wp/4.0_wp)*qp1 &
           + (37.0_wp/8.0_wp)*qp2 + (-9.0_wp/4.0_wp)*qp3 &
           + (7.0_wp/16.0_wp)*qp4
      D4 = (-5.0_wp/12.0_wp)*q0 + (3.0_wp/2.0_wp)*qp1 &
           + (-2.0_wp)*qp2 + (7.0_wp/6.0_wp)*qp3 &
           + (-1.0_wp/4.0_wp)*qp4
      E4 = (1.0_wp/24.0_wp)*q0 + (-1.0_wp/6.0_wp)*qp1 &
           + (1.0_wp/4.0_wp)*qp2 + (-1.0_wp/6.0_wp)*qp3 &
           + (1.0_wp/24.0_wp)*qp4

      ! -----------------------------------------------------------------------
      ! Simplified 3-term smoothness indicators (centred on each stencil)
      ! Stencil r centre c:
      !   beta_r = (13/12)*(q_{c-1}-2*q_c+q_{c+1})^2
      !          + (1/4)  *(q_{c-1}-q_{c+1})^2
      !          + (1/80) *(q_{c-2}-4*q_{c-1}+6*q_c-4*q_{c+1}+q_{c+2})^2
      ! -----------------------------------------------------------------------
      ! Stencil 0: c = qm2 (index i-2), c±1 = qm3/qm1, c±2 = qm4/q0
      beta0 = c13_12*(qm3 - 2.0_wp*qm2 + qm1)**2 &
              + c1_4*(qm3 - qm1)**2 &
              + c1_80*(qm4 - 4.0_wp*qm3 + 6.0_wp*qm2 - 4.0_wp*qm1 + q0)**2

      ! Stencil 1: c = qm1 (index i-1), c±1 = qm2/q0, c±2 = qm3/qp1
      beta1 = c13_12*(qm2 - 2.0_wp*qm1 + q0)**2 &
              + c1_4*(qm2 - q0)**2 &
              + c1_80*(qm3 - 4.0_wp*qm2 + 6.0_wp*qm1 - 4.0_wp*q0 + qp1)**2

      ! Stencil 2: c = q0 (index i), c±1 = qm1/qp1, c±2 = qm2/qp2
      beta2 = c13_12*(qm1 - 2.0_wp*q0 + qp1)**2 &
              + c1_4*(qm1 - qp1)**2 &
              + c1_80*(qm2 - 4.0_wp*qm1 + 6.0_wp*q0 - 4.0_wp*qp1 + qp2)**2

      ! Stencil 3: c = qp1 (index i+1), c±1 = q0/qp2, c±2 = qm1/qp3
      beta3 = c13_12*(q0 - 2.0_wp*qp1 + qp2)**2 &
              + c1_4*(q0 - qp2)**2 &
              + c1_80*(qm1 - 4.0_wp*q0 + 6.0_wp*qp1 - 4.0_wp*qp2 + qp3)**2

      ! Stencil 4: c = qp2 (index i+2), c±1 = qp1/qp3, c±2 = q0/qp4
      beta4 = c13_12*(qp1 - 2.0_wp*qp2 + qp3)**2 &
              + c1_4*(qp1 - qp3)**2 &
              + c1_80*(q0 - 4.0_wp*qp1 + 6.0_wp*qp2 - 4.0_wp*qp3 + qp4)**2

      tau9 = abs(beta0 - beta4)
      alpha0 = d0_w*(1.0_wp + (tau9/(eps + beta0))**2)
      alpha1 = d1_w*(1.0_wp + (tau9/(eps + beta1))**2)
      alpha2 = d2_w*(1.0_wp + (tau9/(eps + beta2))**2)
      alpha3 = d3_w*(1.0_wp + (tau9/(eps + beta3))**2)
      alpha4 = d4_w*(1.0_wp + (tau9/(eps + beta4))**2)
      alpha_sum = alpha0 + alpha1 + alpha2 + alpha3 + alpha4
      omega0 = alpha0/alpha_sum
      omega1 = alpha1/alpha_sum
      omega2 = alpha2/alpha_sum
      omega3 = alpha3/alpha_sum
      omega4 = alpha4/alpha_sum

      ! Quartic swept-average over [1/2 - sigma, 1/2]:
      !   q_face = p(1/2) - sigma/2 * p'(1/2) + sigma^2/6 * p''(1/2)
      !            - sigma^3/24 * p'''(1/2) + sigma^4/120 * p''''(1/2)
      ! where for p = A + B*xi + C*xi^2 + D*xi^3 + E*xi^4:
      !   p(1/2)    = A + B/2 + C/4 + D/8 + E/16
      !   p'(1/2)   = B + C   + 3D/4 + E/2
      !   p''(1/2)  = 2C + 3D + 3E
      !   p'''(1/2) = 6D + 12E
      !   p''''(1/2)= 24E

      ! Stencil 0
      p_half = A0 + 0.5_wp*B0 + 0.25_wp*C0 + 0.125_wp*D0 + 0.0625_wp*E0
      pp_half = B0 + C0 + 0.75_wp*D0 + 0.5_wp*E0
      ppp_half = 2.0_wp*C0 + 3.0_wp*D0 + 3.0_wp*E0
      pppp_half = 6.0_wp*D0 + 12.0_wp*E0
      ppppp_half = 24.0_wp*E0
      face0 = p_half - 0.5_wp*sigma*pp_half + sigma**2/6.0_wp*ppp_half &
              - sigma**3/24.0_wp*pppp_half + sigma**4/120.0_wp*ppppp_half

      ! Stencil 1
      p_half = A1 + 0.5_wp*B1 + 0.25_wp*C1 + 0.125_wp*D1 + 0.0625_wp*E1
      pp_half = B1 + C1 + 0.75_wp*D1 + 0.5_wp*E1
      ppp_half = 2.0_wp*C1 + 3.0_wp*D1 + 3.0_wp*E1
      pppp_half = 6.0_wp*D1 + 12.0_wp*E1
      ppppp_half = 24.0_wp*E1
      face1 = p_half - 0.5_wp*sigma*pp_half + sigma**2/6.0_wp*ppp_half &
              - sigma**3/24.0_wp*pppp_half + sigma**4/120.0_wp*ppppp_half

      ! Stencil 2
      p_half = A2 + 0.5_wp*B2 + 0.25_wp*C2 + 0.125_wp*D2 + 0.0625_wp*E2
      pp_half = B2 + C2 + 0.75_wp*D2 + 0.5_wp*E2
      ppp_half = 2.0_wp*C2 + 3.0_wp*D2 + 3.0_wp*E2
      pppp_half = 6.0_wp*D2 + 12.0_wp*E2
      ppppp_half = 24.0_wp*E2
      face2 = p_half - 0.5_wp*sigma*pp_half + sigma**2/6.0_wp*ppp_half &
              - sigma**3/24.0_wp*pppp_half + sigma**4/120.0_wp*ppppp_half

      ! Stencil 3
      p_half = A3 + 0.5_wp*B3 + 0.25_wp*C3 + 0.125_wp*D3 + 0.0625_wp*E3
      pp_half = B3 + C3 + 0.75_wp*D3 + 0.5_wp*E3
      ppp_half = 2.0_wp*C3 + 3.0_wp*D3 + 3.0_wp*E3
      pppp_half = 6.0_wp*D3 + 12.0_wp*E3
      ppppp_half = 24.0_wp*E3
      face3 = p_half - 0.5_wp*sigma*pp_half + sigma**2/6.0_wp*ppp_half &
              - sigma**3/24.0_wp*pppp_half + sigma**4/120.0_wp*ppppp_half

      ! Stencil 4
      p_half = A4 + 0.5_wp*B4 + 0.25_wp*C4 + 0.125_wp*D4 + 0.0625_wp*E4
      pp_half = B4 + C4 + 0.75_wp*D4 + 0.5_wp*E4
      ppp_half = 2.0_wp*C4 + 3.0_wp*D4 + 3.0_wp*E4
      pppp_half = 6.0_wp*D4 + 12.0_wp*E4
      ppppp_half = 24.0_wp*E4
      face4 = p_half - 0.5_wp*sigma*pp_half + sigma**2/6.0_wp*ppp_half &
              - sigma**3/24.0_wp*pppp_half + sigma**4/120.0_wp*ppppp_half

      face = omega0*face0 + omega1*face1 + omega2*face2 + omega3*face3 + omega4*face4
   end function weno9_face_swept

end module rdb_recon_weno