subroutine spline(x,y,n,yp1,ypn,y2,u)
c
c=======================================================================
c                                                                    ===
c  Given X, Y of length N containing a tabulated function,  Y=f(X),  ===
c  with the  X  in ascending order,  and given values  Yp1 and  Ypn  ===
c  for the first derivative of the interpolating function at points  ===
c  1 and N, respectively this routine returns an array Y2 of length  ===
c  N  which contains the  second  derivatives of the  interpolating  ===
c  function at the tabulated points X.  If Yp1 and/or Ypn are equal  ===
c  to  1.0E+30  or larger,  the routine  is  signalled  to  set the  ===
c  corresponding boundary condition for a natural spline, with zero  ===
c  second derivative on that boundary.                               ===
c                                                                    ===
c  Reference :                                                       ===
c                                                                    ===
c  Press, W.H, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,   ===
c        1986: Numerical Recipes, the art of scientific computing.   ===
c        Cambridge University Press.                                 ===
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <pconst.h>
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
      integer i,k,n
      FLOAT
     *      p,qn,sig,un,ypn,yp1
      FLOAT
     *      x(n),y(n),y2(n),u(n)
c
c=======================================================================
c  Begin excutable code.
c=======================================================================
c
c-----------------------------------------------------------------------
c  The lower boundary condition is set either to be "natural" or else
c  to have a specified first derivative.
c-----------------------------------------------------------------------
c
      if(yp1.gt.c9p9e29) then
        y2(1)=c0
        u(1)=c0
      else
        y2(1)=mp5
        u(1)=(c3/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      end if
c
c-----------------------------------------------------------------------
c  This is the decomposition loop of the tridiagonal algorithm. Y2 and
c  U are used for temporary storage of the decomposition factors.
c-----------------------------------------------------------------------
c
      do 10 i=2,n-1
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
        p=sig*y2(i-1)+c2
        y2(i)=(sig-c1)/p
        u(i)=(c6*((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
c
c-----------------------------------------------------------------------
c  The upper boundary condition is set either to be "natural" or else
c  to have a specified first derivative.
c-----------------------------------------------------------------------
c
      if(ypn.gt.c9p9e29) then
        qn=c0
        un=c0
      else
        qn=p5
        un=(c3/(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)+c1)
c
c-----------------------------------------------------------------------
c  This is the back-substitution loop of the tridiagonal algorithm.
c-----------------------------------------------------------------------
c
      do 20 k=n-1,1,-1
        y2(k)=y2(k)*y2(k+1)+u(k)
 20   continue
c
      return
      end