SUBROUTINE spline(x,y,n,yp1,ypn,y2) 
IMPLICIT NONE
INTEGER(4),PARAMETER :: NMAX=500
INTEGER(4) :: n,i,k
REAL(8) :: yp1,ypn,x(n),y(n),y2(n),p,qn,sig,un,u(NMAX)
!Given arrays x(1:n) and y(1:n) containing a tabulated function, i.e., 
!yi = f(xi), with x1 < x2 < .. . < xN, and given values yp1 and ypn for 
!the first derivative of the interpolating function at points 1 and n, 
!respectively, this routine returns an array y2(1:n) of length n which 
!contains the second derivatives of the interpolating function at the 
!tabulated points xi. If yp1 and/or ypn are equal to 1 × 1030 or larger, 
!the routine is signaled to set the corresponding boundary condition for 
!a natural spline, with zero second derivative on that boundary. 
!Parameter: NMAX is the largest anticipated value of n. 
if (yp1>.99e30) then 
   y2(1)=0.  ! The lower boundary condition is set either to be  natural
   u(1)=0. 
else !or else to have a specified first derivative. 
   y2(1)=-0.5 
   u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1) 
endif 

do i=2,n-1 
!This is the decomposition loop of the tridiagonal algorithm. y2 and u 
!are used for temporary storage of the decomposed factors. 
   sig=(x(i)-x(i-1))/(x(i+1)-x(i-1)) 
   p=sig*y2(i-1)+2. 
   y2(i)=(sig-1.)/p 
   u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) &
         /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p 
enddo

if (ypn>.99e30) then 
   qn=0. !The upper boundary condition is set either to be  natural
   un=0. 
else 
   qn=0.5 
   un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) 
endif 

y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.) 

do k=n-1,1,-1 
   y2(k)=y2(k)*y2(k+1)+u(k) 
enddo

return 
END