CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DOUBLE PRECISION FUNCTION SPLINE(X,Y,NPOINTS,XIN)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccC
CCC   INPUT: X  x-values with X(i) < X(i+1)
CCC          Y  y-values
CCC          NPOINT   number of points
CCC          XIN  x-value for which y will be returned
CCC --------------------------------------------------------------
CCC   WARNING: will crash if X(i) = X(i+1); exeption not caught
CCC --------------------------------------------------------------
CCC   Function returns y value for a corresponding x value, 
CCC   based on cubic spline.
CCC   Will never oscillates or overshoot. No need to solve matrix.
CCC   Also calculate constants for cubic in case needed (for integration).
CCC --------------------------------------------------------------
CCC   adopted from VBasic routine 
CCC   Developer: C Kruger, Guildford, UK     
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccC
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Y(NPOINTS),X(NPOINTS)
      
      IF( NPOINTS.LE.0 )THEN
         SPLINE = 0.
         RETURN
      ELSEIF( NPOINTS.EQ.1 )THEN
         SPLINE = Y(1)
         RETURN
      ENDIF
      
C1) Find LineNumber or segment/leg. Linear extrapolate if outside range.
      IF (XIN .GE. X(NPOINTS)) THEN
         B=( Y(NPOINTS)-Y(NPOINTS-1) )/( X(NPOINTS)-X(NPOINTS-1) )
         A= Y(NPOINTS) - B*X(NPOINTS)

         SPLINE = A + B*XIN
*         SPLINE = Y(NPOINTS)
         RETURN
      ELSEIF (XIN .LE. X(1)) THEN
         B=( Y(2)-Y(1) )/( X(2)-X(1) )
         A= Y(2) - B*X(2)

         SPLINE = A + B*XIN
*         SPLINE= Y(1)
         RETURN      
      ELSEIF( NPOINTS.EQ.2 )THEN
         B=( Y(2)-Y(1) )/( X(2)-X(1) )
         A= Y(2) - B*X(2)

         SPLINE = A + B*XIN
         RETURN     
      ELSE
         LEG = NPOINTS
         DO I=2,NPOINTS
            IF (XIN .LE. X(I)) THEN
               LEG=I
               GOTO 10
            ENDIF
         ENDDO
10    CONTINUE
      ENDIF
      
C2) Calc first derivative (slope) for intermediate points
C3) Reset first derivative (slope) at first and last point
      IF (LEG .EQ. 2) THEN
*        First point has 0 2nd derivative  
         FPI   = FPRIME(X(LEG-1),Y(LEG-1))    
         FPIM1 = (3.*( Y(LEG)-Y(LEG-1) )/( X(LEG)-X(LEG-1) ) - FPI)/2.
      ELSEIF (LEG .EQ. NPOINTS) THEN
*        Last point has 0 2nd derivative
         FPIM1 = FPRIME(X(LEG-2),Y(LEG-2))         
         FPI = (3.*( Y(LEG)-Y(LEG-1) )/( X(LEG)-X(LEG-1) ) -  FPIM1)/2.
      ELSE
         FPIM1 = FPRIME(X(LEG-2),Y(LEG-2))
         FPI   = FPRIME(X(LEG-1),Y(LEG-1))
      ENDIF


C4) Calc second derivative at points
      F2IM1 = -2D0*( FPI + 2D0*FPIM1 )/( X(LEG)-X(LEG-1) )
     $       + 6D0*( Y(LEG)-Y(LEG-1) )/( X(LEG)-X(LEG-1) )**2
 
      F2I   =  2D0*( 2D0*FPI + FPIM1 )/( X(LEG)-X(LEG-1) )
     $       - 6D0*( Y(LEG)-Y(LEG-1) )/( X(LEG)-X(LEG-1) )**2


C5) Calc constants for cubic
      D = ( F2I-F2IM1 )/( 6D0*(X(LEG)-X(LEG-1)) ) 
      C = ( X(LEG)*F2IM1 - X(LEG-1)*F2I )/( 2D0*( X(LEG)-X(LEG-1 )) )
      B = ( Y(LEG)-Y(LEG-1) 
     $    - C*( X(LEG)**2-X(LEG-1)**2 ) - D*( X(LEG)**3-X(LEG-1)**3 ) ) 
     $    / ( X(LEG)-X(LEG-1) )
      A = Y(LEG-1) - B*X(LEG-1) - C*X(LEG-1)**2 - D*X(LEG-1)**3

      SPLINE = A + B*XIN + C*XIN**2 + D*XIN**3
      END
      DOUBLE PRECISION FUNCTION FPRIME(X,Y)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(3),Y(3)          
      
      DX1 = X(2)-X(1)
      DX2 = X(3)-X(2)
      DY1 = Y(2)-Y(1)
      DY2 = Y(3)-Y(2)
      
      
      IF( DY1.EQ.0. .OR. DY2.EQ.0. )THEN
* Prevent divide by 0
         FPRIME = 0.
      ELSEIF( DY1*DY2 .LT.0. )THEN
* Pos AND neg slope, assume slope = 0 to prevent overshoot
         FPRIME = 0.         
      ELSE
         FPRIME =  2D0/( DX1/DY1 + DX2/DY2 )
      ENDIF
      END