      SUBROUTINE SPLINE
      implicit double precision (A-H,O-Z)
      include 'SIZES'
      include 'MOP6FILE'
      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)
      save CLOSE, BIG, HUGE, USTEP, DSTEP
      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
