SUBROUTINE splint(xa,ya,y2a,n,x,y) 
IMPLICIT NONE

INTEGER(4) :: n,k,khi,klo
REAL(8) :: x,y,xa(n),y2a(n),ya(n), a,b,h

!Given the arrays xa(1:n) and ya(1:n) of length n, which tabulate xai s 
!in order), and given the array y2a(1:n), which is the output and given 
!a value of x, this routine returns a cubic-spline interpolated 
! 
!We will find the right place in the table by This is optimal if sequential 
! calls to this routine values of x. If sequential calls are in spaced, one 
!would do better to store klo and khi and test if they remain next call. 

klo=1
khi=n

1 continue
 
if ((khi-klo)>1) then 
   k=(khi+klo)/2 
   if(xa(k)>x)then 
      khi=k 
   else 
      klo=k 
   endif 
   goto 1 
endif
 
!klo and khi now bracket the input value 

h=xa(khi)-xa(klo) 

if (h.eq.0.) pause  'bad xa input in splint' 
!he xa s must be distinct. 
a=(xa(khi)-x)/h 
!Cubic spline polynomial is now evaluated. 
b=(x-xa(klo))/h 
y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6. 

return 
END