module maths
!-----------------------------------------------------------------------
! NAME
!     maths
!
! DESCRIPTION
!     The module contains all the mathematical subroutines used in the PEM.
!
! AUTHORS & DATE
!     L. Lange
!     JB Clement, 2023-2025
!
! NOTES
!     Adapted from Schorghofer MSIM (N.S, Icarus 2010)
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use numerics, only: dp, di, k4, eps

! DECLARATION
! -----------
implicit none

! PARAMETERS
! ----------
real(dp),    parameter :: pi = 4._dp*atan(1._dp) ! PI = 3.14159...
integer(di), parameter :: limiter_minmod  = 1
integer(di), parameter :: limiter_MC      = 2
integer(di), parameter :: limiter_vanLeer = 3

contains
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!=======================================================================
SUBROUTINE deriv1(z,nz,y,y0,ybot,dzY)
!-----------------------------------------------------------------------
! NAME
!     deriv1
!
! DESCRIPTION
!     Compute the first derivative of a function y(z) on an irregular grid.
!
! AUTHORS & DATE
!     N.S (Icarus 2010)
!     L. Lange
!
! NOTES
!     Upper boundary conditions: y(0) = y0
!     Lower boundary condition: yp = ybottom
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
integer(di),             intent(in)  :: nz       ! number of layer
real(dp), dimension(nz), intent(in)  :: z        ! depth layer
real(dp), dimension(nz), intent(in)  :: y        ! function which needs to be derived
real(dp),                intent(in)  :: y0, ybot ! boundary conditions
real(dp), dimension(nz), intent(out) :: dzY      ! derivative of y w.r.t depth

! LOCAL VARIABLES
! ---------------
integer(di) :: j
real(dp)    :: hm, hp, c1, c2, c3

! CODE
! ----
hp = z(2) - z(1)
c1 = z(1)/(hp*z(2))
c2 = 1/z(1) - 1/(z(2) - z(1))
c3 = -hp/(z(1)*z(2))
dzY(1) = c1*y(2) + c2*y(1) + c3*y0
do j = 2,nz - 1
    hp = z(j + 1) - z(j)
    hm = z(j) - z(j - 1)
    c1 = +hm/(hp*(z(j + 1) - z(j - 1)))
    c2 = 1._dp/hm - 1._dp/hp
    c3 = -hp/(hm*(z(j + 1) - z(j - 1)))
    dzY(j) = c1*y(j + 1) + c2*y(j) + c3*y(j - 1)
end do
dzY(nz) = (ybot - y(nz - 1))/(2._dp*(z(nz) - z(nz - 1)))

END SUBROUTINE deriv1
!=======================================================================

!=======================================================================
SUBROUTINE deriv2_simple(z,nz,y,y0,yNp1,yp2)
!-----------------------------------------------------------------------
! NAME
!     deriv2_simple
!
! DESCRIPTION
!     Compute the second derivative of a function y(z) on an irregular grid.
!
! AUTHORS & DATE
!     N.S (Icarus 2010)
!     L. Lange
!
! NOTES
!     Boundary conditions: y(0) = y0, y(nz + 1) = yNp1
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
integer(di),             intent(in)  :: nz
real(dp), dimension(nz), intent(in)  :: z, y
real(dp),                intent(in)  :: y0, yNp1
real(dp), dimension(nz), intent(out) :: yp2

! LOCAL VARIABLES
! ---------------
integer(di) :: j
real(dp)    :: hm, hp, c1, c2, c3

! CODE
! ----
c1 = +2._dp/((z(2) - z(1))*z(2))
c2 = -2._dp/((z(2) - z(1))*z(1))
c3 = +2._dp/(z(1)*z(2))
yp2(1) = c1*y(2) + c2*y(1) + c3*y0
do j = 2,nz - 1
    hp = z(j + 1) - z(j)
    hm = z(j) - z(j - 1)
    c1 = +2._dp/(hp*(z(j + 1) - z(j - 1)))
    c2 = -2._dp/(hp*hm)
    c3 = +2._dp/(hm*(z(j + 1) - z(j - 1)))
    yp2(j) = c1*y(j + 1) + c2*y(j) + c3*y(j - 1)
end do
yp2(nz) = (yNp1 - 2._dp*y(nz) + y(nz - 1))/(z(nz) - z(nz - 1))**2

END SUBROUTINE deriv2_simple
!=======================================================================

!=======================================================================
SUBROUTINE  deriv1_onesided(j,z,nz,y,dy_zj)
!-----------------------------------------------------------------------
! NAME
!     deriv1_onesided
!
! DESCRIPTION
!     First derivative of function y(z) at z(j) one-sided derivative
!     on irregular grid.
!
! AUTHORS & DATE
!     N.S (Icarus 2010)
!     L. Lange
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
integer(di),             intent(in)  :: nz, j
real(dp), dimension(nz), intent(in)  :: z, y
real(dp),                intent(out) :: dy_zj

! LOCAL VARIABLES
! ---------------
real(dp) :: h1, h2, c1, c2, c3

! CODE
! ----
if (j < 1 .or. j > nz - 2) then
    dy_zj = -1._dp
else
    h1 = z(j + 1) - z(j)
    h2 = z(j + 2)- z(j + 1)
    c1 = -(2._dp*h1 + h2)/(h1*(h1 + h2))
    c2 = (h1 + h2)/(h1*h2)
    c3 = -h1/(h2*(h1 + h2))
    dy_zj = c1*y(j) + c2*y(j + 1) + c3*y(j + 2)
end if

END SUBROUTINE deriv1_onesided
!=======================================================================

!=======================================================================
PURE SUBROUTINE colint(y,z,nz,i1,i2,integral)
!-----------------------------------------------------------------------
! NAME
!     colint
!
! DESCRIPTION
!     Column integrates y on irregular grid.
!
! AUTHORS & DATE
!     N.S (Icarus 2010)
!     L. Lange
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
integer(di),             intent(in)  :: nz, i1, i2
real(dp), dimension(nz), intent(in)  :: y, z
real(dp),                intent(out) :: integral

! LOCAL VARIABLES
! ---------------
integer(di)             :: i
real(dp), dimension(nz) :: dz

! CODE
! ----
dz(1) = (z(2) - 0.)/2
do i = 2,nz - 1
    dz(i) = (z(i + 1) - z(i - 1))/2.
end do
dz(nz) = z(nz) - z(nz - 1)
integral = sum(y(i1:i2)*dz(i1:i2))

END SUBROUTINE colint
!=======================================================================

!=======================================================================
SUBROUTINE findroot(y1,y2,z1,z2,zr)
!-----------------------------------------------------------------------
! NAME
!     findroot
!
! DESCRIPTION
!     Compute the root zr, between two values y1 and y2 at depth z1,z2.
!
! AUTHORS & DATE
!     L. Lange
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
real(dp), intent(in)  :: y1, y2
real(dp), intent(in)  :: z1, z2
real(dp), intent(out) :: zr

! CODE
! ----
zr = (y1*z2 - y2*z1)/(y1 - y2)

END SUBROUTINE findroot
!=======================================================================

!=======================================================================
SUBROUTINE solve_tridiag(a,b,c,d,n,x,error)
!-----------------------------------------------------------------------
! NAME
!     solve_tridiag
!
! DESCRIPTION
!     Solve a tridiagonal system Ax = d using the Thomas' algorithm
!     a: sub-diagonal
!     b: main diagonal
!     c: super-diagonal
!     d: right-hand side
!     x: solution
!
! AUTHORS & DATE
!     JB Clement, 2025
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
integer(di),                intent(in)  :: n
real(dp), dimension(n),     intent(in)  :: b, d
real(dp), dimension(n - 1), intent(in)  :: a, c
real(dp), dimension(n),     intent(out) :: x
integer(di),                intent(out) :: error

! LOCAL VARIABLES
! ---------------
integer(di)            :: i
real(dp)               :: m
real(dp), dimension(n) :: c_p, d_p

! CODE
! ----
! Check stability: diagonally dominant condition
error = 0
if (abs(b(1)) < abs(c(1))) then
    error = 1
    return
end if
do i = 2,n - 1
    if (abs(b(i)) < abs(a(i - 1)) + abs(c(i))) then
        error = 1
        return
    end if
end do
if (abs(b(n)) < abs(a(n - 1))) then
    error = 1
    return
end if

! Initialization
c_p(1) = c(1)/b(1)
d_p(1) = d(1)/b(1)

! Forward sweep
do i = 2,n - 1
    m = b(i) - a(i - 1)*c_p(i - 1)
    c_p(i) = c(i)/m
    d_p(i) = (d(i) - a(i - 1)*d_p(i - 1))/m
end do
m = b(n) - a(n - 1)*c_p(n - 1)
d_p(n) = (d(n) - a(n - 1)*d_p(n - 1))/m

! Backward substitution
x(n) = d_p(n)
do i = n - 1,1,-1
    x(i) = d_p(i) - c_p(i)*x(i + 1)
end do

END SUBROUTINE solve_tridiag
!=======================================================================

!=======================================================================
SUBROUTINE solve_steady_heat(n,z,mz,kappa,mkappa,T_left,q_right,T)
!-----------------------------------------------------------------------
! NAME
!     solve_steady_heat
!
! DESCRIPTION
!     Solve 1D steady-state heat equation with space-dependent thermal diffusivity.
!     Uses Thomas algorithm to solve tridiagonal system.
!     Left boundary: prescribed temperature (Dirichlet).
!     Right boundary: prescribed thermal flux (Neumann).
!
! AUTHORS & DATE
!     JB Clement, 2025
!
! NOTES
!     Grid points at z, mid-grid points at mz.
!     Thermal diffusivity specified at both grids (kappa, mkappa).
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use stoppage, only: stop_clean

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
integer(di),                intent(in)  :: n
real(dp), dimension(n),     intent(in)  :: z, kappa
real(dp), dimension(n - 1), intent(in)  :: mz, mkappa
real(dp),                   intent(in)  :: T_left, q_right
real(dp), dimension(n),     intent(out) :: T

! LOCAL VARIABLES
! ---------------
integer(di)                :: i, error
real(dp), dimension(n)     :: b, d
real(dp), dimension(n - 1) :: a, c

! CODE
! ----
! Initialization
a = 0._dp; b = 0._dp; c = 0._dp; d = 0._dp

! Left boundary condition (Dirichlet: prescribed temperature)
b(1) = 1._dp
d(1) = T_left

! Internal points
do i = 2,n - 1
    a(i - 1) = -mkappa(i - 1)/((mz(i) - mz(i - 1))*(z(i) - z(i - 1)))
    c(i) = -mkappa(i)/((mz(i) - mz(i - 1))*(z(i + 1) - z(i)))
    b(i) = -(a(i - 1) + c(i))
end do

! Right boundary condition (Neumann: prescribed temperature)
a(n - 1) = kappa(n - 1)/(z(n) - z(n - 1))
b(n) = -kappa(n)/(z(n) - z(n - 1))
d(n) = q_right

! Solve the tridiagonal linear system with the Thomas' algorithm
call solve_tridiag(a,b,c,d,n,T,error)
if (error /= 0) call stop_clean(__FILE__,__LINE__,"unstable solving!",1)

END SUBROUTINE solve_steady_heat
!=======================================================================

!=======================================================================
FUNCTION clamp(x,xmin,xmax) RESULT(y)
!-----------------------------------------------------------------------
! NAME
!     clamp
!
! DESCRIPTION
!     Clamp a scalar value into an admissible interval.
!
! AUTHORS & DATE
!     JB Clement, 04/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
real(dp), intent(in) :: x
real(dp), intent(in) :: xmin, xmax

! LOCAL VARIABLES
! ---------------
real(dp) :: y

! CODE
! ----
if (xmin > xmax) then
    y = x
else
    y = max(xmin,min(xmax,x))
end if

END FUNCTION clamp
!=======================================================================

!=======================================================================
SUBROUTINE solve_linear_system(a,b,x)
!-----------------------------------------------------------------------
! NAME
!     solve_linear_system
!
! DESCRIPTION
!     Solve a dense linear system by Gaussian elimination.
!
! AUTHORS & DATE
!     JB Clement, 04/2026
!
! NOTES
!     Uses partial pivoting and returns ok = .false. if singular.
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use stoppage, only: stop_clean

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
real(dp), dimension(:,:), intent(in)  :: a
real(dp), dimension(:),   intent(in)  :: b
real(dp), dimension(:),   intent(out) :: x

! LOCAL VARIABLES
! ---------------
integer(di)                           :: n, i, j, k, ipiv
real(dp)                              :: pivot, fac
real(dp), dimension(:),   allocatable :: rhs, rowtmp
real(dp), dimension(:,:), allocatable :: m

! CODE
! ----
! Initialization
n = size(b)
x(:) = 0._dp
if (size(a,1) /= n .or. size(a,2) /= n) call stop_clean(__FILE__,__LINE__,"size of a must be n x n!",1)
if (size(x) /= n) call stop_clean(__FILE__,__LINE__,"size of x must match size of b!",1)

allocate(m(n,n),rhs(n),rowtmp(n))
m(:,:) = a(:,:)
rhs(:) = b(:)

! Gaussian elimination with partial pivoting
do k = 1,n - 1
    ipiv = k
    pivot = abs(m(k,k))
    do i = k + 1,n
        if (abs(m(i,k)) > pivot) then
            pivot = abs(m(i,k))
            ipiv = i
        end if
    end do
    if (pivot <= eps) call stop_clean(__FILE__,__LINE__,"matrix is singular!",1)

    if (ipiv /= k) then
        rowtmp(:) = m(k,:)
        m(k,:) = m(ipiv,:)
        m(ipiv,:) = rowtmp(:)
        fac = rhs(k)
        rhs(k) = rhs(ipiv)
        rhs(ipiv) = fac
    end if

    do i = k + 1,n
        fac = m(i,k)/m(k,k)
        m(i,k:n) = m(i,k:n) - fac*m(k,k:n)
        rhs(i) = rhs(i) - fac*rhs(k)
    end do
end do

if (abs(m(n,n)) <= eps) call stop_clean(__FILE__,__LINE__,"matrix is singular!",1)

do i = n,1,-1
    x(i) = rhs(i)
    do j = i + 1,n
        x(i) = x(i) - m(i,j)*x(j)
    end do
    if (abs(m(i,i)) <= eps) call stop_clean(__FILE__,__LINE__,"matrix is singular!",1)
    x(i) = x(i)/m(i,i)
end do

END SUBROUTINE solve_linear_system
!=======================================================================

!=======================================================================
SUBROUTINE coef2prof(coef,x,fx)
!-----------------------------------------------------------------------
! NAME
!     coef2prof
!
! DESCRIPTION
!     Evaluate polynomial coefficients on a 1D profile.
!
! AUTHORS & DATE
!     JB Clement, 04/2026
!
! NOTES
!     Coefficients are ordered from degree 0 to degree n-1.
!-----------------------------------------------------------------------

! DEPENDENCIES
! -----------
use stoppage, only: stop_clean

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
real(dp), dimension(:), intent(in)  :: coef
real(dp), dimension(:), intent(in)  :: x
real(dp), dimension(:), intent(out) :: fx

! LOCAL VARIABLES
! ---------------
integer(di) :: i, j, ncoef
real(dp)    :: yloc

! CODE
! ----
! Initialization
fx(:) = 0._dp
if (size(fx) /= size(x)) call stop_clean(__FILE__,__LINE__,"size of fx must match size of x!",1)

ncoef = int(size(coef),di)

! Evaluate the polynomial at each point using Horner's method for numerical stability
do i = 1,size(x)
    yloc = coef(ncoef)
    do j = ncoef - 1,1,-1
        yloc = yloc*x(i) + coef(j)
    end do
    fx(i) = yloc
end do

END SUBROUTINE coef2prof
!=======================================================================

!=======================================================================
SUBROUTINE prof2coef(x,fx,coef)
!-----------------------------------------------------------------------
! NAME
!     prof2coef
!
! DESCRIPTION
!     Fit polynomial coefficients from a 1D profile by least squares.
!
! AUTHORS & DATE
!     JB Clement, 04/2026
!
! NOTES
!     Coefficients are ordered from degree 0 to degree n-1.
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use stoppage, only: stop_clean

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
real(dp), dimension(:), intent(in)  :: x
real(dp), dimension(:), intent(in)  :: fx
real(dp), dimension(:), intent(out) :: coef

! LOCAL VARIABLES
! ---------------
integer(di)                           :: i, j, k, ncoef
real(dp), dimension(:),   allocatable :: rhs, basis
real(dp), dimension(:,:), allocatable :: normal

! CODE
! ----
! Initialization
coef(:) = 0._dp
ncoef = int(size(coef),di)
if (size(fx) /= size(x)) call stop_clean(__FILE__,__LINE__,"size of fx must match size of x!",1)
if (size(x) < ncoef) call stop_clean(__FILE__,__LINE__,"size of x must be at least as large as the number of coefficients!",1)

allocate(rhs(ncoef),basis(ncoef),normal(ncoef,ncoef))
normal(:,:) = 0._dp
rhs(:) = 0._dp

! Set up the normal equations for least squares fitting
do i = 1,size(x)
    basis(1) = 1._dp
    do j = 2,ncoef
        basis(j) = basis(j - 1)*x(i)
    end do

    do j = 1,ncoef
        rhs(j) = rhs(j) + basis(j)*fx(i)
        do k = 1,ncoef
            normal(j,k) = normal(j,k) + basis(j)*basis(k)
        end do
    end do
end do

! Solve the linear system
call solve_linear_system(normal,rhs,coef)

END SUBROUTINE prof2coef
!=======================================================================

!=======================================================================
FUNCTION interp_profile(x_nodes,y_nodes,x_eval) RESULT(y_eval)
!-----------------------------------------------------------------------
! NAME
!     interp_profile
!
! DESCRIPTION
!     Interpolate a profile given at nodes.
!
! AUTHORS & DATE
!     JB Clement, 04/2026
!
! NOTES
!     Performs piecewise linear interpolation.
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
real(dp), dimension(:), intent(in) :: x_nodes, y_nodes
real(dp),               intent(in) :: x_eval

! LOCAL VARIABLES
! ---------------
real(dp)    :: y_eval, weight
integer(di) :: i, n

! CODE
! ----
! Initialization
n = int(size(x_nodes),di)
y_eval = 0._dp
if (n <= 0_di .or. size(y_nodes) /= n) return

! Handle out-of-bounds cases by clamping to the nearest node value
if (x_eval <= x_nodes(1)) then
    y_eval = y_nodes(1)
    return
end if
if (x_eval >= x_nodes(n)) then
    y_eval = y_nodes(n)
    return
end if

! Find the interval containing x_eval and perform linear interpolation
do i = 1,n - 1
    if (x_eval <= x_nodes(i + 1)) then
        if (abs(x_nodes(i + 1) - x_nodes(i)) <= eps) then ! Avoid division by zero if nodes are too close
            y_eval = y_nodes(i)
        else
            weight = (x_eval - x_nodes(i))/(x_nodes(i + 1) - x_nodes(i))
            y_eval = weight*(y_nodes(i + 1) - y_nodes(i)) + y_nodes(i)
        end if
        return
    end if
end do

y_eval = y_nodes(n)

END FUNCTION interp_profile
!=======================================================================

!=======================================================================
FUNCTION limit_slope(dqL,dqR,limiter) RESULT(limited_slope)
!-----------------------------------------------------------------------
! NAME
!     limit_slope
!
! DESCRIPTION
!     Applies a slope limiter based on the state on the left and on the
!     right.
!
! AUTHORS & DATE
!     JB Clement, 02/2026
!
! NOTES
!     Properties: TVD, monotone, linearity-preserving and compatible
!     with positivity rescaling.
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
real(dp),    intent(in) :: dqL, dqR
integer(di), intent(in) :: limiter

! LOCAL VARIABLES
! ---------------
real(dp) :: limited_slope

! CODE
! ----
if (dqL*dqR <= 0._dp) then
    limited_slope = 0._dp
    return
end if

select case (limiter)
    case(limiter_minmod) ! minmod: safe but strong diffusion
        limited_slope = sign(1._dp,dqL)*min(abs(dqL),abs(dqR))

    case(limiter_MC) ! MC: safe and less diffusion
        limited_slope = sign(1._dp,dqL)*min(2._dp*abs(dqL),0.5_dp*abs(dqL + dqR),2._dp*abs(dqR))

    case(limiter_vanLeer) ! van Leer: less diffusive, less grid sensitive but a bit more expensive
        limited_slope = (2._dp*dqL*dqR)/(dqL + dqR)

    case default ! Default is the MC limiter as it is a general-purpose choice
        limited_slope = sign(1._dp,dqL)*min(2._dp*abs(dqL),0.5_dp*abs(dqL + dqR),2._dp*abs(dqR))
end select

END FUNCTION limit_slope
!=======================================================================

!=======================================================================
FUNCTION enforce_positivity_slope(Qk,slope,dx_top,dx_bot) RESULT(enforced_slope)
!-----------------------------------------------------------------------
! NAME
!     enforce_positivity_slope
!
! DESCRIPTION
!     Enforce positivity through the slope.
!
! AUTHORS & DATE
!     JB Clement, 02/2026
!
! NOTES
!     For a linear function, the extrema inside the cell occur at the
!     interfaces.
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
real(dp), intent(in) :: Qk, slope, dx_top, dx_bot
real(dp)             :: enforced_slope

! LOCAL VARIABLES
! ---------------
real(dp) :: q1, q2, alpha

! CODE
! ----
q1 = Qk + slope*dx_top
q2 = Qk + slope*dx_bot

alpha = 1._dp
if (q1 < 0._dp) alpha = min(alpha,Qk/(Qk - q1))
if (q2 < 0._dp) alpha = min(alpha,Qk/(Qk - q2))

enforced_slope = alpha*slope

END FUNCTION enforce_positivity_slope
!=======================================================================

!=======================================================================
FUNCTION enforce_bounds_slope(Qk,slope,dxL,dxR,qmin,qmax) RESULT(enforced_slope)
!-----------------------------------------------------------------------
! NAME
!     enforce_bounds_slope
!
! DESCRIPTION
!     Enforce boundedness through the slope.
!
! AUTHORS & DATE
!     JB Clement, 02/s2026
!
! NOTES
!     For a linear function, the extrema inside the cell occur at the
!     interfaces.
!-----------------------------------------------------------------------

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
real(dp), intent(in) :: Qk, slope, dxL, dxR, qmin, qmax
real(dp)             :: enforced_slope

! LOCAL VARIABLES
! ---------------
real(dp) :: qL, qR, alpha

! CODE
! ----
qL = Qk + slope*dxL
qR = Qk + slope*dxR

alpha = 1._dp
if (qL > qmax) alpha = min(alpha,(qmax - Qk)/(qL - Qk))
if (qL < qmin) alpha = min(alpha,(qmin - Qk)/(qL - Qk))
if (qR > qmax) alpha = min(alpha,(qmax - Qk)/(qR - Qk))
if (qR < qmin) alpha = min(alpha,(qmin - Qk)/(qR - Qk))

! Ensure alpha is in [0,1]
if (alpha < 0._dp) alpha = 0._dp
if (alpha > 1._dp) alpha = 1._dp

enforced_slope = alpha*slope

END FUNCTION enforce_bounds_slope
!=======================================================================

!=======================================================================
SUBROUTINE remap_conserv_1d(Q_old,Q_new,x_old_i,x_new_i,n_old, n_new,order,limiter,enforce_pos,enforce_bounds,qmin,qmax)
!-----------------------------------------------------------------------
! NAME
!     remap_conserv_1d
!
! DESCRIPTION
!     1D conservative remapper, optional 1st/2nd order, with slope limiters,
!     positivity and boundedness enforcements.
!
! AUTHORS & DATE
!     JB Clement, 02/2026
!
! NOTES
!
!-----------------------------------------------------------------------

! DEPENDENCIES
! ------------
use stoppage, only: stop_clean

! DECLARATION
! -----------
implicit none

! ARGUMENTS
! ---------
integer(di),                    intent(in)           :: n_old, n_new
real(dp), dimension(n_old),     intent(in)           :: Q_old
real(dp), dimension(n_old + 1), intent(in)           :: x_old_i
real(dp), dimension(n_new + 1), intent(in)           :: x_new_i
integer(di),                    intent(in), optional :: order          ! 1 or 2 (default = 2)
integer(di),                    intent(in), optional :: limiter        ! See module parameters
logical(k4),                    intent(in), optional :: enforce_pos    ! Positivity (≥ 0)
logical(k4),                    intent(in), optional :: enforce_bounds ! Boundedness (qmin..qmax)
real(dp),                       intent(in), optional :: qmin, qmax     ! Required if enforce_bounds = .true.
real(dp), dimension(n_new),     intent(out)          :: Q_new

! LOCAL VARIABLES
! ---------------
integer(di)                :: ko, kn, order_loc, limiter_loc
logical(k4)                :: pos_loc, bounds_loc
real(dp)                   :: qmin_loc, qmax_loc, dqL, dqR, x1, x2, zL, zR, dx_overlap, dx_new, sumQ
real(dp), dimension(n_old) :: x_mid, dx_old, slope

! CODE
! ----
! Defaults
order_loc = 2
limiter_loc = limiter_MC
pos_loc = .false.
bounds_loc = .false.
qmin_loc = 0._dp
qmax_loc = 1._dp

if (present(order)) order_loc = order
if (present(limiter)) limiter_loc = limiter
if (present(enforce_pos)) pos_loc = enforce_pos
if (present(enforce_bounds)) bounds_loc = enforce_bounds
if (order_loc < 1) order_loc = 1
if (order_loc > 2) order_loc = 2
if (bounds_loc) then
    if (.not. present(qmin) .or. .not. present(qmax)) call stop_clean(__FILE__,__LINE__,"'enforce_bounds = .true.' requires qmin and qmax!",1)
    qmin_loc = qmin
    qmax_loc = qmax
end if

! Old grid geometry (midpoints and absolute layer widths)
do ko = 1, n_old
     x_mid(ko) = 0.5*(x_old_i(ko) + x_old_i(ko + 1))
     dx_old(ko) = abs(x_old_i(ko + 1) - x_old_i(ko))
end do

! Initialize slopes (zero for first order)
slope = 0._dp
if (order_loc == 2) then
    do ko = 1,n_old
        ! Centered difference estimates (avoid div by zero)
        if (ko == 1) then
            dqL = 0._dp
            if (abs(x_mid(2) - x_mid(1)) > eps) then
                dqR = (Q_old(2) - Q_old(1))/(x_mid(2) - x_mid(1))
            else
                dqR = 0._dp
            end if
        else if (ko == n_old) then
            if (abs(x_mid(n_old) - x_mid(n_old - 1)) > eps) then
                dqL = (Q_old(n_old) - Q_old(n_old - 1))/(x_mid(n_old) - x_mid(n_old - 1))
            else
                dqL = 0._dp
            end if
            dqR = 0._dp
        else
            if (abs(x_mid(ko) - x_mid(ko - 1)) > eps) then
                dqL = (Q_old(ko) - Q_old(ko - 1))/(x_mid(ko) - x_mid(ko - 1))
            else
                dqL = 0._dp
            end if
            if (abs(x_mid(ko + 1) - x_mid(ko)) > eps) then
                dqR = (Q_old(ko + 1) - Q_old(ko))/(x_mid(ko + 1) - x_mid(ko))
            else
                dqR = 0._dp
            end if
        end if

        ! Apply the chosen limiter
        slope(ko) = limit_slope(dqL,dqR,limiter_loc)

        ! Optional positivity (≥ 0)
        if (pos_loc) slope(ko) = enforce_positivity_slope(Q_old(ko),slope(ko),x_old_i(ko) - x_mid(ko),x_old_i(ko + 1) - x_mid(ko))

        ! Optional boundedness [qmin_loc,qmax_loc]
        if (bounds_loc) slope(ko) = enforce_bounds_slope(Q_old(ko),slope(ko),x_old_i(ko) - x_mid(ko),x_old_i(ko + 1) - x_mid(ko),qmin_loc,qmax_loc)
    end do
end if

! Conservative remap (exact integration of piecewise-linear reconstruction)
do kn = 1,n_new
    x1 = x_new_i(kn)
    x2 = x_new_i(kn + 1)
    dx_new = abs(x2 - x1)
    sumQ = 0._dp

    do ko = 1,n_old
        ! Compute overlap interval [zL,zR] robustly independent of coordinate ordering
        zL = max(min(x1,x2),min(x_old_i(ko),x_old_i(ko + 1)))
        zR = min(max(x1,x2),max(x_old_i(ko),x_old_i(ko + 1)))

        dx_overlap = zR - zL

        ! Integrate Q_old + slope*(x - x_mid) over [zL,zR]:
        if (dx_overlap > 0._dp) sumQ = sumQ + Q_old(ko)*dx_overlap + slope(ko)*0.5_dp*((zR - x_mid(ko))**2 - (zL - x_mid(ko))**2)
    end do

    ! Avoid division by zero (should not happen)
    if (abs(dx_new) > eps) then
        Q_new(kn) = sumQ/dx_new
    else
        Q_new(kn) = 0._dp
    end if
end do

END SUBROUTINE remap_conserv_1d
!=======================================================================

end module maths