SUBROUTINE SPLINE(X,Y,N,YP1,YPN,Y2,U)
C
C     GIVEN N TABULATED FUNCTIONS X, Y AND FIRST DERIVATIVES
C     AT FIRST AND LAST POINTS YP1, YPN, THIS SUBROUTINE RETURNS
C     Y2, THE 2ND DERIVATIVES AT N POINTS. IF YP1 AND/OR YPN ARE
C     UNKNOWN, SET THEM TO 10**30 OR LARGER AND THEY ARE CALCULATED
C     ON RETURN. U IS A SCRATCH ARRAY OF DIMENSION N.
C----------------------------------------------------------------
      IMPLICIT NONE
      INTEGER N,I,K,nmax
      LOGICAL ORDER
      parameter (nmax=1000)
      double precision X(nmax),Y(nmax),U(nmax),Y2(nmax),YP1,YPN,
     &       P,ONE,TWO,SIX,XB,XF,YB,YF
c     REAL*8 X(N),Y(N),U(N),Y2(N),YP1,YPN,SIG,P,
c    &       ONE,TWO,SIX,XB,XF,YB,YF
      PARAMETER (ONE=1.0,TWO=2.0,SIX=6.0)
      IF (N.LT.4) THEN
         WRITE(6,*)'ERROR: SPLINE FITTING DATA MUST > 3'
         STOP
      END IF
      ORDER=.TRUE.
C     === REORDER DATA IF X NOT IN ASCENDING ORDER
      DO 5 I=2,N
         IF (X(I).EQ.X(I-1)) THEN
            WRITE(6,*)'ERROR: IDENTICAL DATA FOR PAIR ',
     &                I-1,I
            STOP
         ELSE IF (X(I).LT.X(I-1)) THEN
            ORDER=.FALSE.
         END IF
5     CONTINUE
      IF (.NOT. ORDER) THEN
         WRITE(6,*)'*** WARNING: DATA ARE NOT IN ASCENDING ORDER.',
     &             ' REORDER THE DATA ...'
         CALL SORT(X,Y,N)
      END IF
C     === CAL. SPLINE COEFFCIENTS Y2
      XF=X(3)-X(2)
      XB=X(2)-X(1)
      YF=Y(3)-Y(2)
      YB=Y(2)-Y(1)
      IF (YP1.GT..99E30) THEN
        P=TWO*XF+XB
        Y2(2)=(XB-XF)/P
        U(2)=SIX*(YF-YB*XF/XB)/((XF+XB)*P)
      ELSE
        Y2(1)=-0.5
        U(1)=(3./XB)*(YB/XB-YP1)
        P=TWO*(XF+XB)+XB*Y2(1)
        Y2(2)=-XF/P
        U(2)=(SIX*(YF/XF-YB/XB)-XB*U(1))/P
      ENDIF
      DO 11 I=3,N-1
        XF=X(I+1)-X(I)
        XB=X(I)-X(I-1)
        YF=Y(I+1)-Y(I)
        YB=Y(I)-Y(I-1)
        P=TWO*(XF+XB)+XB*Y2(I-1)
        Y2(I)=-XF/P
        U(I)=(SIX*(YF/XF-YB/XB)-XB*U(I-1))/P
11    CONTINUE
      IF (YPN.GT..99E30) THEN
        P=TWO*XF+XB
        Y2(N)=(SIX*(YF-YB*XF/XB)/(XF+XB)-P*U(N-1))
     *        /(XF-XB+P*Y2(N-1))
      ELSE
        Y2(N)=(SIX*(YPN-YF/XF)/XF-U(N-1))/(TWO+Y2(N-1))
      ENDIF
      DO 12 K=N-1,2,-1
        Y2(K)=Y2(K)*Y2(K+1)+U(K)
12    CONTINUE
      IF (YP1.GT..99E30) THEN
        Y2(1)=Y2(2)-(X(2)-X(1))*(Y2(3)-Y2(2))/(X(3)-X(2))
        YP1=(Y(2)-Y(1))/(X(2)-X(1))
     &     -(X(2)-X(1))*(Y2(1)/3D0+Y2(2)/6D0)
      ELSE
        Y2(1)=Y2(1)*Y2(2)+U(1)
      END IF
      IF (YPN.GT..99E30) THEN
        YPN=YF/XF+XF*(Y2(N-1)/6D0+Y2(N)/3D0)
      END IF
         RETURN
            END


      SUBROUTINE SPLINT(XA,YA,Y2A,N,X,Y,KIND)
C---------------------------------------------------------------
C     CUBIC SPLINE FITTING PROGRAM
C
C     GIVEN N TABULATED FUNCTIONS XA, YA AND 2ND DERIVATIVES Y2A,
C     THIS SUBROUTINE RETURNS AT X
C              FUNCTION VALUE FOR KIND=0
C              1ST DERIVATIVE FOR KIND=1
C              2nd DERIVATIVE FOR KIND=2
C----------------------------------------------------------------
      IMPLICIT NONE
      INTEGER N,K,KHI,KLO,KIND,nmax
      parameter (nmax=1000)
      double precision XA(nmax),YA(nmax),Y2A(nmax),A,B,H,X,y
c     REAL*8 XA(N),YA(N),Y2A(N),A,B,H,X,Y
      KLO=1
      KHI=N
1     IF (KHI-KLO.GT.1) THEN
        K=(KHI+KLO)/2
        IF(XA(K).GT.X)THEN
          KHI=K
        ELSE
          KLO=K
        ENDIF
      GOTO 1
      ENDIF
      H=XA(KHI)-XA(KLO)
      IF (H.EQ.0.) PAUSE 'Bad XA input.'
      A=(XA(KHI)-X)/H
      B=(X-XA(KLO))/H
      IF (KIND.EQ.0) THEN
         Y=A*YA(KLO)+B*YA(KHI)+
     *      ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6.
      ELSE IF (KIND.EQ.1) THEN
         Y=(YA(KHI)-YA(KLO))/H+
     *      ((1.-3.*A*A)*Y2A(KLO)+(3.*B*B-1.)*Y2A(KHI))*H/6.
      ELSE
         Y= A*Y2A(KLO)+B*Y2A(KHI)
      END IF
      RETURN
      END