c_________________________________________________________________________
c**********************************************************************
*SBR SPLINE: FITS CUBIC SPLINE CURVE TO TIME SERIES *******************
c**********************************************************************
*
*  DERIVED FROM IMSL LIBRARY ROUTINES BY EDWARD R COOK,
*  LAMONT-DOHERTY GEOLOGICAL OBSERVATORY, AND MODIFIED BY
*  RICHARD L HOLMES, TREE-RING LABORATORY, UNIVERSITY OF ARIZONA.
*
*  MODIFIED BY R L HOLMES 10 JUN 1988
*  MODIFIED BY E R COOK TO FULL DOUBLE PRECISION 4 APR 1991
*
********************************************************************
*
* CALCULATES A CUBIC SPLINE FUNCTION TO MODEL A TIME SERIES.
* 'STIFFNESS' OF SPLINE IS SPECIFIED BY PARAMETER 'LS'
*
* PARAMETERS:
*
* N:    NUMBER OF VALUES IN TIME SERIES
* Y:    TIME SERIES ARRAY TO BE MODELED
* F:    ARRAY OF CUBIC SPLINE FUNCTION VALUES COMPUTED
* P:    VALUE COMPUTED IN THIS ROUTINE, USED TO CALCULATE SPLINE
* LS:   LENGTH (STIFFNESS) OF SPLINE TO BE USED TO MODEL SERIES
* A:    WORKING ARRAY OF DIMENSION (N,4)
*
********************************************************************
*
      SUBROUTINE spline(maxdat,N,Y,F,LS)
      integer maxdat,ls
c      IMPLICIT REAL*8 (A-H,O-Z)
      integer mxy
      PARAMETER (MXY=10000)
      real*8 A(MXY,4),C1(4),C2(3)
      real*8 pct,pi,v,arg,rn,d1,d2,sum
      REAL*4 Y(maxdat),F(maxdat)
      DATA C1 /1.D0,-4.D0,6.D0,-2.D0/
      DATA C2 /0.D0,.33333333333333D0,1.33333333333333D0/
      DATA PCT,PI / .50D0, 3.1415926535897935D0 /
      IF(N .LT. 4)THEN
      F(1)=-9998.
      RETURN
      ENDIF
      IF(N.GT.MXY)THEN
        WRITE(*,'(/T2,''*** SERIES TOO LONG FOR SPLINE SUBROUTINE ***''/
     *           T2,'' *** ENLARGE MXY PARAMETER IN SUBROUTINE ***''/)')
        STOP
      ENDIF
      NM2=N-2
      V=DBLE(LS)
      ARG=(2.D0*PI)/V
      P=(6.D0*(DCOS(ARG)-1.D0)**2.D0)/(DCOS(ARG)+2.D0)
      P1=((1.D0/(1.D0-PCT)-1.D0)*6.D0*(DCOS(PI*2.D0/V)-1.D0)**2)/
     +     (DCOS(PI*2.D0/V)+2.D0)
      P2=(6.D0*(DCOS(ARG)-1.D0)**2.D0)/((2.D0**0.5D0-1.D0)*
     *		(DCOS(ARG)+2.D0))
      DO 1 I=1,NM2
      DO 1 J=1,3
      A(I,J)=C1(J)+P*C2(J)
    1 A(I,4)=DBLE(Y(I))+C1(4)*DBLE(Y(I+1))+DBLE(Y(I+2))
      A(1,1)=C2(1)
      A(1,2)=C2(1)
      A(2,1)=C2(1)
      NC=2
* BEGIN LUDAPB
      RN=DBLE(1.D0/(DFLOAT(NM2)*16.D0))
      D1=1.D0
      D2=0.D0
      NCP1=NC+1
      IF(NC .NE. 0)THEN
* INITIALIZE ZERO ELEMENTS
      DO 5 I=1,NC
      DO 5 J=I,NC
      K=NCP1-J
    5 A(I,K)=0.D0
* 'I' IS ROW INDEX OF ELEMENT BEING COMPUTED
* 'J' IS COLUMN INDEX OF ELEMENT BEING COMPUTED
* 'L' IS ROW INDEX OF PREVIOUSLY COMPUTED VECTOR
*       BEING USED TO COMPUTE INNER PRODUCT
* 'M' IS COLUMN INDEX
      ENDIF
      DO 60 I=1,NM2
      IMNCP1=I-NCP1
      I1=MAX0(1,1-IMNCP1)
      DO 60 J=I1,NCP1
      L=IMNCP1+J
      I2=NCP1-J
      SUM=A(I,J)
      JM1=J-1
      IF(JM1 .GT. 0)THEN
      DO 25 K=1,JM1
      M=I2+K
   25 SUM=SUM-A(I,K)*A(L,M)
      ENDIF
      IF(J .EQ. NCP1)THEN
      IF(A(I,J)+SUM*RN .LE. A(I,J))THEN
      F(1)=-9999.
      WRITE(*,'(''>> Sbr SPLINE: Matrix not positive definite'')')
      RETURN
      ENDIF
      A(I,J)=1.D0/DSQRT(SUM)
* UPDATE DETERMINANT
      D1=D1*SUM
   35 IF(DABS(D1) .GT. 1.D0)THEN
      D1=D1*.0625D0
      D2=D2+4.D0
      GOTO 35
      ENDIF
   47 IF(DABS(D1) .LE. .0625D0)THEN
      D1=D1*16.D0
      D2=D2-4.D0
      GOTO 47
      ELSE
      GOTO 60
      ENDIF
      ENDIF
      A(I,J)=SUM*A(L,NCP1)
   60 CONTINUE
* END LUDAPB / BEGIN LUELPB
* SOLUTION LY = B
      NC1=NC+1
      IW=0
      L=0
      DO 15 I=1,NM2
      SUM=A(I,4)
      IF(NC .GT. 0)THEN
      IF(IW .NE. 0)THEN
      L=L+1
      IF(L .GT. NC) L=NC
      K=NC1-L
      KL=I-L
      DO 45 J=K,NC
      SUM=SUM-A(KL,4)*A(I,J)
   45 KL=KL+1
      ELSE
      IF(SUM .NE. 0.D0)IW=1
      ENDIF
      ENDIF
   15 A(I,4)=SUM*A(I,NC1)
* SOLUTION UX = Y
      A(NM2,4)=A(NM2,4)*A(NM2,NC1)
      N1=NM2+1
      DO 30 I=2,NM2
      K=N1-I
      SUM=A(K,4)
      IF(NC .GT. 0)THEN
      KL=K+1
      K1=MIN0(NM2,K+NC)
      L=1
      DO 50 J=KL,K1
      SUM=SUM-A(J,4)*A(J,NC1-L)
   50 L=L+1
      ENDIF
   30 A(K,4)=SUM*A(K,NC1)
* END LUELPB
      DO 2 I=3,NM2
    2 F(I)=A(I-2,4)+C1(4)*A(I-1,4)+A(I,4)
      F(1)=A(1,4)
      F(2)=C1(4)*A(1,4)+A(2,4)
      F(N-1)=A(NM2-1,4)+C1(4)*A(NM2,4)
      F(N)=A(NM2,4)
      DO 3 I=1,N
    3 F(I)=Y(I)-F(I)
      RETURN
      END
