function spline(u,xx)

!
!  the function spline returns a cubic spline interpolated value of a 
!   sequence of four real numbers u, interpolated between the center two 
!   values.  The interpolation point is defined by xx, which ranges from 0 to 1
!   corresponding to the samples from u(2) to u(3).  As an example, the call
!
!       y=spline(u,0.25)
!
!   returns the cubic spline estimate of the value of the function defined by
!   u(1) u(2) u(3) u(4) interpolated to a point 0.25 of the way from u(2) to
!   u(3).  The u(*) are any set of real values assumed to be evenly spaced
!   along the x axis.  u(x) are the data points measured at each location x.
!

real x(4),u(4),a(4),s(4)

data x/-1.,0.,1.,2./

q(u1,x1,u2,x2)=(u1/x1**2-u2/x2**2)/(1.0/x1-1.0/x2)

q1=q(u(2)-u(1),x(2)-x(1),u(3)-u(1),x(3)-x(1))
qn=q(u(3)-u(4),x(3)-x(4),u(2)-u(4),x(2)-x(4))

s(1)=6.0*((u(2)-u(1))/(x(2)-x(1)) - q1)

do i=2,3
   s(i)= (u(i-1)/(x(i)-x(i-1)) - u(i)*(1.0/(x(i)-x(i-1))+ 1.0/(x(i+1)-x(i))) + u(i+1)/(x(i+1)-x(i)))*6.0
enddo
s(4)=6.0*(qn + (u(3)-u(4))/(x(4)-x(3)))
a(1)=2.0*(x(2)-x(1))
a(2)=1.5*(x(2)-x(1)) + 2.0*(x(3)-x(2))
s(2)=s(2) - 0.5*s(1)
do i=3,3
   c=(x(i)-x(i-1))/a(i-1)
   a(i)=2.0*(x(i+1)-x(i-1)) - c*(x(i)-x(i-1))
   s(i)=s(i) - c*s(i-1)
enddo
c=(x(4)-x(3))/a(3)
a(4)=(2.0-c)*(x(4)-x(3))
s(4)=s(4) - c*s(3)
!  back substitiute
s(4)= s(4)/a(4)
do j=1,3
   i=4-j
   s(i) =(s(i) - (x(i+1)-x(i))*s(i+1))/a(i)
enddo

dt = xx - x(2)
h = x(3)-x(2)

spline = u(2) &
     + dt * ( ( u(3)-u(2)) / h &
     - ( s(3) / 6.0E+00 + s(2) / 3.0E+00 ) * h &
     + dt * ( 0.5E+00 * s(2) &
     + dt * ( ( s(3)-s(2) ) / ( 6.0E+00 * h ) ) ) )
           
return

end function spline
