CALL PCS( X, Y, N, X0, FX0, FPX0, INTRVL, WORK ) --  This subroutine
 will interpolate a set of points using a piecewise-cubic spline.

 X - Independent array.
 Y - Dependent array.
 N - Number of points in the above arrays.
 X0 - Value at which to calculate the interpolate.
 FX0 - Interpolated value.
 FPX0 - Derivative of the interpolate at X0.
 INTRVL - Equals interval number of where X0 is.
          Set to 0 or N if extrapolation, must equal -1
          on first call to initialize interpolation.
 WORK - Working array, dimensioned at least 3*N.