      SUBROUTINE SPLINE
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'mopac.inc'
      LOGICAL SKIP1,SKIP2
C
C     FIT F(X) BY A CUBIC SPLINE GIVEN VALUES OF THE FUNCTION
C     AND ITS FIRST DERIVATIVE AT N PNTS.
C     SUBROUTINE RETURNS VALUES OF XMIN,FMIN, AND DFMIN
C     AND MAY REORDER THE DATA.
C     CALLING PROGRAM SUPPLIES ALL OTHER VALUES IN THE
C     COMMON BLOCK.
C     XLOW AND XHIGH SET LIMITS ON THE INTERVAL WITHIN WHICH
C     TO SEARCH.  SUBROUTINE MAY FURTHER REDUCE THIS INTERVAL.
C
      COMMON/FIT/N,IDUM2,XLOW,XHIGH,XMIN,FMIN,DFMIN,X(12),F(12),DF(12)
      DATA CLOSE, BIG, HUGE, USTEP, DSTEP/1.0E-8,500.0,1.0E+10,1.0,2.0/
C
C     SUBROUTINE ASSUMES THAT THE FIRST N-1 DATA PNTS HAVE BEEN
C     PREVIOUSLY ORDERED,  X(I).LT.X(I+1) FOR I=1,2,...,N-2
C     NOW MOVE NTH POINT TO ITS PROPER PLACE.
C
      XMIN=X(N)
      FMIN=F(N)
      DFMIN=DF(N)
      N1=N-1
      K=N1
   10 IF(X(K).LT.XMIN) GO TO 20
      X(K+1)=X(K)
      F(K+1)=F(K)
      DF(K+1)=DF(K)
      K=K-1
      IF(K.GT.0) GO TO 10
   20 X(K+1)=XMIN
      F(K+1)=FMIN
      DF(K+1)=DFMIN
C
C     DEFINE THE INTERVAL WITHIN WHICH WE TRUST THE SPLINE FIT.
C     USTEP =  UP HILL STEP SIZE FACTOR
C     DSTEP = DOWN HILL STEP SIZE FACTOR
C
      IF(DF(1).GT.0.0) STEP=DSTEP
      IF(DF(1).LE.0.0) STEP=USTEP
      XSTART=X(1)-STEP*(X(2)-X(1))
      XSTART=MAX(XSTART,XLOW)
      IF(DF(N).GT.0.0) STEP=USTEP
      IF(DF(N).LE.0.0) STEP=DSTEP
      XSTOP=X(N)+STEP*(X(N)-X(N1))
      XSTOP=MIN(XSTOP,XHIGH)
C
C     SEARCH FOR MINIMUM
C
      DO 110 K=1,N1
         SKIP1=K.NE.1
         SKIP2=K.NE.N1
         IF(F(K).GE.FMIN) GO TO 30
         XMIN=X(K)
         FMIN=F(K)
         DFMIN=DF(K)
   30    DX=X(K+1)-X(K)
C
C     SKIP INTERVAL IF PNTS ARE TOO CLOSE TOGETHER
C
         IF(DX.LE.CLOSE) GO TO 110
         X1=0.0
         IF(K.EQ.1) X1=XSTART-X(1)
         X2=DX
         IF(K.EQ.N1) X2=XSTOP-X(N1)
C
C     (A,B,C)=COEF OF (CUBIC,QUADRATIC,LINEAR) TERMS
C
         DUM=(F(K+1)-F(K))/DX
         A=(DF(K)+DF(K+1)-DUM-DUM)/(DX*DX)
         B=(DUM+DUM+DUM-DF(K)-DF(K)-DF(K+1))/DX
         C=DF(K)
C
C     XK = X-X(K) AT THE MINIMUM WITHIN THE KTH SUBINTERVAL
C     TEST FOR PATHOLOGICAL CASES.
C
         BB=B*B
         AC3=(A+A+A)*C
         IF(BB.LT.AC3) GO TO 90
         IF( B.GT.0.0) GO TO 40
         IF(ABS(B).GT.HUGE*ABS(A)) GO TO 90
         GO TO 50
   40    IF(BB.GT.BIG*ABS(AC3)) GO TO 60
C
C     WELL BEHAVED CUBIC
C
   50    XK=(-B+SQRT(BB-AC3))/(A+A+A)
         GO TO 70
C
C     CUBIC IS DOMINATED BY QUADRATIC TERM
C
   60    R=AC3/BB
         XK=-(((0.039063*R+0.0625)*R+0.125)*R+0.5)*C/B
   70    IF(XK.LT.X1.OR.XK.GT.X2) GO TO 90
   80    FM=((A*XK+B)*XK+C)*XK+F(K)
         IF(FM.GT.FMIN) GO TO 90
         XMIN=XK+X(K)
         FMIN=FM
         DFMIN=((A+A+A)*XK+B+B)*XK+C
C
C     EXTRAPOLATE TO END OF INTERVAL IF K=1 AND/OR K=N1
C
   90    IF(SKIP1) GO TO 100
         SKIP1=.TRUE.
         XK=X1
         GO TO 80
  100    IF(SKIP2) GO TO 110
         SKIP2=.TRUE.
         XK=X2
         GO TO 80
  110 CONTINUE
      RETURN
      END
