SUBROUTINE PCS( X, Y, N, X0, FX0, FPX0, INTRVL, WORK )
C
C	This subroutine will interpolate a set of points using a
C	piecewise-cubic spline.
C
C	X	Independent array.
C	Y	Dependent array.
C	N	Number of points in the above arrays.
C	X0	Value at which to calculate the interpolate.
C	FX0	Interpolated value.
C	FPX0	Derivative of the interpolate at X0.
C	INTRVL	Equals interval number of where X0 is.
C		Set to 0 or N if extrapolation, must equal -1
C		on first call to initialize interpolation.
C	WORK	Working array, dimensioned at least 3*N.
C
	DIMENSION X(1), Y(1), WORK(1)
C
C - INTERPOLATION ALREADY CALCULATED.
	IF( INTRVL .GE. 0 ) GO TO 5100
C
C - CONSTABLES.
	N2 = 0
	N3 = N
	N4 = N + N3
	ND = N3
	NS = N4
C
C - CALCULATE SLOPES FOR END POINTS.
	CALL ORTHO( X, Y, 4, WORK(1), 3, WORK(5), WORK(11) )
	S1 = WORK(2) + X(1)*( 2.*WORK(3) + 
	1  3.*X(1)*WORK(4) )
C
	CALL ORTHO( X( N-3 ), Y( N-3 ), 4, WORK(1), 3, WORK(5),
	1  WORK(11) )
	WORK( N2+N ) = WORK(2) + X(N)*( 2.*WORK(3) +
	1  3.*X(N)*WORK(4) )
	WORK( N2+1 ) = S1
	WORK( ND+1 ) = 0.
	WORK( NS+1 ) = 1.
C
C - CALCULATE DIFFERENCES AND SLOPES.
	DO 1000 I=2,N
	WORK( ND+I ) = X( I ) - X( I-1 )
1000	WORK( NS+I ) = ( Y( I ) - Y( I-1 ) )/WORK( ND+I )
C
	DO 2000 I=2,N-1
	WORK( N2+I ) = 3.*( WORK( ND+I)*WORK( NS+I+1 ) +
	1  WORK( ND+I+1 )*WORK( NS+I ) )
2000	WORK( NS+I ) = 2.*( WORK( ND+I ) + WORK( ND+I+1 ) )
C
	DO 3000 I=2,N-1
	G = -WORK( ND+I+1)/WORK( NS+I-1 )
	WORK( NS+I ) = WORK( NS+I ) + G*WORK( ND+I-1 )
3000	WORK( N2+I ) = WORK( N2+I ) + G*WORK( N2+I-1 )
C
	DO 4000 I=N-1,1,-1
4000	WORK( N2+I ) = ( WORK( N2+I ) - WORK( ND+I )*WORK( N2+I+1 ) )/
	1  WORK( NS+I )
C
	DO 5000 I=1,N-1
	DX = X( I+1 ) - X( I )
	F1 = ( Y( I+1 ) - Y( I ) )/DX
	F3 = WORK( N2+I ) + WORK( N2+I+1 ) - 2.*F1
	WORK( N3+I ) = ( F1 - WORK( N2+I ) - F3 )/DX
5000	WORK( N4+I ) = F3/( DX*DX )
	INTRVL = 1
C
C - FIND WHICH INTERVAL POINT IS WITHIN.
5100	IF( X0 .GT. X(1) ) GO TO 5200
	INTRVL = 0
	ISEG = 1
	GO TO 5500
5200	IF( X0 .LT. X(N) ) GO TO 5300
	INTRVL = N
	ISEG = N-1
	GO TO 5500
C
C - MOVE TO PROPER SEGMENT.
5300	INTRVL = MIN0( N-1, MAX0( 1, INTRVL ) )
	IF( ( X(INTRVL) .LE. X0 ) .AND. ( X0 .LE. X(INTRVL+1) ) )
	1  GO TO 5400
	INTRVL = INTRVL + SIGN( 1., X0 - X(INTRVL) )
	GO TO 5300
C
C - CALCULATE INTEROPLATE AND DERIVATIVE.
5400	ISEG = INTRVL
5500	DX = X0 - X( ISEG )
	FX0 = Y( ISEG ) + DX*( WORK( N2+ISEG ) +
	1  DX*( WORK( N3+ISEG ) + DX*WORK( N4+ISEG ) ) )
	FPX0 = WORK( N2+ISEG ) + DX*( 2.*WORK( N3+ISEG ) +
	1  DX*3.*WORK( N4+ISEG ) )
C
	RETURN
	END