module fvn_interpol
  use fvn_common
  implicit none
  
  
  ! Utility procedure find interval
  interface fvn_find_interval
      module procedure fvn_s_find_interval,fvn_d_find_interval
  end interface fvn_find_interval
  
  ! Quadratic 1D interpolation
  interface fvn_quad_interpol
      module procedure fvn_s_quad_interpol,fvn_d_quad_interpol
  end interface fvn_quad_interpol
  
  ! Quadratic 2D interpolation
  interface fvn_quad_2d_interpol
      module procedure fvn_s_quad_2d_interpol,fvn_d_quad_2d_interpol
  end interface fvn_quad_2d_interpol
  
  ! Quadratic 3D interpolation
  interface fvn_quad_3d_interpol
      module procedure fvn_s_quad_3d_interpol,fvn_d_quad_3d_interpol
  end interface fvn_quad_3d_interpol
  
  ! Akima interpolation
  interface fvn_akima
      module procedure fvn_s_akima,fvn_d_akima
  end interface fvn_akima
  
  ! Akima evaluation
  interface fvn_spline_eval
      module procedure fvn_s_spline_eval,fvn_d_spline_eval
  end interface fvn_spline_eval
  
  contains
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! Quadratic interpolation of tabulated function of 1,2 or 3 variables
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
  subroutine fvn_s_find_interval(x,i,xdata,n)
        implicit none
        ! This routine find the indice i where xdata(i) <= x < xdata(i+1)
        ! xdata(n) must contains a set of increasingly ordered values 
        ! if x < xdata(1) i=0 is returned
        ! if x > xdata(n) i=n is returned
        ! special case is where x=xdata(n) then n-1 is returned so 
        ! we will not exclude the upper limit
        ! a simple dichotomy method is used
  
        real(kind=4), intent(in) :: x

        integer(kind=4), intent(in) :: n

        real(kind=4), intent(in), dimension(n) :: xdata

        integer(kind=4), intent(out) :: i
  
        integer(kind=4) :: imin,imax,imoyen
  
        ! special case is where x=xdata(n) then n-1 is returned so 
        ! we will not exclude the upper limit
        if (x == xdata(n)) then
              i=n-1
              return
        end if
  
        ! if x < xdata(1) i=0 is returned
        if (x < xdata(1)) then
              i=0
              return
        end if
  
        ! if x > xdata(n) i=n is returned
        if (x > xdata(n)) then
              i=n
              return
        end if
  
        ! here xdata(1) <= x <= xdata(n)
        imin=0
        imax=n+1
  
        do while((imax-imin) > 1)
              imoyen=(imax+imin)/2
              if (x >= xdata(imoyen)) then
                    imin=imoyen
              else
                    imax=imoyen
              end if
        end do
  
        i=imin
  
  end subroutine
  
  
  subroutine fvn_d_find_interval(x,i,xdata,n)
        implicit none
        ! This routine find the indice i where xdata(i) <= x < xdata(i+1)
        ! xdata(n) must contains a set of increasingly ordered values 
        ! if x < xdata(1) i=0 is returned
        ! if x > xdata(n) i=n is returned
        ! special case is where x=xdata(n) then n-1 is returned so 
        ! we will not exclude the upper limit
        ! a simple dichotomy method is used
  
        real(kind=8), intent(in) :: x

        integer(kind=4), intent(in) :: n

        real(kind=8), intent(in), dimension(n) :: xdata

        integer(kind=4), intent(out) :: i
  
        integer(kind=4) :: imin,imax,imoyen
  
        ! special case is where x=xdata(n) then n-1 is returned so 
        ! we will not exclude the upper limit
        if (x == xdata(n)) then
              i=n-1
              return
        end if
  
        ! if x < xdata(1) i=0 is returned
        if (x < xdata(1)) then
              i=0
              return
        end if
  
        ! if x > xdata(n) i=n is returned
        if (x > xdata(n)) then
              i=n
              return
        end if
  
        ! here xdata(1) <= x <= xdata(n)
        imin=0
        imax=n+1
  
        do while((imax-imin) > 1)
              imoyen=(imax+imin)/2
              if (x >= xdata(imoyen)) then
                    imin=imoyen
              else
                    imax=imoyen
              end if
        end do
  
        i=imin
  
  end subroutine
  
  
  function fvn_s_quad_interpol(x,n,xdata,ydata)
        implicit none
        ! This function evaluate the value of a function defined by a set of points
        ! and values, using a quadratic interpolation
        ! xdata must be increasingly ordered
        ! x must be within xdata(1) and xdata(n) to actually do interpolation
        ! otherwise extrapolation is done
        integer(kind=4), intent(in) :: n
        real(kind=4), intent(in), dimension(n) :: xdata,ydata
        real(kind=4), intent(in) :: x
        real(kind=4) ::  fvn_s_quad_interpol
  
        integer(kind=4) :: iinf,base,i,j
        real(kind=4) :: p
  
        call fvn_s_find_interval(x,iinf,xdata,n)
  
        ! Settings for extrapolation
        if (iinf==0) then
              ! TODO -> Lower bound extrapolation warning
              iinf=1
        end if
  
        if (iinf==n) then
              ! TODO -> Higher bound extrapolation warning
              iinf=n-1
        end if
  
        ! The three points we will use are iinf-1,iinf and iinf+1 with the
        ! exception of the first interval, where iinf=1 we will use 1,2 and 3
        if (iinf==1) then
              base=0
        else
              base=iinf-2
        end if
  
        ! The three points we will use are :
        ! xdata/ydata(base+1),xdata/ydata(base+2),xdata/ydata(base+3)
  
        ! Straight forward Lagrange polynomial
        fvn_s_quad_interpol=0.
        do i=1,3
              !  polynome i
              p=ydata(base+i)
              do j=1,3
                    if (j /= i) then
                          p=p*(x-xdata(base+j))/(xdata(base+i)-xdata(base+j))
                    end if
              end do
              fvn_s_quad_interpol=fvn_s_quad_interpol+p
        end do
  
  end function
  
  
  function fvn_d_quad_interpol(x,n,xdata,ydata)
        implicit none
        ! This function evaluate the value of a function defined by a set of points
        ! and values, using a quadratic interpolation
        ! xdata must be increasingly ordered
        ! x must be within xdata(1) and xdata(n) to actually do interpolation
        ! otherwise extrapolation is done
        integer(kind=4), intent(in) :: n
        real(kind=8), intent(in), dimension(n) :: xdata,ydata
        real(kind=8), intent(in) :: x
        real(kind=8) ::  fvn_d_quad_interpol
  
        integer(kind=4) :: iinf,base,i,j
        real(kind=8) :: p
  
        call fvn_d_find_interval(x,iinf,xdata,n)
  
        ! Settings for extrapolation
        if (iinf==0) then
              ! TODO -> Lower bound extrapolation warning
              iinf=1
        end if
  
        if (iinf==n) then
              ! TODO Higher bound extrapolation warning
              iinf=n-1
        end if
  
        ! The three points we will use are iinf-1,iinf and iinf+1 with the
        ! exception of the first interval, where iinf=1 we will use 1,2 and 3
        if (iinf==1) then
              base=0
        else
              base=iinf-2
        end if
  
        ! The three points we will use are :
        ! xdata/ydata(base+1),xdata/ydata(base+2),xdata/ydata(base+3)
  
        ! Straight forward Lagrange polynomial
        fvn_d_quad_interpol=0.
        do i=1,3
              !  polynome i
              p=ydata(base+i)
              do j=1,3
                    if (j /= i) then
                          p=p*(x-xdata(base+j))/(xdata(base+i)-xdata(base+j))
                    end if
              end do
              fvn_d_quad_interpol=fvn_d_quad_interpol+p
        end do
  
  end function
  
  
  function fvn_s_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata)
        implicit none
        ! This function evaluate the value of a two variable function defined by a
        ! set of points and values, using a quadratic interpolation
        ! xdata and ydata must be increasingly ordered
        ! the couple (x,y) must be as x within xdata(1) and xdata(nx) and
        ! y within ydata(1) and ydata(ny) to actually do interpolation
        ! otherwise extrapolation is done
        integer(kind=4), intent(in) :: nx,ny
        real(kind=4), intent(in) :: x,y
        real(kind=4), intent(in), dimension(nx) :: xdata
        real(kind=4), intent(in), dimension(ny) :: ydata
        real(kind=4), intent(in), dimension(nx,ny) :: zdata
        real(kind=4) :: fvn_s_quad_2d_interpol
  
        integer(kind=4) :: ixinf,iyinf,basex,basey,i
        real(kind=4),dimension(3) :: ztmp
        !real(kind=4), external :: fvn_s_quad_interpol
  
        call fvn_s_find_interval(x,ixinf,xdata,nx)
        call fvn_s_find_interval(y,iyinf,ydata,ny)
  
        ! Settings for extrapolation
        if (ixinf==0) then
              ! TODO -> Lower x  bound extrapolation warning
              ixinf=1
        end if
  
        if (ixinf==nx) then
              ! TODO -> Higher x bound extrapolation warning
              ixinf=nx-1
        end if
  
        if (iyinf==0) then
              ! TODO -> Lower y  bound extrapolation warning
              iyinf=1
        end if
  
        if (iyinf==ny) then
              ! TODO -> Higher y bound extrapolation warning
              iyinf=ny-1
        end if
  
        ! The three points we will use are iinf-1,iinf and iinf+1 with the
        ! exception of the first interval, where iinf=1 we will use 1,2 and 3
        if (ixinf==1) then
              basex=0
        else
              basex=ixinf-2
        end if
  
        if (iyinf==1) then
              basey=0
        else
              basey=iyinf-2
        end if
  
        ! First we make 3 interpolations for x at y(base+1),y(base+2),y(base+3)
        ! stored in ztmp(1:3)
        do i=1,3
              ztmp(i)=fvn_s_quad_interpol(x,nx,xdata,zdata(:,basey+i))
        end do
  
        ! Then we make an interpolation for y using previous interpolations
        fvn_s_quad_2d_interpol=fvn_s_quad_interpol(y,3,ydata(basey+1:basey+3),ztmp)
  end function
  
  
  function fvn_d_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata)
        implicit none
        ! This function evaluate the value of a two variable function defined by a
        ! set of points and values, using a quadratic interpolation
        ! xdata and ydata must be increasingly ordered
        ! the couple (x,y) must be as x within xdata(1) and xdata(nx) and
        ! y within ydata(1) and ydata(ny) to actually do interpolation
        ! otherwise extrapolation is done
        integer(kind=4), intent(in) :: nx,ny
        real(kind=8), intent(in) :: x,y
        real(kind=8), intent(in), dimension(nx) :: xdata
        real(kind=8), intent(in), dimension(ny) :: ydata
        real(kind=8), intent(in), dimension(nx,ny) :: zdata
        real(kind=8) :: fvn_d_quad_2d_interpol
  
        integer(kind=4) :: ixinf,iyinf,basex,basey,i
        real(kind=8),dimension(3) :: ztmp
        !real(kind=8), external :: fvn_d_quad_interpol
  
        call fvn_d_find_interval(x,ixinf,xdata,nx)
        call fvn_d_find_interval(y,iyinf,ydata,ny)
  
        ! Settings for extrapolation
        if (ixinf==0) then
              ! TODO -> Lower x  bound extrapolation warning
              ixinf=1
        end if
  
        if (ixinf==nx) then
              ! TODO -> Higher x bound extrapolation warning
              ixinf=nx-1
        end if
  
        if (iyinf==0) then
              ! TODO -> Lower y  bound extrapolation warning
              iyinf=1
        end if
  
        if (iyinf==ny) then
              ! TODO -> Higher y bound extrapolation warning
              iyinf=ny-1
        end if
  
        ! The three points we will use are iinf-1,iinf and iinf+1 with the
        ! exception of the first interval, where iinf=1 we will use 1,2 and 3
        if (ixinf==1) then
              basex=0
        else
              basex=ixinf-2
        end if
  
        if (iyinf==1) then
              basey=0
        else
              basey=iyinf-2
        end if
  
        ! First we make 3 interpolations for x at y(base+1),y(base+2),y(base+3)
        ! stored in ztmp(1:3)
        do i=1,3
              ztmp(i)=fvn_d_quad_interpol(x,nx,xdata,zdata(:,basey+i))
        end do
  
        ! Then we make an interpolation for y using previous interpolations
        fvn_d_quad_2d_interpol=fvn_d_quad_interpol(y,3,ydata(basey+1:basey+3),ztmp)
  end function
  
  
  function fvn_s_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata)
        implicit none
        ! This function evaluate the value of a 3 variables function defined by a
        ! set of points and values, using a quadratic interpolation
        ! xdata, ydata and zdata must be increasingly ordered
        ! The triplet (x,y,z) must be within xdata,ydata and zdata to actually 
        ! perform an interpolation, otherwise extrapolation is done
        integer(kind=4), intent(in) :: nx,ny,nz
        real(kind=4), intent(in) :: x,y,z
        real(kind=4), intent(in), dimension(nx) :: xdata
        real(kind=4), intent(in), dimension(ny) :: ydata
        real(kind=4), intent(in), dimension(nz) :: zdata
        real(kind=4), intent(in), dimension(nx,ny,nz) :: tdata
        real(kind=4) :: fvn_s_quad_3d_interpol
  
        integer(kind=4) :: ixinf,iyinf,izinf,basex,basey,basez,i,j
        !real(kind=4), external :: fvn_s_quad_interpol,fvn_s_quad_2d_interpol
        real(kind=4),dimension(3,3) :: ttmp
  
        call fvn_s_find_interval(x,ixinf,xdata,nx)
        call fvn_s_find_interval(y,iyinf,ydata,ny)
        call fvn_s_find_interval(z,izinf,zdata,nz)
  
        ! Settings for extrapolation
        if (ixinf==0) then
              ! TODO -> Lower x  bound extrapolation warning
              ixinf=1
        end if
  
        if (ixinf==nx) then
              ! TODO -> Higher x bound extrapolation warning
              ixinf=nx-1
        end if
  
        if (iyinf==0) then
              ! TODO -> Lower y  bound extrapolation warning
              iyinf=1
        end if
  
        if (iyinf==ny) then
              ! TODO -> Higher y bound extrapolation warning
              iyinf=ny-1
        end if
  
        if (izinf==0) then
              ! TODO -> Lower z bound extrapolation warning
              izinf=1
        end if
  
        if (izinf==nz) then
              ! TODO -> Higher z bound extrapolation warning
              izinf=nz-1
        end if
  
        ! The three points we will use are iinf-1,iinf and iinf+1 with the
        ! exception of the first interval, where iinf=1 we will use 1,2 and 3
        if (ixinf==1) then
              basex=0
        else
              basex=ixinf-2
        end if
  
        if (iyinf==1) then
              basey=0
        else
              basey=iyinf-2
        end if
  
        if (izinf==1) then
              basez=0
        else
              basez=izinf-2
        end if
  
        ! We first make 9 one dimensional interpolation on variable x.
        ! results are stored in ttmp
        do i=1,3
              do j=1,3
                    ttmp(i,j)=fvn_s_quad_interpol(x,nx,xdata,tdata(:,basey+i,basez+j))
              end do
        end do
  
        ! We then make a 2 dimensionnal interpolation on variables y and z
        fvn_s_quad_3d_interpol=fvn_s_quad_2d_interpol(y,z, &
              3,ydata(basey+1:basey+3),3,zdata(basez+1:basez+3),ttmp)
  end function
  
  
  function fvn_d_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata)
        implicit none
        ! This function evaluate the value of a 3 variables function defined by a
        ! set of points and values, using a quadratic interpolation
        ! xdata, ydata and zdata must be increasingly ordered
        ! The triplet (x,y,z) must be within xdata,ydata and zdata to actually 
        ! perform an interpolation, otherwise extrapolation is done
        integer(kind=4), intent(in) :: nx,ny,nz
        real(kind=8), intent(in) :: x,y,z
        real(kind=8), intent(in), dimension(nx) :: xdata
        real(kind=8), intent(in), dimension(ny) :: ydata
        real(kind=8), intent(in), dimension(nz) :: zdata
        real(kind=8), intent(in), dimension(nx,ny,nz) :: tdata
        real(kind=8) :: fvn_d_quad_3d_interpol
  
        integer(kind=4) :: ixinf,iyinf,izinf,basex,basey,basez,i,j
        !real(kind=8), external :: fvn_d_quad_interpol,fvn_d_quad_2d_interpol
        real(kind=8),dimension(3,3) :: ttmp
  
        call fvn_d_find_interval(x,ixinf,xdata,nx)
        call fvn_d_find_interval(y,iyinf,ydata,ny)
        call fvn_d_find_interval(z,izinf,zdata,nz)
  
        ! Settings for extrapolation
        if (ixinf==0) then
              ! TODO -> Lower x  bound extrapolation warning
              ixinf=1
        end if
  
        if (ixinf==nx) then
              ! TODO -> Higher x bound extrapolation warning
              ixinf=nx-1
        end if
  
        if (iyinf==0) then
              ! TODO -> Lower y  bound extrapolation warning
              iyinf=1
        end if
  
        if (iyinf==ny) then
              ! TODO -> Higher y bound extrapolation warning
              iyinf=ny-1
        end if
  
        if (izinf==0) then
              ! TODO -> Lower z bound extrapolation warning
              izinf=1
        end if
  
        if (izinf==nz) then
              ! TODO -> Higher z bound extrapolation warning
              izinf=nz-1
        end if
  
        ! The three points we will use are iinf-1,iinf and iinf+1 with the
        ! exception of the first interval, where iinf=1 we will use 1,2 and 3
        if (ixinf==1) then
              basex=0
        else
              basex=ixinf-2
        end if
  
        if (iyinf==1) then
              basey=0
        else
              basey=iyinf-2
        end if
  
        if (izinf==1) then
              basez=0
        else
              basez=izinf-2
        end if
  
        ! We first make 9 one dimensional interpolation on variable x.
        ! results are stored in ttmp
        do i=1,3
              do j=1,3
                    ttmp(i,j)=fvn_d_quad_interpol(x,nx,xdata,tdata(:,basey+i,basez+j))
              end do
        end do
  
        ! We then make a 2 dimensionnal interpolation on variables y and z
        fvn_d_quad_3d_interpol=fvn_d_quad_2d_interpol(y,z, &
              3,ydata(basey+1:basey+3),3,zdata(basez+1:basez+3),ttmp)
  end function
  
  
  
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! Akima spline interpolation and spline evaluation
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
  ! Single precision
  subroutine fvn_s_akima(n,x,y,br,co)
      implicit none
      integer, intent(in)  :: n
      real, intent(in) :: x(n)
      real, intent(in) :: y(n)
      real, intent(out) :: br(n)
      real, intent(out) :: co(4,n)
      
      real, allocatable :: var(:),z(:)
      real :: wi_1,wi
      integer :: i
      real :: dx,a,b
  
      ! br is just a copy of x
      br(:)=x(:)
      
      allocate(var(n+3))
      allocate(z(n))
      ! evaluate the variations
      do i=1, n-1
          var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))
      end do
      var(n+2)=2.e0*var(n+1)-var(n)
      var(n+3)=2.e0*var(n+2)-var(n+1)
      var(2)=2.e0*var(3)-var(4)
      var(1)=2.e0*var(2)-var(3)
    
      do i = 1, n
      wi_1=abs(var(i+3)-var(i+2))
      wi=abs(var(i+1)-var(i))
      if ((wi_1+wi).eq.0.e0) then
          z(i)=(var(i+2)+var(i+1))/2.e0
      else
          z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi)
      end if
      end do
      
      do i=1, n-1
          dx=x(i+1)-x(i)
          a=(z(i+1)-z(i))*dx                      ! coeff intermediaires pour calcul wd
          b=y(i+1)-y(i)-z(i)*dx                   ! coeff intermediaires pour calcul wd
          co(1,i)=y(i)
          co(2,i)=z(i)
          !co(3,i)=-(a-3.*b)/dx**2                ! méthode wd
          !co(4,i)=(a-2.*b)/dx**3                 ! méthode wd
          co(3,i)=(3.e0*var(i+2)-2.e0*z(i)-z(i+1))/dx   ! méthode JP Moreau
          co(4,i)=(z(i)+z(i+1)-2.e0*var(i+2))/dx**2  !
          ! les coefficients donnés par imsl sont co(3,i)*2 et co(4,i)*6
          ! etrangement la fonction csval corrige et donne la bonne valeur ...
      end do
      co(1,n)=y(n)
      co(2,n)=z(n)
      co(3,n)=0.e0
      co(4,n)=0.e0
  
      deallocate(z)
      deallocate(var)
  
  end subroutine
  
  ! Double precision
  subroutine fvn_d_akima(n,x,y,br,co)
  
      implicit none
      integer, intent(in)  :: n
      double precision, intent(in) :: x(n)
      double precision, intent(in) :: y(n)
      double precision, intent(out) :: br(n)
      double precision, intent(out) :: co(4,n)
      
      double precision, allocatable :: var(:),z(:)
      double precision :: wi_1,wi
      integer :: i
      double precision :: dx,a,b
      
      ! br is just a copy of x
      br(:)=x(:)
  
      allocate(var(n+3))
      allocate(z(n))
      ! evaluate the variations
      do i=1, n-1
          var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))
      end do
      var(n+2)=2.d0*var(n+1)-var(n)
      var(n+3)=2.d0*var(n+2)-var(n+1)
      var(2)=2.d0*var(3)-var(4)
      var(1)=2.d0*var(2)-var(3)
    
      do i = 1, n
      wi_1=dabs(var(i+3)-var(i+2))
      wi=dabs(var(i+1)-var(i))
      if ((wi_1+wi).eq.0.d0) then
          z(i)=(var(i+2)+var(i+1))/2.d0
      else
          z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi)
      end if
      end do
      
      do i=1, n-1
          dx=x(i+1)-x(i)
          a=(z(i+1)-z(i))*dx                      ! coeff intermediaires pour calcul wd
          b=y(i+1)-y(i)-z(i)*dx                   ! coeff intermediaires pour calcul wd
          co(1,i)=y(i)
          co(2,i)=z(i)
          !co(3,i)=-(a-3.*b)/dx**2                ! méthode wd
          !co(4,i)=(a-2.*b)/dx**3                 ! méthode wd
          co(3,i)=(3.d0*var(i+2)-2.d0*z(i)-z(i+1))/dx   ! méthode JP Moreau
          co(4,i)=(z(i)+z(i+1)-2.d0*var(i+2))/dx**2  !
          ! les coefficients donnés par imsl sont co(3,i)*2 et co(4,i)*6
          ! etrangement la fonction csval corrige et donne la bonne valeur ...
      end do
      co(1,n)=y(n)
      co(2,n)=z(n)
      co(3,n)=0.d0
      co(4,n)=0.d0
  
      deallocate(z)
      deallocate(var)
  
  end subroutine
  
  !
  ! Single precision spline evaluation
  !
  function fvn_s_spline_eval(x,n,br,co)
      implicit none
      real, intent(in) :: x           ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated
      integer, intent(in) :: n        ! number of intervals
      real, intent(in) :: br(n+1)     ! breakpoints
      real, intent(in) :: co(4,n+1)   ! spline coeeficients
      real :: fvn_s_spline_eval
      
      integer :: i
      real :: dx
      
      if (x<=br(1)) then
          i=1
      else if (x>=br(n+1)) then
          i=n
      else
      i=1
      do while(x>=br(i))
          i=i+1
      end do
      i=i-1
      end if
      dx=x-br(i)
      fvn_s_spline_eval=co(1,i)+co(2,i)*dx+co(3,i)*dx**2+co(4,i)*dx**3
  
  end function
  
  ! Double precision spline evaluation
  function fvn_d_spline_eval(x,n,br,co)
      implicit none
      double precision, intent(in) :: x           ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated
      integer, intent(in) :: n        ! number of intervals
      double precision, intent(in) :: br(n+1)     ! breakpoints
      double precision, intent(in) :: co(4,n+1)   ! spline coeeficients
      double precision :: fvn_d_spline_eval
      
      integer :: i
      double precision :: dx
      
      
      if (x<=br(1)) then
          i=1
      else if (x>=br(n+1)) then
          i=n
      else
      i=1
      do while(x>=br(i))
          i=i+1
      end do
      i=i-1
      end if
      
      dx=x-br(i)
      fvn_d_spline_eval=co(1,i)+co(2,i)*dx+co(3,i)*dx**2+co(4,i)*dx**3
  
  end function
  
  
  end module fvn_interpol