c-----------------------------------------------
      subroutine splinev(n,x,y,m,t,yp1,ypn,sy)
C     PROGRAM FOR COMPLETE OR NATURAL CUBIC SPLINE INTERPOLATION
C     SPLINEV EVALUATES CUBIC SPLINE INTERPOLATION AT POINTS T(M).
C     INPUT: N,X(N),Y(N),M,T(M),YP1,YPN
C       N: NUMBER OF INTERPOLATION KNOTS
C       X: ARRAY OF INTERPOLATION KNOTS. X MUST BE AN ORDERED ARRAY, 
C          I.E., X(1)<X(2)<...<X(N).
C       Y: FUNCTION VALUES AT KNOTS IN X
C       M: NUMBER OF KNOTS AT WHICH VALUES OF SPLINE ARE FOUND 
C       T: ARRAY OF POINTS AT WHICH VALUES OF SPLINE ARE FOUND
C       YP1: VALUE OF FIRST DERIVATVE AT ENDPOINT X(1)
C       YPN: VALUE OF FIRST ORDER DERIVATVE AT ENDPOINT X(N)
C         NOTE:(1)THE CUBIC SPLINE FUNCTION IS CALLED COMPLETELY CUBIC 
C              SPLINE FUNCTION IF S'(X(1))=YP1 AND S'(X(N))=YPN.
C          (2)THE SUBROUTINE USES THE NATURAL SPLINE FUNCTION WHEN 
C             BOTH YP1 AND YPN ARE GIVEN NUMBERS > OR = 1.E30. 
C             THAT IS S''(X(1))=S''(X(N))=0.
C     OUTPUT: SY(M)
C       SY: VALUE OF SPLINE AT T.
C     WORK ARRAY: Y2(N)
C       Y2: VALUES OF SECOND DERITAVE OF SPLINE FUNCTION AT KNOTS IN X
C     Modified By Dr. Jianping Li on October 17, 1998. 
      dimension x(n),y(n),t(m),sy(m)
      dimension y2(n)
      call spline01(n,x,y,yp1,ypn,y2)
	do 100 i=1,m
      klo=1
      khi=n
   1  if((khi-klo).gt.1)then
        k=(khi+klo)/2
        if(x(k).gt.t(i))then
          khi=k
        else
          klo=k
        endif
        goto 1
      endif
      h=x(khi)-x(klo)
      if(h.eq.0.)pause'BAD X INPUT'
      a=(x(khi)-t(i))/h
      b=(t(i)-x(klo))/h
      sy(i)=a*y(klo)+b*y(khi)+((a**3-a)*y2(klo)
     *      +(b**3-b)*y2(khi))*(h**2)/6.
 100  continue
      return
      end
c----------------------------------------------------
      subroutine spline01(n,x,y,yp1,ypn,y2)
C     SPLINE01 DETERMINES CUBIC SPLINE INTEROPLATION
C     INPUT: N,X,Y,YP1,YPN 
C       N: NUMBER OF INTERPOLATION POINTS
C       X: ARRAY OF INTERPOLATION POINTS
C       Y: FUNCTION VALUES AT KNOTS IN X
C       YP1: VALUE OF FIRST DERIVATVE AT ENDPOINT X(1)
C       YPN: VALUE OF FIRST ORDER DERIVATVE AT ENDPOINT X(N)
C         NOTE:(1)THE CUBIC SPLINE FUNCTION IS CALLED COMPLETELY CUBIC 
C              SPLINE FUNCTION IF S'(X(1))=YP1 AND S'(X(N))=YPN.
C          (2)THE SUBROUTINE USES THE NATURAL SPLINE FUNCTION WHEN 
C             BOTH YP1 AND YPN ARE GIVEN NUMBERS > OR = 1.E30. 
C             THAT IS S''(X(1))=S''(X(N))=0.
C     OUTPUT: Y2
C       Y2: ARRAY OF VALUES OF SECOND DERATIVE AT KNOTS X(I)
C     U: WORK ARRAY
C     Modified By Dr. Jianping Li on October 17, 1998. 
      dimension x(n),y(n),u(100000),y2(n)
      if(yp1.gt.0.99e30)then
        y2(1)=0.
        u(1)=0.
      else
        y2(1)=-0.5
        u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      endif
      do 10 i=2,n-1
        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 
  10  continue
      if(ypn.gt.0.99e30)then
        qn=0.
        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 20 k=n-1,1,-1
        y2(k)=y2(k)*y2(k+1)+u(k)
  20  continue
      return
      end