        SUBROUTINE SPLINE(X,Y,N,YP1,YPN,Y2)
C
C CUBIC SPLINE ROUTINE FROM NUMERICAL RECIPES, THE ART OF SCIENTIFIC
C COMPUTING. PP 88
C
C INPUT PARAMETERS
C       X,Y - TABULAR FUNTION OF LENGTH N
C       N - NUMBER OF POINTS IN X AND Y (100 points max)
C       YP1 - 1ST DIREVATIVE OF SPLINE AT FIRST POINT
C       YPN - 1ST DIREVATIVE OF SPLINE AT LAST POINT
C               IF YP1 AND YPN ARE SET TO 1.0E30 OR LARGER, THAN
C               A NATURAL SPLINE (2ND DIREVATIVE AT END POINTS)
C               WILL BE USED.
C OUTPUT PARAMETERS
C       Y2 - TABLE OF SECOND DIREVATIVES OF THE INTERPOLATING FUNCTION
C               AT EACH POINT I=1 TO N
C
	IMPLICIT NONE
        DOUBLE PRECISION X,Y,Y2,U,YP1,SIG,P,QN,UN,YPN
        INTEGER I,N,K
        DIMENSION X(N),Y(N),Y2(N),U(100)
C
C THE LOWER BOUNDARY NATURAL?
C
        IF(YP1.GT..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
C
C DO DECOMPOSITION LOOP OF THE TRIDIAGONAL ALGORITHM
C
        DO 11 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
11      CONTINUE
C
C UPPER BOUNDARY CONDITION
C
        IF(YPN.GT..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.)
C
C BACKSUBSTITUTION LOOP OF THE TRIDIAGONAL ALGORITHM
C
        DO 12 K=N-1,1,-1
                Y2(K)=Y2(K)*Y2(K+1)+U(K)
12      CONTINUE
        RETURN
        END
