module simplespline

  use const_and_precisions, only : wp_
  implicit none

contains

function spli(cspli,n,k,dx)
  implicit none
  integer, intent(in) :: n, k
  real(wp_), intent(in) :: cspli(n,4), dx
  real(wp_) :: spli
  spli=cspli(k,1)+dx*(cspli(k,2)+dx*(cspli(k,3)+dx*cspli(k,4)))
end function spli

function splid(cspli,n,k,dx)
  implicit none
  integer, intent(in) :: n, k
  real(wp_), intent(in) :: cspli(n,4), dx
  real(wp_) :: splid
  splid=cspli(k,2)+dx*(2.0_wp_*cspli(k,3)+3.0_wp_*dx*cspli(k,4))
end function splid

subroutine difcs(x,y,n,iopt,c,ier)
  implicit none
  integer, intent(in) :: n, iopt
  real(wp_), intent(in) :: x(n), y(n)
  real(wp_), intent(inout) :: c(n*4)
  integer :: ier
  integer :: jmp,iol,ioh,i,ii,j,j1,j2,j3
  real(wp_) :: xb,xc,ya,yb,h,a,r,dya,dyb,dy2
  jmp  =1
  if (n <= 1) return
!
! initialization
!
  xc   =x(1)
  yb   =y(1)
  h    =0.0_wp_
  a    =0.0_wp_
  r    =0.0_wp_
  dyb  =0.0_wp_
!
! iol=0 - given derivative at first point
! ioh=0 - given derivative at last point
!
  iol  =iopt-1
  ioh  =iopt-2
  if (ioh == 1) then
    iol  =0
    ioh  =0
  end if
  dy2  =c(2)
!
! form the system of linear equations
! and eliminate subsequentially
!
  j    =1
  do i=1,n
    j2   =n+i
    j3   =j2+n
    a    =h*(2.0_wp_-a)
    dya  =dyb+h*r
    if (i>=n) then
!
! set derivative dy2 at last point
!
      dyb  =dy2
      h    =0.0_wp_
      if (ioh/=0) then
        dyb  =dya
        goto 13
      end if
    else
      j    =j+jmp
      xb   =xc
      xc   =x(j)
      h    =xc-xb
!
! ii=0 - increasing abscissae
! ii=1 - decreasing abscissae
!
      ii   =0
      if (h==0) return
      if (h<0) ii   =1
      ya   =yb
      yb   =y(j)
      dyb  =(yb-ya)/h
      if (i<=1) then
        j1   =ii
        if (iol/=0) goto 13
        dya  =c(1)
      end if
    end if
    if (j1-ii /= 0) return
    a    =1.0_wp_/(h+h+a)
    13 continue
    r    =a*(dyb-dya)
    c(j3)=r
    a    =h*a
    c(j2)=a
    c(i) =dyb
  end do
!
! back substitution of the system of linear equations
! and computation of the other coefficients
!
  a    =1.0_wp_
  j1   =j3+n+ii-ii*n
  i    =n
  do iol=1,n
    xb   =x(j)
    h    =xc-xb
    xc   =xb
    a    =a+h
    yb   =r
    r    =c(j3)-r*c(j2)
    ya   =r+r
    c(j3)=ya+r
    c(j2)=c(i)-h*(ya+yb)
    c(j1)=(yb-r)/a
    c(i) =y(j)
    a    =0.0_wp_
    j    =j-jmp
    i    =i-1
    j2   =j2-1
    j3   =j3-1
    j1   =j3+n+ii
  end do
  ier  =0
end subroutine difcs

subroutine difcsn(xx,yy,nmx,n,iopt,cc,ier)
!
! same as difcs but with dimension(xx,yy) = nmx > n
!
  implicit none
  integer, intent(in) :: nmx, n, iopt
  real(wp_), intent(in) :: xx(nmx), yy(nmx)
  real(wp_), intent(inout) :: cc(nmx,4)
  integer :: ier
  integer :: jmp,iol,ioh,i,ii,j,j1,j2,j3
  real(wp_) :: x(n),y(n),c(n*4),xb,xc,ya,yb,h,a,r,dya,dyb,dy2
!
  do i=1,n
    x(i)=xx(i)
    y(i)=yy(i)
  end do
  ii=0
  do j=1,4
    do i=1,n
      ii=ii+1
      c(ii)=cc(i,j)
    end do
  end do
!
  jmp  =1
  if (n>1) then
!
! initialization
!
    xc   =x(1)
    yb   =y(1)
    h    =0.0_wp_
    a    =0.0_wp_
    r    =0.0_wp_
    dyb  =0.0_wp_
!
! iol=0 - given derivative at first point
! ioh=0 - given derivative at last point
!
    iol  =iopt-1
    ioh  =iopt-2
    if (ioh==1) then
      iol  =0
      ioh  =0
    end if
    dy2  =c(2)
!
! form the system of linear equations
! and eliminate subsequentially
!
    j    =1
    do i=1,n
      j2   =n+i
      j3   =j2+n
      a    =h*(2.0_wp_-a)
      dya  =dyb+h*r
      if (i>=n) then
!
! set derivative dy2 at last point
!
        dyb  =dy2
        h    =0.0_wp_
        if (ioh/=0) then
          dyb  =dya
          goto 13
        end if
      else
        j    =j+jmp
        xb   =xc
        xc   =x(j)
        h    =xc-xb
!
! ii=0 - increasing abscissae
! ii=1 - decreasing abscissae
!
        ii   =0
        if (h==0) goto 16
        if (h<0) ii   =1
        ya   =yb
        yb   =y(j)
        dyb  =(yb-ya)/h
        if (i<=1) then
          j1   =ii
          if (iol/=0) goto 13
          dya  =c(1)
        end if
      end if
      if (j1/=ii) goto 16
      a    =1.0_wp_/(h+h+a)
      13 continue
      r    =a*(dyb-dya)
      c(j3)=r
      a    =h*a
      c(j2)=a
      c(i) =dyb
    end do
!
! back substitution of the system of linear equations
! and computation of the other coefficients
!
    a    =1.0_wp_
    j1   =j3+n+ii-ii*n
    i    =n
    do iol=1,n
       xb   =x(j)
       h    =xc-xb
       xc   =xb
       a    =a+h
       yb   =r
       r    =c(j3)-r*c(j2)
       ya   =r+r
       c(j3)=ya+r
       c(j2)=c(i)-h*(ya+yb)
       c(j1)=(yb-r)/a
       c(i) =y(j)
       a    =0.0_wp_
       j    =j-jmp
       i    =i-1
       j2   =j2-1
       j3   =j3-1
       j1   =j3+n+ii
    end do
    ier  =0
  end if
!
  16 continue
  ii=0
  do j=1,4
    do i=1,nmx
      if(i<=n) then
        ii=ii+1
        cc(i,j)=c(ii)
       else
         cc(i,j)=0.0_wp_
       end if
    end do
  end do
!
end subroutine difcsn

end module simplespline