      SUBROUTINE CHFDV(X1,X2,F1,F2,D1,D2,NE,XE,FE,DE,NEXT,IERR)
C***BEGIN PROLOGUE  CHFDV
C***DATE WRITTEN   811019   (YYMMDD)
C***REVISION DATE  870707   (YYMMDD)
C***CATEGORY NO.  E3,H1
C***KEYWORDS  LIBRARY=SLATEC(PCHIP),
C             TYPE=SINGLE PRECISION(CHFDV-S DCHFDV-D),
C             CUBIC HERMITE DIFFERENTIATION,CUBIC HERMITE EVALUATION,
C             CUBIC POLYNOMIAL EVALUATION
C***AUTHOR  FRITSCH, F. N., (LLNL)
C             MATHEMATICS AND STATISTICS DIVISION
C             LAWRENCE LIVERMORE NATIONAL LABORATORY
C             P.O. BOX 808  (L-316)
C             LIVERMORE, CA  94550
C             FTS 532-4275, (415) 422-4275
C***PURPOSE  Evaluate a cubic polynomial given in Hermite form and its
C            first derivative at an array of points.  While designed for
C            use by PCHFD, it may be useful directly as an evaluator for
C            a piecewise cubic Hermite function in applications, such as
C            graphing, where the interval is known in advance.
C            If only function values are required, use CHFEV instead.
C***DESCRIPTION
C
C        CHFDV:  Cubic Hermite Function and Derivative Evaluator
C
C     Evaluates the cubic polynomial determined by function values
C     F1,F2 and derivatives D1,D2 on interval (X1,X2), together with
C     its first derivative, at the points  XE(J), J=1(1)NE.
C
C     If only function values are required, use CHFEV, instead.
C
C ----------------------------------------------------------------------
C
C  Calling sequence:
C
C        INTEGER  NE, NEXT(2), IERR
C        REAL  X1, X2, F1, F2, D1, D2, XE(NE), FE(NE), DE(NE)
C
C        CALL  CHFDV (X1,X2, F1,F2, D1,D2, NE, XE, FE, DE, NEXT, IERR)
C
C   Parameters:
C
C     X1,X2 -- (input) endpoints of interval of definition of cubic.
C           (Error return if  X1.EQ.X2 .)
C
C     F1,F2 -- (input) values of function at X1 and X2, respectively.
C
C     D1,D2 -- (input) values of derivative at X1 and X2, respectively.
C
C     NE -- (input) number of evaluation points.  (Error return if
C           NE.LT.1 .)
C
C     XE -- (input) real array of points at which the functions are to
C           be evaluated.  If any of the XE are outside the interval
C           [X1,X2], a warning error is returned in NEXT.
C
C     FE -- (output) real array of values of the cubic function defined
C           by  X1,X2, F1,F2, D1,D2  at the points  XE.
C
C     DE -- (output) real array of values of the first derivative of
C           the same function at the points  XE.
C
C     NEXT -- (output) integer array indicating number of extrapolation
C           points:
C            NEXT(1) = number of evaluation points to left of interval.
C            NEXT(2) = number of evaluation points to right of interval.
C
C     IERR -- (output) error flag.
C           Normal return:
C              IERR = 0  (no errors).
C           "Recoverable" errors:
C              IERR = -1  if NE.LT.1 .
C              IERR = -2  if X1.EQ.X2 .
C                (Output arrays have not been changed in either case.)
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  XERROR
C***END PROLOGUE  CHFDV
C
C ----------------------------------------------------------------------
C
C  Change record:
C     82-08-03   Minor cosmetic changes for release 1.
C
C ----------------------------------------------------------------------
C
C  Programming notes:
C
C     To produce a double precision version, simply:
C        a. Change CHFDV to DCHFDV wherever it occurs,
C        b. Change the real declaration to double precision,
C        c. Change the constant ZERO to double precision, and
C        d. Change the names of the Fortran functions:  AMAX1, AMIN1.
C
C  DECLARE ARGUMENTS.
C
      INTEGER  NE, NEXT(2), IERR
      REAL  X1, X2, F1, F2, D1, D2, XE(NE), FE(NE), DE(NE)
C
C  DECLARE LOCAL VARIABLES.
C
      INTEGER  I
      REAL  C2, C2T2, C3, C3T3, DEL1, DEL2, DELTA, H, X, XMI, XMA, ZERO
      DATA  ZERO /0./
C
C  VALIDITY-CHECK ARGUMENTS.
C
C***FIRST EXECUTABLE STATEMENT  CHFDV
      IF (NE .LT. 1)  GO TO 5001
      H = X2 - X1
      IF (H .EQ. ZERO)  GO TO 5002
C
C  INITIALIZE.
C
      IERR = 0
      NEXT(1) = 0
      NEXT(2) = 0
      XMI = AMIN1(ZERO, H)
      XMA = AMAX1(ZERO, H)
C
C  COMPUTE CUBIC COEFFICIENTS (EXPANDED ABOUT X1).
C
      DELTA = (F2 - F1)/H
      DEL1 = (D1 - DELTA)/H
      DEL2 = (D2 - DELTA)/H
C                                           (DELTA IS NO LONGER NEEDED.)
      C2 = -(DEL1+DEL1 + DEL2)
      C2T2 = C2 + C2
      C3 = (DEL1 + DEL2)/H
C                               (H, DEL1 AND DEL2 ARE NO LONGER NEEDED.)
      C3T3 = C3+C3+C3
C
C  EVALUATION LOOP.
C
      DO 500  I = 1, NE
         X = XE(I) - X1
         FE(I) = F1 + X*(D1 + X*(C2 + X*C3))
         DE(I) = D1 + X*(C2T2 + X*C3T3)
C          COUNT EXTRAPOLATION POINTS.
         IF ( X.LT.XMI )  NEXT(1) = NEXT(1) + 1
         IF ( X.GT.XMA )  NEXT(2) = NEXT(2) + 1
C        (NOTE REDUNDANCY--IF EITHER CONDITION IS TRUE, OTHER IS FALSE.)
  500 CONTINUE
C
C  NORMAL RETURN.
C
      RETURN
C
C  ERROR RETURNS.
C
 5001 CONTINUE
C     NE.LT.1 RETURN.
      IERR = -1
      CALL XERROR ('CHFDV -- NUMBER OF EVALUATION POINTS LESS THAN ONE'
     *           , 50, IERR, 1)
      RETURN
C
 5002 CONTINUE
C     X1.EQ.X2 RETURN.
      IERR = -2
      CALL XERROR ('CHFDV -- INTERVAL ENDPOINTS EQUAL'
     *           , 33, IERR, 1)
      RETURN
C------------- LAST LINE OF CHFDV FOLLOWS ------------------------------
      END
      SUBROUTINE CHFEV(X1,X2,F1,F2,D1,D2,NE,XE,FE,NEXT,IERR)
C***BEGIN PROLOGUE  CHFEV
C***DATE WRITTEN   811019   (YYMMDD)
C***REVISION DATE  870707   (YYMMDD)
C***CATEGORY NO.  E3,H1
C***KEYWORDS  LIBRARY=SLATEC(PCHIP),
C             TYPE=SINGLE PRECISION(CHFEV-S DCHFEV-D),
C             CUBIC HERMITE EVALUATION,CUBIC POLYNOMIAL EVALUATION
C***AUTHOR  FRITSCH, F. N., (LLNL)
C             MATHEMATICS AND STATISTICS DIVISION
C             LAWRENCE LIVERMORE NATIONAL LABORATORY
C             P.O. BOX 808  (L-316)
C             LIVERMORE, CA  94550
C             FTS 532-4275, (415) 422-4275
C***PURPOSE  Evaluate a cubic polynomial given in Hermite form at an
C            array of points.  While designed for use by PCHFE, it may
C            be useful directly as an evaluator for a piecewise cubic
C            Hermite function in applications, such as graphing, where
C            the interval is known in advance.
C***DESCRIPTION
C
C          CHFEV:  Cubic Hermite Function EValuator
C
C     Evaluates the cubic polynomial determined by function values
C     F1,F2 and derivatives D1,D2 on interval (X1,X2) at the points
C     XE(J), J=1(1)NE.
C
C ----------------------------------------------------------------------
C
C  Calling sequence:
C
C        INTEGER  NE, NEXT(2), IERR
C        REAL  X1, X2, F1, F2, D1, D2, XE(NE), FE(NE)
C
C        CALL  CHFEV (X1,X2, F1,F2, D1,D2, NE, XE, FE, NEXT, IERR)
C
C   Parameters:
C
C     X1,X2 -- (input) endpoints of interval of definition of cubic.
C           (Error return if  X1.EQ.X2 .)
C
C     F1,F2 -- (input) values of function at X1 and X2, respectively.
C
C     D1,D2 -- (input) values of derivative at X1 and X2, respectively.
C
C     NE -- (input) number of evaluation points.  (Error return if
C           NE.LT.1 .)
C
C     XE -- (input) real array of points at which the function is to be
C           evaluated.  If any of the XE are outside the interval
C           [X1,X2], a warning error is returned in NEXT.
C
C     FE -- (output) real array of values of the cubic function defined
C           by  X1,X2, F1,F2, D1,D2  at the points  XE.
C
C     NEXT -- (output) integer array indicating number of extrapolation
C           points:
C            NEXT(1) = number of evaluation points to left of interval.
C            NEXT(2) = number of evaluation points to right of interval.
C
C     IERR -- (output) error flag.
C           Normal return:
C              IERR = 0  (no errors).
C           "Recoverable" errors:
C              IERR = -1  if NE.LT.1 .
C              IERR = -2  if X1.EQ.X2 .
C                (The FE-array has not been changed in either case.)
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  XERROR
C***END PROLOGUE  CHFEV
C
C ----------------------------------------------------------------------
C
C  Change record:
C     82-08-03   Minor cosmetic changes for release 1.
C
C ----------------------------------------------------------------------
C
C  Programming notes:
C
C     To produce a double precision version, simply:
C        a. Change CHFEV to DCHFEV wherever it occurs,
C        b. Change the real declaration to double precision,
C        c. Change the constant ZERO to double precision, and
C        d. Change the names of the Fortran functions:  AMAX1, AMIN1.
C
C  DECLARE ARGUMENTS.
C
      INTEGER  NE, NEXT(2), IERR
      REAL  X1, X2, F1, F2, D1, D2, XE(NE), FE(NE)
C
C  DECLARE LOCAL VARIABLES.
C
      INTEGER  I
      REAL  C2, C3, DEL1, DEL2, DELTA, H, X, XMI, XMA, ZERO
      DATA  ZERO /0./
C
C  VALIDITY-CHECK ARGUMENTS.
C
C***FIRST EXECUTABLE STATEMENT  CHFEV
      IF (NE .LT. 1)  GO TO 5001
      H = X2 - X1
      IF (H .EQ. ZERO)  GO TO 5002
C
C  INITIALIZE.
C
      IERR = 0
      NEXT(1) = 0
      NEXT(2) = 0
      XMI = AMIN1(ZERO, H)
      XMA = AMAX1(ZERO, H)
C
C  COMPUTE CUBIC COEFFICIENTS (EXPANDED ABOUT X1).
C
      DELTA = (F2 - F1)/H
      DEL1 = (D1 - DELTA)/H
      DEL2 = (D2 - DELTA)/H
C                                           (DELTA IS NO LONGER NEEDED.)
      C2 = -(DEL1+DEL1 + DEL2)
      C3 = (DEL1 + DEL2)/H
C                               (H, DEL1 AND DEL2 ARE NO LONGER NEEDED.)
C
C  EVALUATION LOOP.
C
      DO 500  I = 1, NE
         X = XE(I) - X1
         FE(I) = F1 + X*(D1 + X*(C2 + X*C3))
C          COUNT EXTRAPOLATION POINTS.
         IF ( X.LT.XMI )  NEXT(1) = NEXT(1) + 1
         IF ( X.GT.XMA )  NEXT(2) = NEXT(2) + 1
C        (NOTE REDUNDANCY--IF EITHER CONDITION IS TRUE, OTHER IS FALSE.)
  500 CONTINUE
C
C  NORMAL RETURN.
C
      RETURN
C
C  ERROR RETURNS.
C
 5001 CONTINUE
C     NE.LT.1 RETURN.
      IERR = -1
      CALL XERROR ('CHFEV -- NUMBER OF EVALUATION POINTS LESS THAN ONE'
     *           , 50, IERR, 1)
      RETURN
C
 5002 CONTINUE
C     X1.EQ.X2 RETURN.
      IERR = -2
      CALL XERROR ('CHFEV -- INTERVAL ENDPOINTS EQUAL'
     *           , 33, IERR, 1)
      RETURN
C------------- LAST LINE OF CHFEV FOLLOWS ------------------------------
      END
      REAL FUNCTION CHFIV(X1,X2,F1,F2,D1,D2,A,B,IERR)
C***BEGIN PROLOGUE  CHFIV
C***REFER TO  PCHIA
C***ROUTINES CALLED  XERROR
C***REVISION DATE  870707   (YYMMDD)
C***DESCRIPTION
C
C          CHFIV:  Cubic Hermite Function Integral Evaluator.
C
C     Called by  PCHIA  to evaluate the integral of a single cubic (in
C     Hermite form) over an arbitrary interval (A,B).
C
C ----------------------------------------------------------------------
C
C  Calling sequence:
C
C        INTEGER  IERR
C        REAL  X1, X2, F1, F2, D1, D2, A, B
C        REAL  VALUE, CHFIV
C
C        VALUE = CHFIV (X1, X2, F1, F2, D1, D2, A, B, IERR)
C
C   Parameters:
C
C     VALUE -- (output) VALUE of the requested integral.
C
C     X1,X2 -- (input) endpoints if interval of definition of cubic.
C           (Must be distinct.  Error return if not.)
C
C     F1,F2 -- (input) function values at the ends of the interval.
C
C     D1,D2 -- (input) derivative values at the ends of the interval.
C
C     A,B -- (input) endpoints of interval of integration.
C
C     IERR -- (output) error flag.
C           Normal return:
C              IERR = 0 (no errors).
C           "Recoverable errors":
C              IERR = -1  if X1.EQ.X2 .
C                (VALUE has not been set in this case.)
C
C***END PROLOGUE  CHFIV
C
C ----------------------------------------------------------------------
C
C  Programmed by:  Fred N. Fritsch,  FTS 532-4275, (415) 422-4275,
C                  Mathematics and Statistics Division,
C                  Lawrence Livermore National Laboratory.
C
C  Change record:
C     82-08-05   Converted to SLATEC library version.
C
C ----------------------------------------------------------------------
C
C  Programming notes:
C
C     To produce a double precision version, simply:
C        a. Change CHFIV to DCHFIV wherever it occurs,
C        b. Change the real declarations to double precision, and
C        c. Change the constants HALF, TWO, ... to double precision.
C
C  DECLARE ARGUMENTS.
C
      INTEGER  IERR
      REAL  X1, X2, F1, F2, D1, D2, A, B
C
C  DECLARE LOCAL VARIABLES.
C
      REAL  DTERM, FOUR, FTERM, H, HALF, PHIA1, PHIA2, PHIB1, PHIB2,
     *      PSIA1, PSIA2, PSIB1, PSIB2, TA1, TA2, TB1, TB2, THREE, TWO,
     *      UA1, UA2, UB1, UB2
C
C  INITIALIZE.
C
      DATA  HALF /0.5/,  TWO /2./,  THREE /3./,  FOUR /4./,  SIX /6./
C
C  VALIDITY CHECK INPUT.
C
C***FIRST EXECUTABLE STATEMENT  CHFIV
      IF (X1 .EQ. X2)  GO TO 5001
      IERR = 0
C
C  COMPUTE INTEGRAL.
C
      H = X2 - X1
      TA1 = (A - X1) / H
      TA2 = (X2 - A) / H
      TB1 = (B - X1) / H
      TB2 = (X2 - B) / H
C
      UA1 = TA1**3
      PHIA1 = UA1 * (TWO - TA1)
      PSIA1 = UA1 * (THREE*TA1 - FOUR)
      UA2 = TA2**3
      PHIA2 =  UA2 * (TWO - TA2)
      PSIA2 = -UA2 * (THREE*TA2 - FOUR)
C
      UB1 = TB1**3
      PHIB1 = UB1 * (TWO - TB1)
      PSIB1 = UB1 * (THREE*TB1 - FOUR)
      UB2 = TB2**3
      PHIB2 =  UB2 * (TWO - TB2)
      PSIB2 = -UB2 * (THREE*TB2 - FOUR)
C
      FTERM =   F1*(PHIA2 - PHIB2) + F2*(PHIB1 - PHIA1)
      DTERM = ( D1*(PSIA2 - PSIB2) + D2*(PSIB1 - PSIA1) )*(H/SIX)
C
C  RETURN VALUE.
C
      CHFIV = (HALF*H) * (FTERM + DTERM)
      RETURN
C
C  ERROR RETURN.
C
 5001 CONTINUE
      IERR = -1
      CALL XERROR ('CHFIV -- X1 EQUAL TO X2'
     *           , 23, IERR, 1)
      RETURN
C------------- LAST LINE OF CHFIV FOLLOWS ------------------------------
      END
      INTEGER FUNCTION CHFMC(D1,D2,DELTA)
C***BEGIN PROLOGUE  CHFMC
C***REFER TO  PCHMC
C***ROUTINES CALLED  R1MACH
C***REVISION DATE  870707   (YYMMDD)
C***DESCRIPTION
C
C          CHFMC:  Cubic Hermite Function Monotonicity Checker.
C
C    Called by  PCHMC  to determine the monotonicity properties of the
C    cubic with boundary derivative values D1,D2 and chord slope DELTA.
C
C ----------------------------------------------------------------------
C
C  Calling sequence:
C
C        REAL  D1, D2, DELTA
C        INTEGER  ISMON, CHFMC
C
C        ISMON = CHFMC (D1, D2, DELTA)
C
C   Parameters:
C
C     D1,D2 -- (input) derivative values at the ends of an interval.
C
C     DELTA -- (input) data slope over that interval.
C
C     ISMON -- (output) integer function value, indicating the monoto-
C           nicity of the cubic segment:
C             ISMON = -1  if function is strictly decreasing;
C             ISMON =  0  if function is constant;
C             ISMON =  1  if function is strictly increasing;
C             ISMON =  2  if function is non-monotonic;
C             ISMON =  3  if unable to determine.
C
C  Fortran intrinsics used:  SIGN.
C  Other routines used:  R1MACH.
C
C***END PROLOGUE  CHFMC
C
C ----------------------------------------------------------------------
C
C  Programmed by:  Fred N. Fritsch,  FTS 532-4275, (415) 422-4275,
C                  Mathematics and Statistics Division,
C                  Lawrence Livermore National Laboratory.
C
C  Change record:
C     82-08-05   Converted to SLATEC library version.
C     83-12-01   Changed from  ISIGN  to SIGN  to correct bug that
C                produced wrong sign when -1 .LT. DELTA .LT. 0 .
C
C ----------------------------------------------------------------------
C
C  Programming notes:
C
C     To produce a double precision version, simply:
C        a. Change CHFMC to DCHFMC wherever it occurs,
C        b. Change the real declarations to double precision, and
C        c. Change the constants ZERO, ONE, ... to double precision.
C
C  DECLARE ARGUMENTS.
C
      REAL  D1, D2, DELTA
C
C  DECLARE LOCAL VARIABLES.
C
      INTEGER  ISMON, ITRUE
      REAL  A, B, EPS, FOUR, ONE, PHI, TEN, THREE, TWO, ZERO
C
C  INITIALIZE.
C
      DATA  ZERO /0./,  ONE /1.0/,  TWO /2./,  THREE /3./,  FOUR /4./,
     1      TEN /10./
C
C        MACHINE-DEPENDENT PARAMETER -- SHOULD BE ABOUT 10*UROUND.
C***FIRST EXECUTABLE STATEMENT  CHFMC
      EPS = TEN*R1MACH(4)
C
C  MAKE THE CHECK.
C
      IF (DELTA .EQ. ZERO)  THEN
C        CASE OF CONSTANT DATA.
         IF ((D1.EQ.ZERO) .AND. (D2.EQ.ZERO))  THEN
            ISMON = 0
         ELSE
            ISMON = 2
         ENDIF
      ELSE
C        DATA IS NOT CONSTANT -- PICK UP SIGN.
         ITRUE = SIGN (ONE, DELTA)
         A = D1/DELTA
         B = D2/DELTA
         IF ((A.LT.ZERO) .OR. (B.LT.ZERO))  THEN
            ISMON = 2
         ELSE IF ((A.LE.THREE-EPS) .AND. (B.LE.THREE-EPS))  THEN
C           INSIDE SQUARE (0,3)X(0,3)  IMPLIES   OK.
            ISMON = ITRUE
         ELSE IF ((A.GT.FOUR+EPS) .AND. (B.GT.FOUR+EPS))  THEN
C           OUTSIDE SQUARE (0,4)X(0,4)  IMPLIES   NONMONOTONIC.
            ISMON = 2
         ELSE
C           MUST CHECK AGAINST BOUNDARY OF ELLIPSE.
            A = A - TWO
            B = B - TWO
            PHI = ((A*A + B*B) + A*B) - THREE
            IF (PHI .LT. -EPS)  THEN
               ISMON = ITRUE
            ELSE IF (PHI .GT. EPS)  THEN
               ISMON = 2
            ELSE
C              TO CLOSE TO BOUNDARY TO TELL,
C                  IN THE PRESENCE OF ROUND-OFF ERRORS.
               ISMON = 3
            ENDIF
         ENDIF
      ENDIF
C
C  RETURN VALUE.
C
      CHFMC = ISMON
      RETURN
C------------- LAST LINE OF CHFMC FOLLOWS ------------------------------
      END
      SUBROUTINE PCHCE(IC,VC,N,X,H,SLOPE,D,INCFD,IERR)
C***BEGIN PROLOGUE  PCHCE
C***REFER TO  PCHIC
C***ROUTINES CALLED  PCHDF,PCHST,XERROR
C***REVISION DATE  870707   (YYMMDD)
C***DESCRIPTION
C
C          PCHCE:  PCHIC End Derivative Setter.
C
C    Called by PCHIC to set end derivatives as requested by the user.
C    It must be called after interior derivative values have been set.
C                      -----
C
C    To facilitate two-dimensional applications, includes an increment
C    between successive values of the D-array.
C
C ----------------------------------------------------------------------
C
C  Calling sequence:
C
C        PARAMETER  (INCFD = ...)
C        INTEGER  IC(2), N, IERR
C        REAL  VC(2), X(N), H(N), SLOPE(N), D(INCFD,N)
C
C        CALL  PCHCE (IC, VC, N, X, H, SLOPE, D, INCFD, IERR)
C
C   Parameters:
C
C     IC -- (input) integer array of length 2 specifying desired
C           boundary conditions:
C           IC(1) = IBEG, desired condition at beginning of data.
C           IC(2) = IEND, desired condition at end of data.
C           ( see prologue to PCHIC for details. )
C
C     VC -- (input) real array of length 2 specifying desired boundary
C           values.  VC(1) need be set only if IC(1) = 2 or 3 .
C                    VC(2) need be set only if IC(2) = 2 or 3 .
C
C     N -- (input) number of data points.  (assumes N.GE.2)
C
C     X -- (input) real array of independent variable values.  (the
C           elements of X are assumed to be strictly increasing.)
C
C     H -- (input) real array of interval lengths.
C     SLOPE -- (input) real array of data slopes.
C           If the data are (X(I),Y(I)), I=1(1)N, then these inputs are:
C                  H(I) =  X(I+1)-X(I),
C              SLOPE(I) = (Y(I+1)-Y(I))/H(I),  I=1(1)N-1.
C
C     D -- (input) real array of derivative values at the data points.
C           The value corresponding to X(I) must be stored in
C                D(1+(I-1)*INCFD),  I=1(1)N.
C          (output) the value of D at X(1) and/or X(N) is changed, if
C           necessary, to produce the requested boundary conditions.
C           no other entries in D are changed.
C
C     INCFD -- (input) increment between successive values in D.
C           This argument is provided primarily for 2-D applications.
C
C     IERR -- (output) error flag.
C           Normal return:
C              IERR = 0  (no errors).
C           Warning errors:
C              IERR = 1  if IBEG.LT.0 and D(1) had to be adjusted for
C                        monotonicity.
C              IERR = 2  if IEND.LT.0 and D(1+(N-1)*INCFD) had to be
C                        adjusted for monotonicity.
C              IERR = 3  if both of the above are true.
C
C    -------
C    WARNING:  This routine does no validity-checking of arguments.
C    -------
C
C  Fortran intrinsics used:  ABS, IABS.
C
C***END PROLOGUE  PCHCE
C
C ----------------------------------------------------------------------
C
C  Programmed by:  Fred N. Fritsch,  FTS 532-4275, (415) 422-4275,
C                  Mathematics and Statistics Division,
C                  Lawrence Livermore National Laboratory.
C
C  Change record:
C     82-08-05   Converted to SLATEC library version.
C     87-07-07   Minor corrections made to prologue..
C
C ----------------------------------------------------------------------
C
C  Programming notes:
C
C     1. The function  PCHST(ARG1,ARG2)  is assumed to return zero if
C        either argument is zero, +1 if they are of the same sign, and
C        -1 if they are of opposite sign.
C     2. To produce a double precision version, simply:
C        a. Change PCHCE to DPCHCE wherever it occurs,
C        b. Change PCHDF to DPCHCE wherever it occurs,
C        c. Change PCHST to DPCHST wherever it occurs,
C        d. Change all references to the Fortran intrinsics to their
C           double presision equivalents,
C        e. Change the real declarations to double precision, and
C        f. Change the constants ZERO, HALF, ... to double precision.
C     3. One could reduce the number of arguments and amount of local
C        storage, at the expense of reduced code clarity, by passing in
C        the array WK (rather than splitting it into H and SLOPE) and
C        increasing its length enough to incorporate STEMP and XTEMP.
C     4. The two monotonicity checks only use the sufficient conditions.
C        Thus, it is possible (but unlikely) for a boundary condition to
C        be changed, even though the original interpolant was monotonic.
C        (At least the result is a continuous function of the data.)
C
C  DECLARE ARGUMENTS.
C
      INTEGER  IC(2), N, INCFD, IERR
      REAL  VC(2), X(N), H(N), SLOPE(N), D(INCFD,N)
C
C  DELCARE LOCAL VARIABLES.
C
      INTEGER  IBEG, IEND, IERF, INDEX, J, K
      REAL  HALF, STEMP(3), THREE, TWO, XTEMP(4), ZERO
      REAL  PCHDF, PCHST
C
C  INITIALIZE.
C
      DATA  ZERO /0./,  HALF /0.5/,  TWO /2./,  THREE /3./
C
C***FIRST EXECUTABLE STATEMENT  PCHCE
      IBEG = IC(1)
      IEND = IC(2)
      IERR = 0
C
C  SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL.
C
      IF ( IABS(IBEG).GT.N )  IBEG = 0
      IF ( IABS(IEND).GT.N )  IEND = 0
C
C  TREAT BEGINNING BOUNDARY CONDITION.
C
      IF (IBEG .EQ. 0)  GO TO 2000
      K = IABS(IBEG)
      IF (K .EQ. 1)  THEN
C        BOUNDARY VALUE PROVIDED.
         D(1,1) = VC(1)
      ELSE IF (K .EQ. 2)  THEN
C        BOUNDARY SECOND DERIVATIVE PROVIDED.
         D(1,1) = HALF*( (THREE*SLOPE(1) - D(1,2)) - HALF*VC(1)*H(1) )
      ELSE IF (K .LT. 5)  THEN
C        USE K-POINT DERIVATIVE FORMULA.
C        PICK UP FIRST K POINTS, IN REVERSE ORDER.
         DO 10  J = 1, K
            INDEX = K-J+1
C           INDEX RUNS FROM K DOWN TO 1.
            XTEMP(J) = X(INDEX)
            IF (J .LT. K)  STEMP(J) = SLOPE(INDEX-1)
   10    CONTINUE
C                 -----------------------------
         D(1,1) = PCHDF (K, XTEMP, STEMP, IERF)
C                 -----------------------------
         IF (IERF .NE. 0)  GO TO 5001
      ELSE
C        USE 'NOT A KNOT' CONDITION.
         D(1,1) = ( THREE*(H(1)*SLOPE(2) + H(2)*SLOPE(1))
     *             - TWO*(H(1)+H(2))*D(1,2) - H(1)*D(1,3) ) / H(2)
      ENDIF
C
      IF (IBEG .GT. 0)  GO TO 2000
C
C  CHECK D(1,1) FOR COMPATIBILITY WITH MONOTONICITY.
C
      IF (SLOPE(1) .EQ. ZERO)  THEN
         IF (D(1,1) .NE. ZERO)  THEN
            D(1,1) = ZERO
            IERR = IERR + 1
         ENDIF
      ELSE IF ( PCHST(D(1,1),SLOPE(1)) .LT. ZERO)  THEN
         D(1,1) = ZERO
         IERR = IERR + 1
      ELSE IF ( ABS(D(1,1)) .GT. THREE*ABS(SLOPE(1)) )  THEN
         D(1,1) = THREE*SLOPE(1)
         IERR = IERR + 1
      ENDIF
C
C  TREAT END BOUNDARY CONDITION.
C
 2000 CONTINUE
      IF (IEND .EQ. 0)  GO TO 5000
      K = IABS(IEND)
      IF (K .EQ. 1)  THEN
C        BOUNDARY VALUE PROVIDED.
         D(1,N) = VC(2)
      ELSE IF (K .EQ. 2)  THEN
C        BOUNDARY SECOND DERIVATIVE PROVIDED.
         D(1,N) = HALF*( (THREE*SLOPE(N-1) - D(1,N-1)) +
     *                                           HALF*VC(2)*H(N-1) )
      ELSE IF (K .LT. 5)  THEN
C        USE K-POINT DERIVATIVE FORMULA.
C        PICK UP LAST K POINTS.
         DO 2010  J = 1, K
            INDEX = N-K+J
C           INDEX RUNS FROM N+1-K UP TO N.
            XTEMP(J) = X(INDEX)
            IF (J .LT. K)  STEMP(J) = SLOPE(INDEX)
 2010    CONTINUE
C                 -----------------------------
         D(1,N) = PCHDF (K, XTEMP, STEMP, IERF)
C                 -----------------------------
         IF (IERF .NE. 0)  GO TO 5001
      ELSE
C        USE 'NOT A KNOT' CONDITION.
         D(1,N) = ( THREE*(H(N-1)*SLOPE(N-2) + H(N-2)*SLOPE(N-1))
     *             - TWO*(H(N-1)+H(N-2))*D(1,N-1) - H(N-1)*D(1,N-2) )
     *                                                         / H(N-2)
      ENDIF
C
      IF (IEND .GT. 0)  GO TO 5000
C
C  CHECK D(1,N) FOR COMPATIBILITY WITH MONOTONICITY.
C
      IF (SLOPE(N-1) .EQ. ZERO)  THEN
         IF (D(1,N) .NE. ZERO)  THEN
            D(1,N) = ZERO
            IERR = IERR + 2
         ENDIF
      ELSE IF ( PCHST(D(1,N),SLOPE(N-1)) .LT. ZERO)  THEN
         D(1,N) = ZERO
         IERR = IERR + 2
      ELSE IF ( ABS(D(1,N)) .GT. THREE*ABS(SLOPE(N-1)) )  THEN
         D(1,N) = THREE*SLOPE(N-1)
         IERR = IERR + 2
      ENDIF
C
C  NORMAL RETURN.
C
 5000 CONTINUE
      RETURN
C
C  ERROR RETURN.
C
 5001 CONTINUE
C     ERROR RETURN FROM PCHDF.
C   *** THIS CASE SHOULD NEVER OCCUR ***
      IERR = -1
      CALL XERROR ('PCHCE -- ERROR RETURN FROM PCHDF'
     *           , 32, IERR, 1)
      RETURN
C------------- LAST LINE OF PCHCE FOLLOWS ------------------------------
      END
      SUBROUTINE PCHCI(N,H,SLOPE,D,INCFD)
C***BEGIN PROLOGUE  PCHCI
C***REFER TO  PCHIC
C***ROUTINES CALLED  PCHST
C***REVISION DATE  870707   (YYMMDD)
C***DESCRIPTION
C
C          PCHCI:  PCHIC Initial Derivative Setter.
C
C    Called by PCHIC to set derivatives needed to determine a monotone
C    piecewise cubic Hermite interpolant to the data.
C
C    Default boundary conditions are provided which are compatible
C    with monotonicity.  If the data are only piecewise monotonic, the
C    interpolant will have an extremum at each point where monotonicity
C    switches direction.
C
C    To facilitate two-dimensional applications, includes an increment
C    between successive values of the D-array.
C
C    The resulting piecewise cubic Hermite function should be identical
C    (within roundoff error) to that produced by PCHIM.
C
C ----------------------------------------------------------------------
C
C  Calling sequence:
C
C        PARAMETER  (INCFD = ...)
C        INTEGER  N
C        REAL  H(N), SLOPE(N), D(INCFD,N)
C
C        CALL  PCHCI (N, H, SLOPE, D, INCFD)
C
C   Parameters:
C
C     N -- (input) number of data points.
C           If N=2, simply does linear interpolation.
C
C     H -- (input) real array of interval lengths.
C     SLOPE -- (input) real array of data slopes.
C           If the data are (X(I),Y(I)), I=1(1)N, then these inputs are:
C                  H(I) =  X(I+1)-X(I),
C              SLOPE(I) = (Y(I+1)-Y(I))/H(I),  I=1(1)N-1.
C
C     D -- (output) real array of derivative values at the data points.
C           If the data are monotonic, these values will determine a
C           a monotone cubic Hermite function.
C           The value corresponding to X(I) is stored in
C                D(1+(I-1)*INCFD),  I=1(1)N.
C           No other entries in D are changed.
C
C     INCFD -- (input) increment between successive values in D.
C           This argument is provided primarily for 2-D applications.
C
C    -------
C    WARNING:  This routine does no validity-checking of arguments.
C    -------
C
C  Fortran intrinsics used:  ABS, AMAX1, AMIN1.
C
C***END PROLOGUE  PCHCI
C
C ----------------------------------------------------------------------
C
C  Programmed by:  Fred N. Fritsch,  FTS 532-4275, (415) 422-4275,
C                  Mathematics and Statistics Division,
C                  Lawrence Livermore National Laboratory.
C
C  Change record:
C     82-06-01   Modified end conditions to be continuous functions of
C                data when monotonicity switches in next interval.
C     82-06-02   1. Modified formulas so end conditions are less prone
C                   to over/underflow problems.
C                2. Minor modification to HSUM calculation.
C     82-08-05   Converted to SLATEC library version.
C
C ----------------------------------------------------------------------
C
C  Programming notes:
C
C     1. The function  PCHST(ARG1,ARG2)  is assumed to return zero if
C        either argument is zero, +1 if they are of the same sign, and
C        -1 if they are of opposite sign.
C     2. To produce a double precision version, simply:
C        a. Change PCHCI to DPCHCI wherever it occurs,
C        b. Change PCHST to DPCHST wherever it occurs,
C        c. Change all references to the Fortran intrinsics to their
C           double presision equivalents,
C        d. Change the real declarations to double precision, and
C        e. Change the constants ZERO and THREE to double precision.
C
C  DECLARE ARGUMENTS.
C
      INTEGER  N, INCFD
      REAL  H(N), SLOPE(N), D(INCFD,N)
C
C  DECLARE LOCAL VARIBLES.
C
      INTEGER  I, NLESS1
      REAL  DEL1, DEL2, DMAX, DMIN, DRAT1, DRAT2, HSUM, HSUMT3, THREE,
     *      W1, W2, ZERO
      REAL  PCHST
C
C  INITIALIZE.
C
      DATA  ZERO /0./,  THREE /3./
C***FIRST EXECUTABLE STATEMENT  PCHCI
      NLESS1 = N - 1
      DEL1 = SLOPE(1)
C
C  SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION.
C
      IF (NLESS1 .GT. 1)  GO TO 10
      D(1,1) = DEL1
      D(1,N) = DEL1
      GO TO 5000
C
C  NORMAL CASE  (N .GE. 3).
C
   10 CONTINUE
      DEL2 = SLOPE(2)
C
C  SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE
C     SHAPE-PRESERVING.
C
      HSUM = H(1) + H(2)
      W1 = (H(1) + HSUM)/HSUM
      W2 = -H(1)/HSUM
      D(1,1) = W1*DEL1 + W2*DEL2
      IF ( PCHST(D(1,1),DEL1) .LE. ZERO)  THEN
         D(1,1) = ZERO
      ELSE IF ( PCHST(DEL1,DEL2) .LT. ZERO)  THEN
C        NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES.
         DMAX = THREE*DEL1
         IF (ABS(D(1,1)) .GT. ABS(DMAX))  D(1,1) = DMAX
      ENDIF
C
C  LOOP THROUGH INTERIOR POINTS.
C
      DO 50  I = 2, NLESS1
         IF (I .EQ. 2)  GO TO 40
C
         HSUM = H(I-1) + H(I)
         DEL1 = DEL2
         DEL2 = SLOPE(I)
   40    CONTINUE
C
C        SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC.
C
         D(1,I) = ZERO
         IF ( PCHST(DEL1,DEL2) .LE. ZERO)  GO TO 50
C
C        USE BRODLIE MODIFICATION OF BUTLAND FORMULA.
C
         HSUMT3 = HSUM+HSUM+HSUM
         W1 = (HSUM + H(I-1))/HSUMT3
         W2 = (HSUM + H(I)  )/HSUMT3
         DMAX = AMAX1( ABS(DEL1), ABS(DEL2) )
         DMIN = AMIN1( ABS(DEL1), ABS(DEL2) )
         DRAT1 = DEL1/DMAX
         DRAT2 = DEL2/DMAX
         D(1,I) = DMIN/(W1*DRAT1 + W2*DRAT2)
C
   50 CONTINUE
C
C  SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE
C     SHAPE-PRESERVING.
C
      W1 = -H(N-1)/HSUM
      W2 = (H(N-1) + HSUM)/HSUM
      D(1,N) = W1*DEL1 + W2*DEL2
      IF ( PCHST(D(1,N),DEL2) .LE. ZERO)  THEN
         D(1,N) = ZERO
      ELSE IF ( PCHST(DEL1,DEL2) .LT. ZERO)  THEN
C        NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES.
         DMAX = THREE*DEL2
         IF (ABS(D(1,N)) .GT. ABS(DMAX))  D(1,N) = DMAX
      ENDIF
C
C  NORMAL RETURN.
C
 5000 CONTINUE
      RETURN
C------------- LAST LINE OF PCHCI FOLLOWS ------------------------------
      END
      SUBROUTINE PCHCS(SWITCH,N,H,SLOPE,D,INCFD,IERR)
C***BEGIN PROLOGUE  PCHCS
C***REFER TO  PCHIC
C***ROUTINES CALLED  PCHST,PCHSW
C***DESCRIPTION
C
C         PCHCS:  PCHIC Monotonicity Switch Derivative Setter.
C
C     Called by  PCHIC  to adjust the values of D in the vicinity of a
C     switch in direction of monotonicity, to produce a more "visually
C     pleasing" curve than that given by  PCHIM .
C
C ----------------------------------------------------------------------
C
C  Calling sequence:
C
C        PARAMETER  (INCFD = ...)
C        INTEGER  N, IERR
C        REAL  SWITCH, H(N), SLOPE(N), D(INCFD,N)
C
C        CALL  PCHCS (SWITCH, N, H, SLOPE, D, INCFD, IERR)
C
C   Parameters:
C
C     SWITCH -- (input) indicates the amount of control desired over
C           local excursions from data.
C
C     N -- (input) number of data points.  (assumes N.GT.2 .)
C
C     H -- (input) real array of interval lengths.
C     SLOPE -- (input) real array of data slopes.
C           If the data are (X(I),Y(I)), I=1(1)N, then these inputs are:
C                  H(I) =  X(I+1)-X(I),
C              SLOPE(I) = (Y(I+1)-Y(I))/H(I),  I=1(1)N-1.
C
C     D -- (input) real array of derivative values at the data points,
C           as determined by PCHCI.
C          (output) derivatives in the vicinity of switches in direction
C           of monotonicity may be adjusted to produce a more "visually
C           pleasing" curve.
C           The value corresponding to X(I) is stored in
C                D(1+(I-1)*INCFD),  I=1(1)N.
C           No other entries in D are changed.
C
C     INCFD -- (input) increment between successive values in D.
C           This argument is provided primarily for 2-D applications.
C
C     IERR -- (output) error flag.  should be zero.
C           If negative, trouble in PCHSW.  (should never happen.)
C
C    -------
C    WARNING:  This routine does no validity-checking of arguments.
C    -------
C
C  Fortran intrinsics used:  ABS, AMAX1, AMIN1.
C
C***END PROLOGUE  PCHCS
C
C ----------------------------------------------------------------------
C
C  Programmed by:  Fred N. Fritsch,  FTS 532-4275, (415) 422-4275,
C                  Mathematics and Statistics Division,
C                  Lawrence Livermore National Laboratory.
C
C  Change record:
C     82-06-17   Redesigned to (1) fix  problem with lack of continuity
C                approaching a flat-topped peak (2) be cleaner and
C                easier to verify.
C                Eliminated subroutines PCHSA and PCHSX in the process.
C     82-06-22   1. Limited fact to not exceed one, so computed D is a
C                   convex combination of PCHCI value and PCHSD value.
C                2. Changed fudge from 1 to 4 (based on experiments).
C     82-06-23   Moved PCHSD to an inline function (eliminating MSWTYP).
C     82-08-05   Converted to SLATEC library version.
C
C ----------------------------------------------------------------------
C
C  Programming notes:
C
C     1. The function  PCHST(ARG1,ARG2)  is assumed to return zero if
C        either argument is zero, +1 if they are of the same sign, and
C        -1 if they are of opposite sign.
C     2. To produce a double precision version, simply:
C        a. Change PCHCS to DPCHCS wherever it occurs,
C        b. Change PCHSD to dpchsd wherever it occurs,
C        c. Change PCHST to DPCHST wherever it occurs,
C        d. Change PCHSW to DPCHSW wherever it occurs,
C        e. Change all references to the Fortran intrinsics to their
C           double precision equivalents,
C        f. Change the real declarations to double precision, and
C        g. Change the constants ZERO and ONE to double precision.
C
C  DECLARE ARGUMENTS.
C
      INTEGER  N, INCFD, IERR
      REAL  SWITCH, H(N), SLOPE(N), D(INCFD,N)
C
C  DECLARE LOCAL VARIABLES.
C
      INTEGER  I, INDX, K, NLESS1
      REAL  DEL(3), DEXT, DFLOC, DFMX, FACT, FUDGE, ONE, SLMAX,
     *      WTAVE(2), ZERO
      REAL  PCHST
C
C  DEFINE INLINE FUNCTION FOR WEIGHTED AVERAGE OF SLOPES.
C
      REAL  PCHSD, S1, S2, H1, H2
      PCHSD(S1,S2,H1,H2) = (H2/(H1+H2))*S1 + (H1/(H1+H2))*S2
C
C  INITIALIZE.
C
      DATA  ZERO /0./,  ONE /1./
      DATA  FUDGE /4./
C***FIRST EXECUTABLE STATEMENT  PCHCS
      IERR = 0
      NLESS1 = N - 1
C
C  LOOP OVER SEGMENTS.
C
      DO 900  I = 2, NLESS1
         IF ( PCHST(SLOPE(I-1),SLOPE(I)) )  100, 300, 900
C             --------------------------
C
  100    CONTINUE
C
C....... SLOPE SWITCHES MONOTONICITY AT I-TH POINT .....................
C
C           DO NOT CHANGE D IF 'UP-DOWN-UP'.
            IF (I .GT. 2)  THEN
               IF ( PCHST(SLOPE(I-2),SLOPE(I)) .GT. ZERO)  GO TO 900
C                   --------------------------
            ENDIF
            IF (I .LT. NLESS1)  THEN
               IF ( PCHST(SLOPE(I+1),SLOPE(I-1)) .GT. ZERO)  GO TO 900
C                   ----------------------------
            ENDIF
C
C   ....... COMPUTE PROVISIONAL VALUE FOR D(1,I).
C
            DEXT = PCHSD (SLOPE(I-1), SLOPE(I), H(I-1), H(I))
C
C   ....... DETERMINE WHICH INTERVAL CONTAINS THE EXTREMUM.
C
            IF ( PCHST(DEXT, SLOPE(I-1)) )  200, 900, 250
C                -----------------------
C
  200       CONTINUE
C              DEXT AND SLOPE(I-1) HAVE OPPOSITE SIGNS --
C                        EXTREMUM IS IN (X(I-1),X(I)).
               K = I-1
C              SET UP TO COMPUTE NEW VALUES FOR D(1,I-1) AND D(1,I).
               WTAVE(2) = DEXT
               IF (K .GT. 1)
     *            WTAVE(1) = PCHSD (SLOPE(K-1), SLOPE(K), H(K-1), H(K))
               GO TO 400
C
  250       CONTINUE
C              DEXT AND SLOPE(I) HAVE OPPOSITE SIGNS --
C                        EXTREMUM IS IN (X(I),X(I+1)).
               K = I
C              SET UP TO COMPUTE NEW VALUES FOR D(1,I) AND D(1,I+1).
               WTAVE(1) = DEXT
               IF (K .LT. NLESS1)
     *            WTAVE(2) = PCHSD (SLOPE(K), SLOPE(K+1), H(K), H(K+1))
               GO TO 400
C
  300    CONTINUE
C
C....... AT LEAST ONE OF SLOPE(I-1) AND SLOPE(I) IS ZERO --
C                     CHECK FOR FLAT-TOPPED PEAK .......................
C
            IF (I .EQ. NLESS1)  GO TO 900
            IF ( PCHST(SLOPE(I-1), SLOPE(I+1)) .GE. ZERO)  GO TO 900
C                -----------------------------
C
C           WE HAVE FLAT-TOPPED PEAK ON (X(I),X(I+1)).
            K = I
C           SET UP TO COMPUTE NEW VALUES FOR D(1,I) AND D(1,I+1).
            WTAVE(1) = PCHSD (SLOPE(K-1), SLOPE(K), H(K-1), H(K))
            WTAVE(2) = PCHSD (SLOPE(K), SLOPE(K+1), H(K), H(K+1))
C
  400    CONTINUE
C
C....... AT THIS POINT WE HAVE DETERMINED THAT THERE WILL BE AN EXTREMUM
C        ON (X(K),X(K+1)), WHERE K=I OR I-1, AND HAVE SET ARRAY WTAVE--
C           WTAVE(1) IS A WEIGHTED AVERAGE OF SLOPE(K-1) AND SLOPE(K),
C                    IF K.GT.1
C           WTAVE(2) IS A WEIGHTED AVERAGE OF SLOPE(K) AND SLOPE(K+1),
C                    IF K.LT.N-1
C
         SLMAX = ABS(SLOPE(K))
         IF (K .GT. 1)    SLMAX = AMAX1( SLMAX, ABS(SLOPE(K-1)) )
         IF (K.LT.NLESS1) SLMAX = AMAX1( SLMAX, ABS(SLOPE(K+1)) )
C
         IF (K .GT. 1)  DEL(1) = SLOPE(K-1) / SLMAX
         DEL(2) = SLOPE(K) / SLMAX
         IF (K.LT.NLESS1)  DEL(3) = SLOPE(K+1) / SLMAX
C
         IF ((K.GT.1) .AND. (K.LT.NLESS1))  THEN
C           NORMAL CASE -- EXTREMUM IS NOT IN A BOUNDARY INTERVAL.
            FACT = FUDGE* ABS(DEL(3)*(DEL(1)-DEL(2))*(WTAVE(2)/SLMAX))
            D(1,K) = D(1,K) + AMIN1(FACT,ONE)*(WTAVE(1) - D(1,K))
            FACT = FUDGE* ABS(DEL(1)*(DEL(3)-DEL(2))*(WTAVE(1)/SLMAX))
            D(1,K+1) = D(1,K+1) + AMIN1(FACT,ONE)*(WTAVE(2) - D(1,K+1))
         ELSE
C           SPECIAL CASE K=1 (WHICH CAN OCCUR ONLY IF I=2) OR
C                        K=NLESS1 (WHICH CAN OCCUR ONLY IF I=NLESS1).
            FACT = FUDGE* ABS(DEL(2))
            D(1,I) = AMIN1(FACT,ONE) * WTAVE(I-K+1)
C              NOTE THAT I-K+1 = 1 IF K=I  (=NLESS1),
C                        I-K+1 = 2 IF K=I-1(=1).
         ENDIF
C
C
C....... ADJUST IF NECESSARY TO LIMIT EXCURSIONS FROM DATA.
C
         IF (SWITCH .LE. ZERO)  GO TO 900
C
         DFLOC = H(K)*ABS(SLOPE(K))
         IF (K .GT. 1)    DFLOC = AMAX1( DFLOC, H(K-1)*ABS(SLOPE(K-1)) )
         IF (K.LT.NLESS1) DFLOC = AMAX1( DFLOC, H(K+1)*ABS(SLOPE(K+1)) )
         DFMX = SWITCH*DFLOC
         INDX = I-K+1
C        INDX = 1 IF K=I, 2 IF K=I-1.
C        ---------------------------------------------------------------
         CALL PCHSW (DFMX, INDX, D(1,K), D(1,K+1), H(K), SLOPE(K), IERR)
C        ---------------------------------------------------------------
         IF (IERR .NE. 0)  RETURN
C
C....... END OF SEGMENT LOOP.
C
  900 CONTINUE
C
      RETURN
C------------- LAST LINE OF PCHCS FOLLOWS ------------------------------
      END
      REAL FUNCTION PCHDF(K,X,S,IERR)
C***BEGIN PROLOGUE  PCHDF
C***REFER TO  PCHCE,PCHSP
C***ROUTINES CALLED  XERROR
C***DESCRIPTION
C
C          PCHDF:   PCHIP Finite Difference Formula
C
C     Uses a divided difference formulation to compute a K-point approx-
C     imation to the derivative at X(K) based on the data in X and S.
C
C     Called by  PCHCE  and  PCHSP  to compute 3- and 4-point boundary
C     derivative approximations.
C
C ----------------------------------------------------------------------
C
C     On input:
C        K      is the order of the desired derivative approximation.
C               K must be at least 3 (error return if not).
C        X      contains the K values of the independent variable.
C               X need not be ordered, but the values **MUST** be
C               distinct.  (Not checked here.)
C        S      contains the associated slope values:
C                  S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1.
C               (Note that S need only be of length K-1.)
C
C     On return:
C        S      will be destroyed.
C        IERR   will be set to -1 if K.LT.2 .
C        PCHDF  will be set to the desired derivative approximation if
C               IERR=0 or to zero if IERR=-1.
C
C ----------------------------------------------------------------------
C
C  Reference:  Carl de Boor, A Practical Guide to Splines, Springer-
C              Verlag (New York, 1978), pp. 10-16.
C
C ----------------------------------------------------------------------
C
C  Programmed by:  Fred N. Fritsch,  FTS 532-4275, (415) 422-4275,
C                  Mathematics and Statistics Division,
C                  Lawrence Livermore National Laboratory.
C
C  Change record:
C     82-08-05   Converted to SLATEC library version.
C
C ----------------------------------------------------------------------
C
C  Programming notes:
C
C     To produce a double precision version, simply:
C        a. Change PCHDF to DPCHDF wherever it occurs,
C        b. Change the real declarations to double precision, and
C        c. Change the constant ZERO to double precision.
C***END PROLOGUE  PCHDF
      INTEGER  K, IERR
      REAL  X(K), S(K)
C
C  DECLARE LOCAL VARIABLES.
C
      INTEGER  I, J
      REAL  VALUE, ZERO
      DATA  ZERO /0./
C
C  CHECK FOR LEGAL VALUE OF K.
C
C***FIRST EXECUTABLE STATEMENT  PCHDF
      IF (K .LT. 3)  GO TO 5001
C
C  COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL.
C
      DO 10  J = 2, K-1
         DO 9  I = 1, K-J
            S(I) = (S(I+1)-S(I))/(X(I+J)-X(I))
    9    CONTINUE
   10 CONTINUE
C
C  EVALUATE DERIVATIVE AT X(K).
C
      VALUE = S(1)
      DO 20  I = 2, K-1
         VALUE = S(I) + VALUE*(X(K)-X(I))
   20 CONTINUE
C
C  NORMAL RETURN.
C
      IERR = 0
      PCHDF = VALUE
      RETURN
C
C  ERROR RETURN.
C
 5001 CONTINUE
C     K.LT.3 RETURN.
      IERR = -1
      CALL XERROR ('PCHDF -- K LESS THAN THREE'
     *           , 26, IERR, 1)
      PCHDF = ZERO
      RETURN
C------------- LAST LINE OF PCHDF FOLLOWS ------------------------------
      END
      SUBROUTINE PCHEV(N,X,F,D,NVAL,XVAL,FVAL,DVAL,IERR)
C***BEGIN PROLOGUE  PCHEV
C***DATE WRITTEN   870828   (YYMMDD)
C***REVISION DATE  870828   (YYMMDD)
C***CATEGORY NO.  E3,H1
C***KEYWORDS  CUBIC HERMITE OR SPLINE DIFFERENTIATION,CUBIC HERMITE
C             EVALUATION,EASY TO USE SPLINE OR CUBIC HERMITE EVALUATOR
C***AUTHOR  KAHANER, D.K., (NBS)
C             SCIENTIFIC COMPUTING DIVISION
C             NATIONAL BUREAU OF STANDARDS
C             ROOM A161, TECHNOLOGY BUILDING
C             GAITHERSBURG, MARYLAND 20899
C             (301) 975-3808
C***PURPOSE  Evaluates the function and first derivative of a piecewise
C            cubic Hermite or spline function at an array of points XVAL,
C            easy to use.
C***DESCRIPTION
C
C          PCHEV:  Piecewise Cubic Hermite or Spline Derivative Evaluator,
C                  Easy to Use.
C
C     From the book "Numerical Methods and Software"
C          by  D. Kahaner, C. Moler, S. Nash
C                 Prentice Hall 1988
C
C     Evaluates the function and first derivative of the cubic Hermite
C     or spline function defined by  N, X, F, D, at the array of points XVAL.
C
C     This is an easy to use driver for the routines by F.N. Fritsch
C     described in reference (2) below. Those also have other capabilities.
C
C ----------------------------------------------------------------------
C
C  Calling sequence: CALL  PCHEV (N, X, F, D, NVAL, XVAL, FVAL, DVAL, IERR)
C
C     INTEGER  N, NVAL, IERR
C     REAL  X(N), F(N), D(N), XVAL(NVAL), FVAL(NVAL), DVAL(NVAL)
C
C   Parameters:
C
C     N -- (input) number of data points.  (Error return if N.LT.2 .)
C
C     X -- (input) real array of independent variable values.  The
C           elements of X must be strictly increasing:
C             X(I-1) .LT. X(I),  I = 2(1)N. (Error return if not.)
C
C     F -- (input) real array of function values.  F(I) is
C           the value corresponding to X(I).
C
C     D -- (input) real array of derivative values.  D(I) is
C           the value corresponding to X(I).
C
C  NVAL -- (input) number of points at which the functions are to be
C           evaluated. ( Error return if NVAL.LT.1 )
C
C  XVAL -- (input) real array of points at which the functions are to
C           be evaluated.
C
C          NOTES:
C           1. The evaluation will be most efficient if the elements
C              of XVAL are increasing relative to X;
C              that is,   XVAL(J) .GE. X(I)
C              implies    XVAL(K) .GE. X(I),  all K.GE.J .
C           2. If any of the XVAL are outside the interval [X(1),X(N)],
C              values are extrapolated from the nearest extreme cubic,
C              and a warning error is returned.
C
C  FVAL -- (output) real array of values of the cubic Hermite function
C           defined by  N, X, F, D  at the points  XVAL.
C
C  DVAL -- (output) real array of values of the first derivative of
C           the same function at the points  XVAL.
C
C  IERR -- (output) error flag.
C           Normal return:
C              IERR = 0  (no errors).
C           Warning error:
C              IERR.GT.0  means that extrapolation was performed at
C                 IERR points.
C           "Recoverable" errors:
C              IERR = -1  if N.LT.2 .
C              IERR = -3  if the X-array is not strictly increasing.
C              IERR = -4  if NVAL.LT.1 .
C           (Output arrays have not been changed in any of these cases.)
C               NOTE:  The above errors are checked in the order listed,
C                   and following arguments have **NOT** been validated.
C              IERR = -5  if an error has occurred in the lower-level
C                         routine CHFDV.  NB: this should never happen.
C                         Notify the author **IMMEDIATELY** if it does.
C
C ----------------------------------------------------------------------
C***REFERENCES  1. F.N.FRITSCH AND R.E.CARLSON, 'MONOTONE PIECEWISE
C                 CUBIC INTERPOLATION,' SIAM J.NUMER.ANAL. 17, 2 (APRIL
C                 1980), 238-246.
C               2. F.N.FRITSCH, 'PIECEWISE CUBIC HERMITE INTERPOLATION
C                 PACKAGE, FINAL SPECIFICATIONS', LAWRENCE LIVERMORE
C                 NATIONAL LABORATORY, COMPUTER DOCUMENTATION UCID-30194,
C                 AUGUST 1982.
C***ROUTINES CALLED  PCHFD
C***END PROLOGUE  PCHEV
      INTEGER  N, NVAL, IERR
      REAL  X(N), F(N), D(N), XVAL(NVAL), FVAL(NVAL), DVAL(NVAL)
C
C  DECLARE LOCAL VARIABLES.
C
      INTEGER INCFD
      LOGICAL SKIP
      DATA SKIP /.TRUE./
      DATA INCFD /1/

C
C
C***FIRST EXECUTABLE STATEMENT  PCHEV
C
      CALL PCHFD(N,X,F,D,INCFD,SKIP,NVAL,XVAL,FVAL,DVAL,IERR)
C
C
 5000 CONTINUE
      RETURN
C
C------------- LAST LINE OF PCHEV FOLLOWS ------------------------------
      END
      SUBROUTINE PCHEZ(N,X,F,D,SPLINE,WK,LWK,IERR)
C***BEGIN PROLOGUE  PCHEZ
C***DATE WRITTEN   870821   (YYMMDD)
C***REVISION DATE  870908   (YYMMDD)
C***CATEGORY NO.  E1B
C***KEYWORDS  CUBIC HERMITE MONOTONE INTERPOLATION, SPLINE
C             INTERPOLATION, EASY TO USE PIECEWISE CUBIC INTERPOLATION
C***AUTHOR  KAHANER, D.K., (NBS)
C             SCIENTIFIC COMPUTING DIVISION
C             NATIONAL BUREAU OF STANDARDS
C             GAITHERSBURG, MARYLAND 20899
C             (301) 975-3808
C***PURPOSE  Easy to use spline or cubic Hermite interpolation.
C***DESCRIPTION
C
C          PCHEZ:  Piecewise Cubic Interpolation, Easy to Use.
C
C     From the book "Numerical Methods and Software"
C          by  D. Kahaner, C. Moler, S. Nash
C               Prentice Hall 1988
C
C     Sets derivatives for spline (two continuous derivatives) or
C     Hermite cubic (one continuous derivative) interpolation.
C     Spline interpolation is smoother, but may not "look" right if the
C     data contains both "steep" and "flat" sections.  Hermite cubics
C     can produce a "visually pleasing" and monotone interpolant to
C     monotone data. This is an easy to use driver for the routines
C     by F. N. Fritsch in reference (4) below. Various boundary
C     conditions are set to default values by PCHEZ. Many other choices
C     are available in the subroutines PCHIC, PCHIM and PCHSP.
C
C     Use PCHEV to evaluate the resulting function and its derivative.
C
C ----------------------------------------------------------------------
C
C  Calling sequence:   CALL  PCHEZ (N, X, F, D, SPLINE, WK, LWK, IERR)
C
C     INTEGER  N, IERR,  LWK
C     REAL  X(N), F(N), D(N), WK(*)
C     LOGICAL SPLINE
C
C   Parameters:
C
C     N -- (input) number of data points.  (Error return if N.LT.2 .)
C           If N=2, simply does linear interpolation.
C
C     X -- (input) real array of independent variable values.  The
C           elements of X must be strictly increasing:
C                X(I-1) .LT. X(I),  I = 2(1)N.
C           (Error return if not.)
C
C     F -- (input) real array of dependent variable values to be inter-
C           polated.  F(I) is value corresponding to X(I).
C
C     D -- (output) real array of derivative values at the data points.
C
C     SPLINE -- (input) logical variable to specify if the interpolant
C           is to be a spline with two continuous derivaties
C           (set SPLINE=.TRUE.) or a Hermite cubic interpolant with one
C           continuous derivative (set SPLINE=.FALSE.).
C        Note: If SPLINE=.TRUE. the interpolating spline satisfies the
C           default "not-a-knot" boundary condition, with a continuous
C           third derivative at X(2) and X(N-1). See reference (3).
C              If SPLINE=.FALSE. the interpolating Hermite cubic will be
C           monotone if the input data is monotone. Boundary conditions are
C           computed from the derivative of a local quadratic unless this
C           alters monotonicity.
C
C     WK -- (scratch) real work array, which must be declared by the calling
C           program to be at least 2*N if SPLINE is .TRUE. and not used
C           otherwise.
C
C     LWK -- (input) length of work array WK. (Error return if
C           LWK.LT.2*N and SPLINE is .TRUE., not checked otherwise.)
C
C     IERR -- (output) error flag.
C           Normal return:
C              IERR = 0  (no errors).
C           Warning error:
C              IERR.GT.0  (can only occur when SPLINE=.FALSE.) means that
C                 IERR switches in the direction of monotonicity were detected.
C                 When SPLINE=.FALSE.,  PCHEZ guarantees that if the input
C                 data is monotone, the interpolant will be too. This warning
C                 is to alert you to the fact that the input data was not
C                 monotone.
C           "Recoverable" errors:
C              IERR = -1  if N.LT.2 .
C              IERR = -3  if the X-array is not strictly increasing.
C              IERR = -7  if LWK is less than 2*N and SPLINE is .TRUE.
C             (The D-array has not been changed in any of these cases.)
C               NOTE:  The above errors are checked in the order listed,
C                   and following arguments have **NOT** been validated.
C
C ----------------------------------------------------------------------
C***REFERENCES  1. F.N.FRITSCH AND R.E.CARLSON, 'MONOTONE PIECEWISE
C                 CUBIC INTERPOLATION,' SIAM J.NUMER.ANAL. 17, 2 (APRIL
C                 1980), 238-246.
C               2. F.N.FRITSCH AND J.BUTLAND, 'A METHOD FOR CONSTRUCTING
C                 LOCAL MONOTONE PIECEWISE CUBIC INTERPOLANTS,' LLNL
C                 PREPRINT UCRL-87559 (APRIL 1982).
C               3. CARL DE BOOR, A PRACTICAL GUIDE TO SPLINES, SPRINGER-
C                 VERLAG (NEW YORK, 1978).  (ESP. CHAPTER IV, PP.49-62.)
C               4. F.N.FRITSCH, 'PIECEWISE CUBIC HERMITE INTERPOLATION
C                 PACKAGE, FINAL SPECIFICATIONS', LAWRENCE LIVERMORE
C                 NATIONAL LABORATORY, COMPUTER DOCUMENTATION UCID-30194,
C                 AUGUST 1982.
C***ROUTINES CALLED  PCHIM,PCHSP
C***END PROLOGUE  PCHEZ
      INTEGER  N, LWK, IERR
      REAL  X(N), F(N), D(N), WK(LWK)
      LOGICAL SPLINE
C
C  DECLARE LOCAL VARIABLES.
C
      INTEGER IC(2), INCFD
      REAL  VC(2)
      DATA IC(1) /0/
      DATA IC(2) /0/
      DATA INCFD /1/
C
C
C***FIRST EXECUTABLE STATEMENT  PCHEZ
C
      IF ( SPLINE ) THEN
        CALL  PCHSP (IC, VC, N, X, F, D, INCFD, WK, LWK, IERR)
      ELSE
        CALL  PCHIM (N, X, F, D, INCFD, IERR)
      ENDIF
C
C  ERROR CONDITIONS ALREADY CHECKED IN PCHSP OR PCHIM

      RETURN
C------------- LAST LINE OF PCHEZ FOLLOWS ------------------------------
      END
      SUBROUTINE PCHFD(N,X,F,D,INCFD,SKIP,NE,XE,FE,DE,IERR)
C***BEGIN PROLOGUE  PCHFD
C***DATE WRITTEN   811020   (YYMMDD)
C***REVISION DATE  870707   (YYMMDD)
C***CATEGORY NO.  E3,H1
C***KEYWORDS  LIBRARY=SLATEC(PCHIP),
C             TYPE=SINGLE PRECISION(PCHFD-S DPCHFD-D),
C             CUBIC HERMITE DIFFERENTIATION,CUBIC HERMITE EVALUATION,
C             HERMITE INTERPOLATION,PIECEWISE CUBIC EVALUATION
C***AUTHOR  FRITSCH, F. N., (LLNL)
C             MATHEMATICS AND STATISTICS DIVISION
C             LAWRENCE LIVERMORE NATIONAL LABORATORY
C             P.O. BOX 808  (L-316)
C             LIVERMORE, CA  94550
C             FTS 532-4275, (415) 422-4275
C***PURPOSE  Evaluate a piecewise cubic hermite function and its first
C            derivative at an array of points.  May be used by itself
C            for Hermite interpolation, or as an evaluator for PCHIM
C            or PCHIC.  If only function values are required, use
C            PCHFE instead.
C***DESCRIPTION
C
C          PCHFD:  Piecewise Cubic Hermite Function and Derivative
C                  evaluator
C
C     Evaluates the cubic Hermite function defined by  N, X, F, D,  to-
C     gether with its first derivative, at the points  XE(J), J=1(1)NE.
C
C     If only function values are required, use PCHFE, instead.
C
C     To provide compatibility with PCHIM and PCHIC, includes an
C     increment between successive values of the F- and D-arrays.
C
C ----------------------------------------------------------------------
C
C  Calling sequence:
C
C        PARAMETER  (INCFD = ...)
C        INTEGER  N, NE, IERR
C        REAL  X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE), DE(NE)
C        LOGICAL  SKIP
C
C        CALL  PCHFD (N, X, F, D, INCFD, SKIP, NE, XE, FE, DE, IERR)
C
C   Parameters:
C
C     N -- (input) number of data points.  (Error return if N.LT.2 .)
C
C     X -- (input) real array of independent variable values.  The
C           elements of X must be strictly increasing:
C                X(I-1) .LT. X(I),  I = 2(1)N.
C           (Error return if not.)
C
C     F -- (input) real array of function values.  F(1+(I-1)*INCFD) is
C           the value corresponding to X(I).
C
C     D -- (input) real array of derivative values.  D(1+(I-1)*INCFD) is
C           the value corresponding to X(I).
C
C     INCFD -- (input) increment between successive values in F and D.
C           (Error return if  INCFD.LT.1 .)
C
C     SKIP -- (input/output) logical variable which should be set to
C           .TRUE. if the user wishes to skip checks for validity of
C           preceding parameters, or to .FALSE. otherwise.
C           This will save time in case these checks have already
C           been performed (say, in PCHIM or PCHIC).
C           SKIP will be set to .TRUE. on normal return.
C
C     NE -- (input) number of evaluation points.  (Error return if
C           NE.LT.1 .)
C
C     XE -- (input) real array of points at which the functions are to
C           be evaluated.
C
C
C          NOTES:
C           1. The evaluation will be most efficient if the elements
C              of XE are increasing relative to X;
C              that is,   XE(J) .GE. X(I)
C              implies    XE(K) .GE. X(I),  all K.GE.J .
C           2. If any of the XE are outside the interval [X(1),X(N)],
C              values are extrapolated from the nearest extreme cubic,
C              and a warning error is returned.
C
C     FE -- (output) real array of values of the cubic Hermite function
C           defined by  N, X, F, D  at the points  XE.
C
C     DE -- (output) real array of values of the first derivative of
C           the same function at the points  XE.
C
C     IERR -- (output) error flag.
C           Normal return:
C              IERR = 0  (no errors).
C           Warning error:
C              IERR.GT.0  means that extrapolation was performed at
C                 IERR points.
C           "Recoverable" errors:
C              IERR = -1  if N.LT.2 .
C              IERR = -2  if INCFD.LT.1 .
C              IERR = -3  if the X-array is not strictly increasing.
C              IERR = -4  if NE.LT.1 .
C           (Output arrays have not been changed in any of these cases.)
C               NOTE:  The above errors are checked in the order listed,
C                   and following arguments have **NOT** been validated.
C              IERR = -5  if an error has occurred in the lower-level
C                         routine CHFDV.  NB: this should never happen.
C                         Notify the author **IMMEDIATELY** if it does.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  CHFDV,XERROR
C***END PROLOGUE  PCHFD
C
C ----------------------------------------------------------------------
C
C  Change record:
C     82-08-03   Minor cosmetic changes for release 1.
C     87-07-07   Minor cosmetic changes to prologue.
C
C ----------------------------------------------------------------------
C
C  Programming notes:
C
C     1. To produce a double precision version, simply:
C        a. Change PCHFD to DPCHFD, and CHFDV to DCHFDV, wherever they
C           occur,
C        b. Change the real declaration to double precision,
C
C     2. Most of the coding between the call to CHFDV and the end of
C        the IR-loop could be eliminated if it were permissible to
C        assume that XE is ordered relative to X.
C
C     3. CHFDV does not assume that X1 is less than X2.  thus, it would
C        be possible to write a version of PCHFD that assumes a strict-
C        ly decreasing X-array by simply running the IR-loop backwards
C        (and reversing the order of appropriate tests).
C
C     4. The present code has a minor bug, which I have decided is not
C        worth the effort that would be required to fix it.
C        If XE contains points in [X(N-1),X(N)], followed by points .LT.
C        X(N-1), followed by points .GT.X(N), the extrapolation points
C        will be counted (at least) twice in the total returned in IERR.
C
C  DECLARE ARGUMENTS.
C
      INTEGER  N, INCFD, NE, IERR
      REAL  X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE), DE(NE)
      LOGICAL  SKIP
C
C  DECLARE LOCAL VARIABLES.
C
      INTEGER  I, IERC, IR, J, JFIRST, NEXT(2), NJ
C
C  VALIDITY-CHECK ARGUMENTS.
C
C***FIRST EXECUTABLE STATEMENT  PCHFD
      IF (SKIP)  GO TO 5
C
      IF ( N.LT.2 )  GO TO 5001
      IF ( INCFD.LT.1 )  GO TO 5002
      DO 1  I = 2, N
         IF ( X(I).LE.X(I-1) )  GO TO 5003
    1 CONTINUE
C
C  FUNCTION DEFINITION IS OK, GO ON.
C
    5 CONTINUE
      IF ( NE.LT.1 )  GO TO 5004
      IERR = 0
      SKIP = .TRUE.
C
C  LOOP OVER INTERVALS.        (   INTERVAL INDEX IS  IL = IR-1  . )
C                              ( INTERVAL IS X(IL).LE.X.LT.X(IR) . )
      JFIRST = 1
      IR = 2
   10 CONTINUE
C
C     SKIP OUT OF LOOP IF HAVE PROCESSED ALL EVALUATION POINTS.
C
         IF (JFIRST .GT. NE)  GO TO 5000
C
C     LOCATE ALL POINTS IN INTERVAL.
C
         DO 20  J = JFIRST, NE
            IF (XE(J) .GE. X(IR))  GO TO 30
   20    CONTINUE
         J = NE + 1
         GO TO 40
C
C     HAVE LOCATED FIRST POINT BEYOND INTERVAL.
C
   30    CONTINUE
         IF (IR .EQ. N)  J = NE + 1
C
   40    CONTINUE
         NJ = J - JFIRST
C
C     SKIP EVALUATION IF NO POINTS IN INTERVAL.
C
         IF (NJ .EQ. 0)  GO TO 50
C
C     EVALUATE CUBIC AT XE(I),  I = JFIRST (1) J-1 .
C
C       ----------------------------------------------------------------
        CALL CHFDV (X(IR-1),X(IR), F(1,IR-1),F(1,IR), D(1,IR-1),D(1,IR),
     *              NJ, XE(JFIRST), FE(JFIRST), DE(JFIRST), NEXT, IERC)
C       ----------------------------------------------------------------
         IF (IERC .LT. 0)  GO TO 5005
C
         IF (NEXT(2) .EQ. 0)  GO TO 42
C        IF (NEXT(2) .GT. 0)  THEN
C           IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(2) TO THE
C           RIGHT OF X(IR).
C
            IF (IR .LT. N)  GO TO 41
C           IF (IR .EQ. N)  THEN
C              THESE ARE ACTUALLY EXTRAPOLATION POINTS.
               IERR = IERR + NEXT(2)
               GO TO 42
   41       CONTINUE
C           ELSE
C              WE SHOULD NEVER HAVE GOTTEN HERE.
               GO TO 5005
C           ENDIF
C        ENDIF
   42    CONTINUE
C
         IF (NEXT(1) .EQ. 0)  GO TO 49
C        IF (NEXT(1) .GT. 0)  THEN
C           IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(1) TO THE
C           LEFT OF X(IR-1).
C
            IF (IR .GT. 2)  GO TO 43
C           IF (IR .EQ. 2)  THEN
C              THESE ARE ACTUALLY EXTRAPOLATION POINTS.
               IERR = IERR + NEXT(1)
               GO TO 49
   43       CONTINUE
C           ELSE
C              XE IS NOT ORDERED RELATIVE TO X, SO MUST ADJUST
C              EVALUATION INTERVAL.
C
C              FIRST, LOCATE FIRST POINT TO LEFT OF X(IR-1).
               DO 44  I = JFIRST, J-1
                  IF (XE(I) .LT. X(IR-1))  GO TO 45
   44          CONTINUE
C              NOTE-- CANNOT DROP THROUGH HERE UNLESS THERE IS AN ERROR
C                     IN CHFDV.
               GO TO 5005
C
   45          CONTINUE
C              RESET J.  (THIS WILL BE THE NEW JFIRST.)
               J = I
C
C              NOW FIND OUT HOW FAR TO BACK UP IN THE X-ARRAY.
               DO 46  I = 1, IR-1
                  IF (XE(J) .LT. X(I)) GO TO 47
   46          CONTINUE
C              NB-- CAN NEVER DROP THROUGH HERE, SINCE XE(J).LT.X(IR-1).
C
   47          CONTINUE
C              AT THIS POINT, EITHER  XE(J) .LT. X(1)
C                 OR      X(I-1) .LE. XE(J) .LT. X(I) .
C              RESET IR, RECOGNIZING THAT IT WILL BE INCREMENTED BEFORE
C              CYCLING.
               IR = MAX0(1, I-1)
C           ENDIF
C        ENDIF
   49    CONTINUE
C
         JFIRST = J
C
C     END OF IR-LOOP.
C
   50 CONTINUE
      IR = IR + 1
      IF (IR .LE. N)  GO TO 10
C
C  NORMAL RETURN.
C
 5000 CONTINUE
      RETURN
C
C  ERROR RETURNS.
C
 5001 CONTINUE
C     N.LT.2 RETURN.
      IERR = -1
      CALL XERROR ('PCHFD -- NUMBER OF DATA POINTS LESS THAN TWO'
     *           , 44, IERR, 1)
      RETURN
C
 5002 CONTINUE
C     INCFD.LT.1 RETURN.
      IERR = -2
      CALL XERROR ('PCHFD -- INCREMENT LESS THAN ONE'
     *           , 32, IERR, 1)
      RETURN
C
 5003 CONTINUE
C     X-ARRAY NOT STRICTLY INCREASING.
      IERR = -3
      CALL XERROR ('PCHFD -- X-ARRAY NOT STRICTLY INCREASING'
     *           , 40, IERR, 1)
      RETURN
C
 5004 CONTINUE
C     NE.LT.1 RETURN.
      IERR = -4
      CALL XERROR ('PCHFD -- NUMBER OF EVALUATION POINTS LESS THAN ONE'
     *           , 50, IERR, 1)
      RETURN
C
 5005 CONTINUE
C     ERROR RETURN FROM CHFDV.
C   *** THIS CASE SHOULD NEVER OCCUR ***
      IERR = -5
      CALL XERROR ('PCHFD -- ERROR RETURN FROM CHFDV -- FATAL'
     *           , 41, IERR, 2)
      RETURN
C------------- LAST LINE OF PCHFD FOLLOWS ------------------------------
      END
      SUBROUTINE PCHFE(N,X,F,D,INCFD,SKIP,NE,XE,FE,IERR)
C***BEGIN PROLOGUE  PCHFE
C***DATE WRITTEN   811020   (YYMMDD)
C***REVISION DATE  870707   (YYMMDD)
C***CATEGORY NO.  E3
C***KEYWORDS  LIBRARY=SLATEC(PCHIP),
C             TYPE=SINGLE PRECISION(PCHFE-S DPCHFE-D),
C             CUBIC HERMITE EVALUATION,HERMITE INTERPOLATION,
C             PIECEWISE CUBIC EVALUATION
C***AUTHOR  FRITSCH, F. N., (LLNL)
C             MATHEMATICS AND STATISTICS DIVISION
C             LAWRENCE LIVERMORE NATIONAL LABORATORY
C             P.O. BOX 808  (L-316)
C             LIVERMORE, CA  94550
C             FTS 532-4275, (415) 422-4275
C***PURPOSE  Evaluate a piecewise cubic Hermite function at an array of
C            points.  May be used by itself for Hermite interpolation,
C            or as an evaluator for PCHIM or PCHIC.
C***DESCRIPTION
C
C          PCHFE:  Piecewise Cubic Hermite Function Evaluator
C
C     Evaluates the cubic Hermite function defined by  N, X, F, D  at
C     the points  XE(J), J=1(1)NE.
C
C     To provide compatibility with PCHIM and PCHIC, includes an
C     increment between successive values of the F- and D-arrays.
C
C ----------------------------------------------------------------------
C
C  Calling sequence:
C
C        PARAMETER  (INCFD = ...)
C        INTEGER  N, NE, IERR
C        REAL  X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE)
C        LOGICAL  SKIP
C
C        CALL  PCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)
C
C   Parameters:
C
C     N -- (input) number of data points.  (Error return if N.LT.2 .)
C
C     X -- (input) real array of independent variable values.  The
C           elements of X must be strictly increasing:
C                X(I-1) .LT. X(I),  I = 2(1)N.
C           (Error return if not.)
C
C     F -- (input) real array of function values.  F(1+(I-1)*INCFD) is
C           the value corresponding to X(I).
C
C     D -- (input) real array of derivative values.  D(1+(I-1)*INCFD) is
C           the value corresponding to X(I).
C
C     INCFD -- (input) increment between successive values in F and D.
C           (Error return if  INCFD.LT.1 .)
C
C     SKIP -- (input/output) logical variable which should be set to
C           .TRUE. if the user wishes to skip checks for validity of
C           preceding parameters, or to .FALSE. otherwise.
C           This will save time in case these checks have already
C           been performed (say, in PCHIM or PCHIC).
C           SKIP will be set to .TRUE. on normal return.
C
C     NE -- (input) number of evaluation points.  (Error return if
C           NE.LT.1 .)
C
C     XE -- (input) real array of points at which the function is to be
C           evaluated.
C
C          NOTES:
C           1. The evaluation will be most efficient if the elements
C              of XE are increasing relative to X;
C              that is,   XE(J) .GE. X(I)
C              implies    XE(K) .GE. X(I),  all K.GE.J .
C           2. If any of the XE are outside the interval [X(1),X(N)],
C              values are extrapolated from the nearest extreme cubic,
C              and a warning error is returned.
C
C     FE -- (output) real array of values of the cubic Hermite function
C           defined by  N, X, F, D  at the points  XE.
C
C     IERR -- (output) error flag.
C           Normal return:
C              IERR = 0  (no errors).
C           Warning error:
C              IERR.GT.0  means that extrapolation was performed at
C                 IERR points.
C           "Recoverable" errors:
C              IERR = -1  if N.LT.2 .
C              IERR = -2  if INCFD.LT.1 .
C              IERR = -3  if the X-array is not strictly increasing.
C              IERR = -4  if NE.LT.1 .
C             (The FE-array has not been changed in any of these cases.)
C               NOTE:  The above errors are checked in the order listed,
C                   and following arguments have **NOT** been validated.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  CHFEV,XERROR
C***END PROLOGUE  PCHFE
C
C ----------------------------------------------------------------------
C
C  Change record:
C     82-08-03   Minor cosmetic changes for release 1.
C     87-07-07   Minor cosmetic changes to prologue.
C
C ----------------------------------------------------------------------
C
C  Programming notes:
C
C     1. To produce a double precision version, simply:
C        a. Change PCHFE to DPCHFE, and CHFEV to DCHFEV, wherever they
C           occur,
C        b. Change the real declaration to double precision,
C
C     2. Most of the coding between the call to CHFEV and the end of
C        the IR-loop could be eliminated if it were permissible to
C        assume that XE is ordered relative to X.
C
C     3. CHFEV does not assume that X1 is less than X2.  thus, it would
C        be possible to write a version of PCHFE that assumes a strict-
C        ly decreasing X-array by simply running the IR-loop backwards
C        (and reversing the order of appropriate tests).
C
C     4. The present code has a minor bug, which I have decided is not
C        worth the effort that would be required to fix it.
C        If XE contains points in [X(N-1),X(N)], followed by points .LT.
C        X(N-1), followed by points .GT.X(N), the extrapolation points
C        will be counted (at least) twice in the total returned in IERR.
C
C  DECLARE ARGUMENTS.
C
      INTEGER  N, INCFD, NE, IERR
      REAL  X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE)
      LOGICAL  SKIP
C
C  DECLARE LOCAL VARIABLES.
C
      INTEGER  I, IERC, IR, J, JFIRST, NEXT(2), NJ
C
C  VALIDITY-CHECK ARGUMENTS.
C
C***FIRST EXECUTABLE STATEMENT  PCHFE
      IF (SKIP)  GO TO 5
C
      IF ( N.LT.2 )  GO TO 5001
      IF ( INCFD.LT.1 )  GO TO 5002
      DO 1  I = 2, N
         IF ( X(I).LE.X(I-1) )  GO TO 5003
    1 CONTINUE
C
C  FUNCTION DEFINITION IS OK, GO ON.
C
    5 CONTINUE
      IF ( NE.LT.1 )  GO TO 5004
      IERR = 0
      SKIP = .TRUE.
C
C  LOOP OVER INTERVALS.        (   INTERVAL INDEX IS  IL = IR-1  . )
C                              ( INTERVAL IS X(IL).LE.X.LT.X(IR) . )
      JFIRST = 1
      IR = 2
   10 CONTINUE
C
C     SKIP OUT OF LOOP IF HAVE PROCESSED ALL EVALUATION POINTS.
C
         IF (JFIRST .GT. NE)  GO TO 5000
C
C     LOCATE ALL POINTS IN INTERVAL.
C
         DO 20  J = JFIRST, NE
            IF (XE(J) .GE. X(IR))  GO TO 30
   20    CONTINUE
         J = NE + 1
         GO TO 40
C
C     HAVE LOCATED FIRST POINT BEYOND INTERVAL.
C
   30    CONTINUE
         IF (IR .EQ. N)  J = NE + 1
C
   40    CONTINUE
         NJ = J - JFIRST
C
C     SKIP EVALUATION IF NO POINTS IN INTERVAL.
C
         IF (NJ .EQ. 0)  GO TO 50
C
C     EVALUATE CUBIC AT XE(I),  I = JFIRST (1) J-1 .
C
C       ----------------------------------------------------------------
        CALL CHFEV (X(IR-1),X(IR), F(1,IR-1),F(1,IR), D(1,IR-1),D(1,IR),
     *              NJ, XE(JFIRST), FE(JFIRST), NEXT, IERC)
C       ----------------------------------------------------------------
         IF (IERC .LT. 0)  GO TO 5005
C
         IF (NEXT(2) .EQ. 0)  GO TO 42
C        IF (NEXT(2) .GT. 0)  THEN
C           IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(2) TO THE
C           RIGHT OF X(IR).
C
            IF (IR .LT. N)  GO TO 41
C           IF (IR .EQ. N)  THEN
C              THESE ARE ACTUALLY EXTRAPOLATION POINTS.
               IERR = IERR + NEXT(2)
               GO TO 42
   41       CONTINUE
C           ELSE
C              WE SHOULD NEVER HAVE GOTTEN HERE.
               GO TO 5005
C           ENDIF
C        ENDIF
   42    CONTINUE
C
         IF (NEXT(1) .EQ. 0)  GO TO 49
C        IF (NEXT(1) .GT. 0)  THEN
C           IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(1) TO THE
C           LEFT OF X(IR-1).
C
            IF (IR .GT. 2)  GO TO 43
C           IF (IR .EQ. 2)  THEN
C              THESE ARE ACTUALLY EXTRAPOLATION POINTS.
               IERR = IERR + NEXT(1)
               GO TO 49
   43       CONTINUE
C           ELSE
C              XE IS NOT ORDERED RELATIVE TO X, SO MUST ADJUST
C              EVALUATION INTERVAL.
C
C              FIRST, LOCATE FIRST POINT TO LEFT OF X(IR-1).
               DO 44  I = JFIRST, J-1
                  IF (XE(I) .LT. X(IR-1))  GO TO 45
   44          CONTINUE
C              NOTE-- CANNOT DROP THROUGH HERE UNLESS THERE IS AN ERROR
C                     IN CHFEV.
               GO TO 5005
C
   45          CONTINUE
C              RESET J.  (THIS WILL BE THE NEW JFIRST.)
               J = I
C
C              NOW FIND OUT HOW FAR TO BACK UP IN THE X-ARRAY.
               DO 46  I = 1, IR-1
                  IF (XE(J) .LT. X(I)) GO TO 47
   46          CONTINUE
C              NB-- CAN NEVER DROP THROUGH HERE, SINCE XE(J).LT.X(IR-1).
C
   47          CONTINUE
C              AT THIS POINT, EITHER  XE(J) .LT. X(1)
C                 OR      X(I-1) .LE. XE(J) .LT. X(I) .
C              RESET IR, RECOGNIZING THAT IT WILL BE INCREMENTED BEFORE
C              CYCLING.
               IR = MAX0(1, I-1)
C           ENDIF
C        ENDIF
   49    CONTINUE
C
         JFIRST = J
C
C     END OF IR-LOOP.
C
   50 CONTINUE
      IR = IR + 1
      IF (IR .LE. N)  GO TO 10
C
C  NORMAL RETURN.
C
 5000 CONTINUE
      RETURN
C
C  ERROR RETURNS.
C
 5001 CONTINUE
C     N.LT.2 RETURN.
      IERR = -1
      CALL XERROR ('PCHFE -- NUMBER OF DATA POINTS LESS THAN TWO'
     *           , 44, IERR, 1)
      RETURN
C
 5002 CONTINUE
C     INCFD.LT.1 RETURN.
      IERR = -2
      CALL XERROR ('PCHFE -- INCREMENT LESS THAN ONE'
     *           , 32, IERR, 1)
      RETURN
C
 5003 CONTINUE
C     X-ARRAY NOT STRICTLY INCREASING.
      IERR = -3
      CALL XERROR ('PCHFE -- X-ARRAY NOT STRICTLY INCREASING'
     *           , 40, IERR, 1)
      RETURN
C
 5004 CONTINUE
C     NE.LT.1 RETURN.
      IERR = -4
      CALL XERROR ('PCHFE -- NUMBER OF EVALUATION POINTS LESS THAN ONE'
     *           , 50, IERR, 1)
      RETURN
C
 5005 CONTINUE
C     ERROR RETURN FROM CHFEV.
C   *** THIS CASE SHOULD NEVER OCCUR ***
      IERR = -5
      CALL XERROR ('PCHFE -- ERROR RETURN FROM CHFEV -- FATAL'
     *           , 41, IERR, 2)
      RETURN
C------------- LAST LINE OF PCHFE FOLLOWS ------------------------------
      END
      REAL FUNCTION PCHIA(N,X,F,D,INCFD,SKIP,A,B,IERR)
C***BEGIN PROLOGUE  PCHIA
C***DATE WRITTEN   820730   (YYMMDD)
C***REVISION DATE  870707   (YYMMDD)
C***CATEGORY NO.  E3,H2A2
C***KEYWORDS  LIBRARY=SLATEC(PCHIP),
C             TYPE=SINGLE PRECISION(PCHIA-S DPCHIA-D),
C             CUBIC HERMITE INTERPOLATION,NUMERICAL INTEGRATION,
C             QUADRATURE
C***AUTHOR  FRITSCH, F. N., (LLNL)
C             MATHEMATICS AND STATISTICS DIVISION
C             LAWRENCE LIVERMORE NATIONAL LABORATORY
C             P.O. BOX 808  (L-316)
C             LIVERMORE, CA  94550
C             FTS 532-4275, (415) 422-4275
C***PURPOSE  Evaluate the definite integral of a piecewise cubic
C            Hermite function over an arbitrary interval.
C***DESCRIPTION
C
C          PCHIA:  Piecewise Cubic Hermite Integrator, Arbitrary limits
C
C     Evaluates the definite integral of the cubic Hermite function
C     defined by  N, X, F, D  over the interval [A, B].
C
C     To provide compatibility with PCHIM and PCHIC, includes an
C     increment between successive values of the F- and D-arrays.
C
C ----------------------------------------------------------------------
C
C  Calling sequence:
C
C        PARAMETER  (INCFD = ...)
C        INTEGER  N, IERR
C        REAL  X(N), F(INCFD,N), D(INCFD,N), A, B
C        LOGICAL  SKIP
C
C        VALUE = PCHIA (N, X, F, D, INCFD, SKIP, A, B, IERR)
C
C   Parameters:
C
C     VALUE -- (output) VALUE of the requested integral.
C
C     N -- (input) number of data points.  (Error return if N.LT.2 .)
C
C     X -- (input) real array of independent variable values.  The
C           elements of X must be strictly increasing:
C                X(I-1) .LT. X(I),  I = 2(1)N.
C           (Error return if not.)
C
C     F -- (input) real array of function values.  F(1+(I-1)*INCFD) is
C           the value corresponding to X(I).
C
C     D -- (input) real array of derivative values.  D(1+(I-1)*INCFD) is
C           the value corresponding to X(I).
C
C     INCFD -- (input) increment between successive values in F and D.
C           (Error return if  INCFD.LT.1 .)
C
C     SKIP -- (input/output) logical variable which should be set to
C           .TRUE. if the user wishes to skip checks for validity of
C           preceding parameters, or to .FALSE. otherwise.
C           This will save time in case these checks have already
C           been performed (say, in PCHIM or PCHIC).
C           SKIP will be set to .TRUE. on return with IERR.GE.0 .
C
C     A,B -- (input) the limits of integration.
C           NOTE:  There is no requirement that [A,B] be contained in
C                  [X(1),X(N)].  However, the resulting integral value
C                  will be highly suspect, if not.
C
C     IERR -- (output) error flag.
C           Normal return:
C              IERR = 0  (no errors).
C           Warning errors:
C              IERR = 1  if  A  is outside the interval [X(1),X(N)].
C              IERR = 2  if  B  is outside the interval [X(1),X(N)].
C              IERR = 3  if both of the above are true.  (Note that this
C                        means that either [A,B] contains data interval
C                        or the intervals do not intersect at all.)
C           "Recoverable" errors:
C              IERR = -1  if N.LT.2 .
C              IERR = -2  if INCFD.LT.1 .
C              IERR = -3  if the X-array is not strictly increasing.
C                (Value has not been computed in any of these cases.)
C               NOTE:  The above errors are checked in the order listed,
C                   and following arguments have **NOT** been validated.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  CHFIV,PCHID,XERROR
C***END PROLOGUE  PCHIA
C
C ----------------------------------------------------------------------
C
C
C  Change record:
C     82-08-04   Converted to SLATEC library version.
C     87-07-07   Corrected double precision conversion instructions.
C
C ----------------------------------------------------------------------
C
C  Programming notes:
C
C     To produce a double precision version, simply:
C        a. Change PCHIA to DPCHIA wherever it occurs,
C        b. Change PCHID to DPCHID wherever it occurs,
C        c. Change CHFIV to DCHFIV wherever it occurs,
C        d. Change the real declarations to double precision,  and
C        e. Change the constant  ZERO  to double precision.
C
C  DECLARE ARGUMENTS.
C
      INTEGER  N, INCFD, IERR
      REAL  X(N), F(INCFD,N), D(INCFD,N), A, B
      LOGICAL  SKIP
C
C  DECLARE LOCAL VARIABLES.
C
      INTEGER  I, IA, IB, IERD, IERV, IL, IR
      REAL  VALUE, XA, XB, ZERO
      REAL  CHFIV, PCHID
C
C  INITIALIZE.
C
      DATA  ZERO /0./
C
C  VALIDITY-CHECK ARGUMENTS.
C
C***FIRST EXECUTABLE STATEMENT  PCHIA
      IF (SKIP)  GO TO 5
C
      IF ( N.LT.2 )  GO TO 5001
      IF ( INCFD.LT.1 )  GO TO 5002
      DO 1  I = 2, N
         IF ( X(I).LE.X(I-1) )  GO TO 5003
    1 CONTINUE
C
C  FUNCTION DEFINITION IS OK, GO ON.
C
    5 CONTINUE
      SKIP = .TRUE.
      IERR = 0
      IF ( (A.LT.X(1)) .OR. (A.GT.X(N)) )  IERR = IERR + 1
      IF ( (B.LT.X(1)) .OR. (B.GT.X(N)) )  IERR = IERR + 2
C
C  COMPUTE INTEGRAL VALUE.
C
      IF (A .EQ. B)  THEN
         VALUE = ZERO
      ELSE
         XA = AMIN1 (A, B)
         XB = AMAX1 (A, B)
         IF (XB .LE. X(2))  THEN
C           INTERVAL IS TO LEFT OF X(2), SO USE FIRST CUBIC.
C                   --------------------------------------------
            VALUE = CHFIV (X(1),X(2), F(1,1),F(1,2),
     *                                D(1,1),D(1,2), A, B, IERV)
C                   --------------------------------------------
            IF (IERV .LT. 0)  GO TO 5004
         ELSE IF (XA .GE. X(N-1))  THEN
C           INTERVAL IS TO RIGHT OF X(N-1), SO USE LAST CUBIC.
C                   -----------------------------------------------
            VALUE = CHFIV(X(N-1),X(N), F(1,N-1),F(1,N),
     *                                 D(1,N-1),D(1,N), A, B, IERV)
C                   -----------------------------------------------
            IF (IERV .LT. 0)  GO TO 5004
         ELSE
C           'NORMAL' CASE -- XA.LT.XB, XA.LT.X(N-1), XB.GT.X(2).
C      ......LOCATE IA AND IB SUCH THAT
C               X(IA-1).LT.XA.LE.X(IA).LE.X(IB).LE.XB.LE.X(IB+1)
            IA = 1
            DO 10  I = 1, N-1
               IF (XA .GT. X(I))  IA = I + 1
   10       CONTINUE
C             IA = 1 IMPLIES XA.LT.X(1) .  OTHERWISE,
C             IA IS LARGEST INDEX SUCH THAT X(IA-1).LT.XA,.
C
            IB = N
            DO 20  I = N, IA, -1
               IF (XB .LT. X(I))  IB = I - 1
   20       CONTINUE
C             IB = N IMPLIES XB.GT.X(N) .  OTHERWISE,
C             IB IS SMALLEST INDEX SUCH THAT XB.LT.X(IB+1) .
C
C     ......COMPUTE THE INTEGRAL.
            IERV = 0
            IF (IB .LT. IA)  THEN
C              THIS MEANS IB = IA-1 AND
C                 (A,B) IS A SUBSET OF (X(IB),X(IA)).
C                      ------------------------------------------------
               VALUE = CHFIV (X(IB),X(IA), F(1,IB),F(1,IA),
     *                                     D(1,IB),D(1,IA), A, B, IERV)
C                      ------------------------------------------------
               IF (IERV .LT. 0)  GO TO 5004
            ELSE
C
C              FIRST COMPUTE INTEGRAL OVER (X(IA),X(IB)).
               IF (IB .EQ. IA)  THEN
                  VALUE = ZERO
               ELSE
C                         ---------------------------------------------
                  VALUE = PCHID (N, X, F, D, INCFD, SKIP, IA, IB, IERD)
C                         ---------------------------------------------
                  IF (IERD .LT. 0)  GO TO 5005
               ENDIF
C
C              THEN ADD ON INTEGRAL OVER (XA,X(IA)).
               IF (XA .LT. X(IA))  THEN
                  IL = MAX0 (1, IA-1)
                  IR = IL + 1
C                                 -------------------------------------
                  VALUE = VALUE + CHFIV (X(IL),X(IR), F(1,IL),F(1,IR),
     *                                D(1,IL),D(1,IR), XA, X(IA), IERV)
C                                 -------------------------------------
                  IF (IERV .LT. 0)  GO TO 5004
               ENDIF
C
C              THEN ADD ON INTEGRAL OVER (X(IB),XB).
               IF (XB .GT. X(IB))  THEN
                  IR = MIN0 (IB+1, N)
                  IL = IR - 1
C                                 -------------------------------------
                  VALUE = VALUE + CHFIV (X(IL),X(IR), F(1,IL),F(1,IR),
     *                                D(1,IL),D(1,IR), X(IB), XB, IERV)
C                                 -------------------------------------
                  IF (IERV .LT. 0)  GO TO 5004
               ENDIF
C
C              FINALLY, ADJUST SIGN IF NECESSARY.
               IF (A .GT. B)  VALUE = -VALUE
            ENDIF
         ENDIF
      ENDIF
C
C  NORMAL RETURN.
C
      PCHIA = VALUE
      RETURN
C
C  ERROR RETURNS.
C
 5001 CONTINUE
C     N.LT.2 RETURN.
      IERR = -1
      CALL XERROR ('PCHIA -- NUMBER OF DATA POINTS LESS THAN TWO'
     *           , 44, IERR, 1)
      RETURN
C
 5002 CONTINUE
C     INCFD.LT.1 RETURN.
      IERR = -2
      CALL XERROR ('PCHIA -- INCREMENT LESS THAN ONE'
     *           , 32, IERR, 1)
      RETURN
C
 5003 CONTINUE
C     X-ARRAY NOT STRICTLY INCREASING.
      IERR = -3
      CALL XERROR ('PCHIA -- X-ARRAY NOT STRICTLY INCREASING'
     *           , 40, IERR, 1)
      RETURN
C
 5004 CONTINUE
C     TROUBLE IN CHFIV.  (SHOULD NEVER OCCUR.)
      IERR = -4
      CALL XERROR ('PCHIA -- TROUBLE IN CHFIV'
     *           , 25, IERR, 1)
      RETURN
C
 5005 CONTINUE
C     TROUBLE IN PCHID.  (SHOULD NEVER OCCUR.)
      IERR = -5
      CALL XERROR ('PCHIA -- TROUBLE IN PCHID'
     *           , 25, IERR, 1)
      RETURN
C------------- LAST LINE OF PCHIA FOLLOWS ------------------------------
      END
      SUBROUTINE PCHIC(IC,VC,SWITCH,N,X,F,D,INCFD,WK,NWK,IERR)
C***BEGIN PROLOGUE  PCHIC
C***DATE WRITTEN   820218   (YYMMDD)
C***REVISION DATE  870707   (YYMMDD)
C***CATEGORY NO.  E1B
C***KEYWORDS  LIBRARY=SLATEC(PCHIP),
C             TYPE=SINGLE PRECISION(PCHIC-S DPCHIC-D),
C             CUBIC HERMITE INTERPOLATION,MONOTONE INTERPOLATION,
C             PIECEWISE CUBIC INTERPOLATION,
C             SHAPE-PRESERVING INTERPOLATION
C***AUTHOR  FRITSCH, F. N., (LLNL)
C             MATHEMATICS AND STATISTICS DIVISION
C             LAWRENCE LIVERMORE NATIONAL LABORATORY
C             P.O. BOX 808  (L-316)
C             LIVERMORE, CA  94550
C             FTS 532-4275, (415) 422-4275
C***PURPOSE  Set derivatives needed to determine a piecewise monotone
C            piecewise cubic Hermite interpolant to given data.
C            User control is available over boundary conditions and/or
C            treatment of points where monotonicity switches direction.
C***DESCRIPTION
C
C         PCHIC:  Piecewise Cubic Hermite Interpolation Coefficients.
C
C     Sets derivatives needed to determine a piecewise monotone piece-
C     wise cubic interpolant to the data given in X and F satisfying the
C     boundary conditions specified by IC and VC.
C
C     The treatment of points where monotonicity switches direction is
C     controlled by argument SWITCH.
C
C     To facilitate two-dimensional applications, includes an increment
C     between successive values of the F- and D-arrays.
C
C     The resulting piecewise cubic Hermite function may be evaluated
C     by PCHFE or PCHFD.
C
C ----------------------------------------------------------------------
C
C  Calling sequence:
C
C        PARAMETER  (INCFD = ...)
C        INTEGER  IC(2), N, NWK, IERR
C        REAL  VC(2), SWITCH, X(N), F(INCFD,N), D(INCFD,N), WK(NWK)
C
C        CALL  PCHIC (IC, VC, SWITCH, N, X, F, D, INCFD, WK, NWK, IERR)
C
C   Parameters:
C
C     IC -- (input) integer array of length 2 specifying desired
C           boundary conditions:
C           IC(1) = IBEG, desired condition at beginning of data.
C           IC(2) = IEND, desired condition at end of data.
C
C           IBEG = 0  for the default boundary condition (the same as
C                     used by PCHIM).
C           If IBEG.NE.0, then its sign indicates whether the boundary
C                     derivative is to be adjusted, if necessary, to be
C                     compatible with monotonicity:
C              IBEG.GT.0  if no adjustment is to be performed.
C              IBEG.LT.0  if the derivative is to be adjusted for
C                     monotonicity.
C
C           Allowable values for the magnitude of IBEG are:
C           IBEG = 1  if first derivative at X(1) is given in VC(1).
C           IBEG = 2  if second derivative at X(1) is given in VC(1).
C           IBEG = 3  to use the 3-point difference formula for D(1).
C                     (Reverts to the default b.c. if N.LT.3 .)
C           IBEG = 4  to use the 4-point difference formula for D(1).
C                     (Reverts to the default b.c. if N.LT.4 .)
C           IBEG = 5  to set D(1) so that the second derivative is con-
C              tinuous at X(2). (Reverts to the default b.c. if N.LT.4.)
C              This option is somewhat analogous to the "not a knot"
C              boundary condition provided by PCHSP.
C
C          NOTES (IBEG):
C           1. An error return is taken if ABS(IBEG).GT.5 .
C           2. Only in case  IBEG.LE.0  is it guaranteed that the
C              interpolant will be monotonic in the first interval.
C              If the returned value of D(1) lies between zero and
C              3*SLOPE(1), the interpolant will be monotonic.  This
C              is **NOT** checked if IBEG.GT.0 .
C           3. If IBEG.LT.0 and D(1) had to be changed to achieve mono-
C              tonicity, a warning error is returned.
C
C           IEND may take on the same values as IBEG, but applied to
C           derivative at X(N).  In case IEND = 1 or 2, the value is
C           given in VC(2).
C
C          NOTES (IEND):
C           1. An error return is taken if ABS(IEND).GT.5 .
C           2. Only in case  IEND.LE.0  is it guaranteed that the
C              interpolant will be monotonic in the last interval.
C              If the returned value of D(1+(N-1)*INCFD) lies between
C              zero and 3*SLOPE(N-1), the interpolant will be monotonic.
C              This is **NOT** checked if IEND.GT.0 .
C           3. If IEND.LT.0 and D(1+(N-1)*INCFD) had to be changed to
C              achieve monotonicity, a warning error is returned.
C
C     VC -- (input) real array of length 2 specifying desired boundary
C           values, as indicated above.
C           VC(1) need be set only if IC(1) = 1 or 2 .
C           VC(2) need be set only if IC(2) = 1 or 2 .
C
C     SWITCH -- (input) indicates desired treatment of points where
C           direction of monotonicity switches:
C           Set SWITCH to zero if interpolant is required to be mono-
C           tonic in each interval, regardless of monotonicity of data.
C             NOTES:
C              1. This will cause D to be set to zero at all switch
C                 points, thus forcing extrema there.
C              2. The result of using this option with the default boun-
C                 dary conditions will be identical to using PCHIM, but
C                 will generally cost more compute time.
C                 This option is provided only to facilitate comparison
C                 of different switch and/or boundary conditions.
C           Set SWITCH nonzero to use a formula based on the 3-point
C              difference formula in the vicinity of switch points.
C           If SWITCH is positive, the interpolant on each interval
C              containing an extremum is controlled to not deviate from
C              the data by more than SWITCH*DFLOC, where DFLOC is the
C              maximum of the change of F on this interval and its two
C              immediate neighbors.
C           If SWITCH is negative, no such control is to be imposed.
C
C     N -- (input) number of data points.  (Error return if N.LT.2 .)
C
C     X -- (input) real array of independent variable values.  The
C           elements of X must be strictly increasing:
C                X(I-1) .LT. X(I),  I = 2(1)N.
C           (Error return if not.)
C
C     F -- (input) real array of dependent variable values to be inter-
C           polated.  F(1+(I-1)*INCFD) is value corresponding to X(I).
C
C     D -- (output) real array of derivative values at the data points.
C           These values will determine a monotone cubic Hermite func-
C           tion on each subinterval on which the data are monotonic,
C           except possibly adjacent to switches in monotonicity.
C           The value corresponding to X(I) is stored in
C                D(1+(I-1)*INCFD),  I=1(1)N.
C           No other entries in D are changed.
C
C     INCFD -- (input) increment between successive values in F and D.
C           This argument is provided primarily for 2-D applications.
C           (Error return if  INCFD.LT.1 .)
C
C     WK -- (scratch) real array of working storage.  The user may wish
C           to know that the returned values are:
C              WK(I)     = H(I)     = X(I+1) - X(I) ;
C              WK(N-1+I) = SLOPE(I) = (F(1,I+1) - F(1,I)) / H(I)
C           for  I = 1(1)N-1.
C
C     NWK -- (input) length of work array.
C           (Error return if  NWK.LT.2*(N-1) .)
C
C     IERR -- (output) error flag.
C           Normal return:
C              IERR = 0  (no errors).
C           Warning errors:
C              IERR = 1  if IBEG.LT.0 and D(1) had to be adjusted for
C                        monotonicity.
C              IERR = 2  if IEND.LT.0 and D(1+(N-1)*INCFD) had to be
C                        adjusted for monotonicity.
C              IERR = 3  if both of the above are true.
C           "Recoverable" errors:
C              IERR = -1  if N.LT.2 .
C              IERR = -2  if INCFD.LT.1 .
C              IERR = -3  if the X-array is not strictly increasing.
C              IERR = -4  if ABS(IBEG).GT.5 .
C              IERR = -5  if ABS(IEND).GT.5 .
C              IERR = -6  if both of the above are true.
C              IERR = -7  if NWK.LT.2*(N-1) .
C             (The D-array has not been changed in any of these cases.)
C               NOTE:  The above errors are checked in the order listed,
C                   and following arguments have **NOT** been validated.
C
C***REFERENCES  1. F.N.FRITSCH AND R.E.CARLSON, 'MONOTONE PIECEWISE
C                 CUBIC INTERPOLATION,' SIAM J.NUMER.ANAL. 17, 2 (APRIL
C                 1980), 238-246.
C               2. F.N.FRITSCH AND J.BUTLAND, 'A METHOD FOR CONSTRUCTING
C                 LOCAL MONOTONE PIECEWISE CUBIC INTERPOLANTS,' LLNL
C                 PREPRINT UCRL-87559 (APRIL 1982).
C               3. F.N.FRITSCH, 'PIECEWISE CUBIC INTERPOLATION PACKAGE,'
C                 LLNL PREPRINT UCRL-87285 (JULY 1982).
C***ROUTINES CALLED  PCHCE,PCHCI,PCHCS,XERROR
C***END PROLOGUE  PCHIC
C
C ----------------------------------------------------------------------
C
C  Change record:
C     82-08-04   Converted to SLATEC library version.
C
C ----------------------------------------------------------------------
C
C  Programming notes:
C
C     To produce a double precision version, simply:
C        a. Change PCHIC to DPCHIC wherever it occurs,
C        b. Change PCHCE to DPCHCE wherever it occurs,
C        c. Change PCHCI to DPCHCI wherever it occurs,
C        d. Change PCHCS to DPCHCS wherever it occurs,
C        e. Change the real declarations to double precision, and
C        f. Change the constant  ZERO  to double precision.
C
C  DECLARE ARGUMENTS.
C
      INTEGER  IC(2), N, INCFD, NWK, IERR
      REAL  VC(2), SWITCH, X(N), F(INCFD,N), D(INCFD,N), WK(NWK)
C
C  DECLARE LOCAL VARIABLES.
C
      INTEGER  I, IBEG, IEND, NLESS1
      REAL  ZERO
      DATA  ZERO /0./
C
C  VALIDITY-CHECK ARGUMENTS.
C
C***FIRST EXECUTABLE STATEMENT  PCHIC
      IF ( N.LT.2 )  GO TO 5001
      IF ( INCFD.LT.1 )  GO TO 5002
      DO 1  I = 2, N
         IF ( X(I).LE.X(I-1) )  GO TO 5003
    1 CONTINUE
C
      IBEG = IC(1)
      IEND = IC(2)
      IERR = 0
      IF (IABS(IBEG) .GT. 5)  IERR = IERR - 1
      IF (IABS(IEND) .GT. 5)  IERR = IERR - 2
      IF (IERR .LT. 0)  GO TO 5004
C
C  FUNCTION DEFINITION IS OK -- GO ON.
C
      NLESS1 = N - 1
      IF ( NWK .LT. 2*NLESS1 )  GO TO 5007
C
C  SET UP H AND SLOPE ARRAYS.
C
      DO 20  I = 1, NLESS1
         WK(I) = X(I+1) - X(I)
         WK(NLESS1+I) = (F(1,I+1) - F(1,I)) / WK(I)
   20 CONTINUE
C
C  SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION.
C
      IF (NLESS1 .GT. 1)  GO TO 1000
      D(1,1) = WK(2)
      D(1,N) = WK(2)
      GO TO 3000
C
C  NORMAL CASE  (N .GE. 3) .
C
 1000 CONTINUE
C
C  SET INTERIOR DERIVATIVES AND DEFAULT END CONDITIONS.
C
C     --------------------------------------
      CALL PCHCI (N, WK(1), WK(N), D, INCFD)
C     --------------------------------------
C
C  SET DERIVATIVES AT POINTS WHERE MONOTONICITY SWITCHES DIRECTION.
C
      IF (SWITCH .EQ. ZERO)  GO TO 3000
C     ----------------------------------------------------
      CALL PCHCS (SWITCH, N, WK(1), WK(N), D, INCFD, IERR)
C     ----------------------------------------------------
      IF (IERR .NE. 0)  GO TO 5008
C
C  SET END CONDITIONS.
C
 3000 CONTINUE
      IF ( (IBEG.EQ.0) .AND. (IEND.EQ.0) )  GO TO 5000
C     -------------------------------------------------------
      CALL PCHCE (IC, VC, N, X, WK(1), WK(N), D, INCFD, IERR)
C     -------------------------------------------------------
      IF (IERR .LT. 0)  GO TO 5009
C
C  NORMAL RETURN.
C
 5000 CONTINUE
      RETURN
C
C  ERROR RETURNS.
C
 5001 CONTINUE
C     N.LT.2 RETURN.
      IERR = -1
      CALL XERROR ('PCHIC -- NUMBER OF DATA POINTS LESS THAN TWO'
     *           , 44, IERR, 1)
      RETURN
C
 5002 CONTINUE
C     INCFD.LT.1 RETURN.
      IERR = -2
      CALL XERROR ('PCHIC -- INCREMENT LESS THAN ONE'
     *           , 32, IERR, 1)
      RETURN
C
 5003 CONTINUE
C     X-ARRAY NOT STRICTLY INCREASING.
      IERR = -3
      CALL XERROR ('PCHIC -- X-ARRAY NOT STRICTLY INCREASING'
     *           , 40, IERR, 1)
      RETURN
C
 5004 CONTINUE
C     IC OUT OF RANGE RETURN.
      IERR = IERR - 3
      CALL XERROR ('PCHIC -- IC OUT OF RANGE'
     *           , 24, IERR, 1)
      RETURN
C
 5007 CONTINUE
C     NWK .LT. 2*(N-1)  RETURN.
      IERR = -7
      CALL XERROR ('PCHIC -- WORK ARRAY TOO SMALL'
     *           , 29, IERR, 1)
      RETURN
C
 5008 CONTINUE
C     ERROR RETURN FROM PCHCS.
      IERR = -8
      CALL XERROR ('PCHIC -- ERROR RETURN FROM PCHCS'
     *           , 32, IERR, 1)
      RETURN
C
 5009 CONTINUE
C     ERROR RETURN FROM PCHCE.
C   *** THIS CASE SHOULD NEVER OCCUR ***
      IERR = -9
      CALL XERROR ('PCHIC -- ERROR RETURN FROM PCHCE'
     *           , 32, IERR, 1)
      RETURN
C------------- LAST LINE OF PCHIC FOLLOWS ------------------------------
      END
      REAL FUNCTION PCHID(N,X,F,D,INCFD,SKIP,IA,IB,IERR)
C***BEGIN PROLOGUE  PCHID
C***DATE WRITTEN   820723   (YYMMDD)
C***REVISION DATE  870707   (YYMMDD)
C***CATEGORY NO.  E1B,H2A2
C***KEYWORDS  LIBRARY=SLATEC(PCHIP),
C             TYPE=SINGLE PRECISION(PCHID-S DPCHID-D),
C             CUBIC HERMITE INTERPOLATION,NUMERICAL INTEGRATION,
C             QUADRATURE
C***AUTHOR  FRITSCH, F. N., (LLNL)
C             MATHEMATICS AND STATISTICS DIVISION
C             LAWRENCE LIVERMORE NATIONAL LABORATORY
C             P.O. BOX 808  (L-316)
C             LIVERMORE, CA  94550
C             FTS 532-4275, (415) 422-4275
C***PURPOSE  Evaluate the definite integral of a piecewise cubic
C            Hermite function over an interval whose endpoints are
C            data points.
C***DESCRIPTION
C
C          PCHID:  Piecewise Cubic Hermite Integrator, Data limits
C
C     Evaluates the definite integral of the cubic Hermite function
C     defined by  N, X, F, D  over the interval [X(IA), X(IB)].
C
C     To provide compatibility with PCHIM and PCHIC, includes an
C     increment between successive values of the F- and D-arrays.
C
C ----------------------------------------------------------------------
C
C  Calling sequence:
C
C        PARAMETER  (INCFD = ...)
C        INTEGER  N, IA, IB, IERR
C        REAL  X(N), F(INCFD,N), D(INCFD,N)
C        LOGICAL  SKIP
C
C        VALUE = PCHID (N, X, F, D, INCFD, SKIP, IA, IB, IERR)
C
C   Parameters:
C
C     VALUE -- (output) VALUE of the requested integral.
C
C     N -- (input) number of data points.  (Error return if N.LT.2 .)
C
C     X -- (input) real array of independent variable values.  The
C           elements of X must be strictly increasing:
C                X(I-1) .LT. X(I),  I = 2(1)N.
C           (Error return if not.)
C
C     F -- (input) real array of function values.  F(1+(I-1)*INCFD) is
C           the value corresponding to X(I).
C
C     D -- (input) real array of derivative values.  D(1+(I-1)*INCFD) is
C           the value corresponding to X(I).
C
C     INCFD -- (input) increment between successive values in F and D.
C           (Error return if  INCFD.LT.1 .)
C
C     SKIP -- (input/output) logical variable which should be set to
C           .TRUE. if the user wishes to skip checks for validity of
C           preceding parameters, or to .FALSE. otherwise.
C           This will save time in case these checks have already
C           been performed (say, in PCHIM or PCHIC).
C           SKIP will be set to .TRUE. on return with IERR = 0 or -4.
C
C     IA,IB -- (input) indices in X-array for the limits of integration.
C           both must be in the range [1,N].  (Error return if not.)
C           No restrictions on their relative values.
C
C     IERR -- (output) error flag.
C           Normal return:
C              IERR = 0  (no errors).
C           "Recoverable" errors:
C              IERR = -1  if N.LT.2 .
C              IERR = -2  if INCFD.LT.1 .
C              IERR = -3  if the X-array is not strictly increasing.
C              IERR = -4  if IA or IB is out of range.
C                (Value has not been computed in any of these cases.)
C               NOTE:  The above errors are checked in the order listed,
C                   and following arguments have **NOT** been validated.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  XERROR
C***END PROLOGUE  PCHID
C
C ----------------------------------------------------------------------
C
C  Change record:
C     82-08-04   Converted to SLATEC library version.
C
C ----------------------------------------------------------------------
C
C  Programming notes:
C
C     To produce a double precision version, simply:
C        a. Change PCHID to DPCHID wherever it occurs,
C        b. Change the real declarations to double precision,  and
C        c. Change the constants ZERO, HALF, SIX to double precision.
C
C  DECLARE ARGUMENTS.
C
      INTEGER  N, INCFD, IA, IB, IERR
      REAL  X(N), F(INCFD,N), D(INCFD,N)
      LOGICAL  SKIP
C
C  DECLARE LOCAL VARIABLES.
C
      INTEGER  I, IUP, LOW
      REAL  H, HALF, SIX, SUM, VALUE, ZERO
C
C  INITIALIZE.
C
      DATA  ZERO /0./,  HALF /0.5/,  SIX /6./
C
C  VALIDITY-CHECK ARGUMENTS.
C
C***FIRST EXECUTABLE STATEMENT  PCHID
      IF (SKIP)  GO TO 5
C
      IF ( N.LT.2 )  GO TO 5001
      IF ( INCFD.LT.1 )  GO TO 5002
      DO 1  I = 2, N
         IF ( X(I).LE.X(I-1) )  GO TO 5003
    1 CONTINUE
C
C  FUNCTION DEFINITION IS OK, GO ON.
C
    5 CONTINUE
      SKIP = .TRUE.
      IF ((IA.LT.1) .OR. (IA.GT.N))  GO TO 5004
      IF ((IB.LT.1) .OR. (IB.GT.N))  GO TO 5004
      IERR = 0
C
C  COMPUTE INTEGRAL VALUE.
C
      IF (IA .EQ. IB)  THEN
         VALUE = ZERO
      ELSE
         LOW = MIN0(IA, IB)
         IUP = MAX0(IA, IB) - 1
         SUM = ZERO
         DO 10  I = LOW, IUP
            H = X(I+1) - X(I)
            SUM = SUM + H*( (F(1,I) + F(1,I+1)) +
     *                      (D(1,I) - D(1,I+1))*(H/SIX) )
   10    CONTINUE
         VALUE = HALF * SUM
         IF (IA .GT. IB)  VALUE = -VALUE
      ENDIF
C
C  NORMAL RETURN.
C
      PCHID = VALUE
      RETURN
C
C  ERROR RETURNS.
C
 5001 CONTINUE
C     N.LT.2 RETURN.
      IERR = -1
      CALL XERROR ('PCHID -- NUMBER OF DATA POINTS LESS THAN TWO'
     *           , 44, IERR, 1)
      RETURN
C
 5002 CONTINUE
C     INCFD.LT.1 RETURN.
      IERR = -2
      CALL XERROR ('PCHID -- INCREMENT LESS THAN ONE'
     *           , 32, IERR, 1)
      RETURN
C
 5003 CONTINUE
C     X-ARRAY NOT STRICTLY INCREASING.
      IERR = -3
      CALL XERROR ('PCHID -- X-ARRAY NOT STRICTLY INCREASING'
     *           , 40, IERR, 1)
      RETURN
C
 5004 CONTINUE
C     IA OR IB OUT OF RANGE RETURN.
      IERR = -4
      CALL XERROR ('PCHID -- IA OR IB OUT OF RANGE'
     *           , 30, IERR, 1)
      RETURN
C------------- LAST LINE OF PCHID FOLLOWS ------------------------------
      END
      SUBROUTINE PCHIM(N,X,F,D,INCFD,IERR)
C***BEGIN PROLOGUE  PCHIM
C***DATE WRITTEN   811103   (YYMMDD)
C***REVISION DATE  870707   (YYMMDD)
C***CATEGORY NO.  E1B
C***KEYWORDS  LIBRARY=SLATEC(PCHIP),
C             TYPE=SINGLE PRECISION(PCHIM-S DPCHIM-D),
C             CUBIC HERMITE INTERPOLATION,MONOTONE INTERPOLATION,
C             PIECEWISE CUBIC INTERPOLATION
C***AUTHOR  FRITSCH, F. N., (LLNL)
C             MATHEMATICS AND STATISTICS DIVISION
C             LAWRENCE LIVERMORE NATIONAL LABORATORY
C             P.O. BOX 808  (L-316)
C             LIVERMORE, CA  94550
C             FTS 532-4275, (415) 422-4275
C***PURPOSE  Set derivatives needed to determine a monotone piecewise
C            cubic Hermite interpolant to given data.  Boundary values
C            are provided which are compatible with monotonicity.  The
C            interpolant will have an extremum at each point where mono-
C            tonicity switches direction.  (See PCHIC if user control is
C            desired over boundary or switch conditions.)
C***DESCRIPTION
C
C          PCHIM:  Piecewise Cubic Hermite Interpolation to
C                  Monotone data.
C
C     Sets derivatives needed to determine a monotone piecewise cubic
C     Hermite interpolant to the data given in X and F.
C
C     Default boundary conditions are provided which are compatible
C     with monotonicity.  (See PCHIC if user control of boundary con-
C     ditions is desired.)
C
C     If the data are only piecewise monotonic, the interpolant will
C     have an extremum at each point where monotonicity switches direc-
C     tion.  (See PCHIC if user control is desired in such cases.)
C
C     To facilitate two-dimensional applications, includes an increment
C     between successive values of the F- and D-arrays.
C
C     The resulting piecewise cubic Hermite function may be evaluated
C     by PCHFE or PCHFD.
C
C ----------------------------------------------------------------------
C
C  Calling sequence:
C
C        PARAMETER  (INCFD = ...)
C        INTEGER  N, IERR
C        REAL  X(N), F(INCFD,N), D(INCFD,N)
C
C        CALL  PCHIM (N, X, F, D, INCFD, IERR)
C
C   Parameters:
C
C     N -- (input) number of data points.  (Error return if N.LT.2 .)
C           If N=2, simply does linear interpolation.
C
C     X -- (input) real array of independent variable values.  The
C           elements of X must be strictly increasing:
C                X(I-1) .LT. X(I),  I = 2(1)N.
C           (Error return if not.)
C
C     F -- (input) real array of dependent variable values to be inter-
C           polated.  F(1+(I-1)*INCFD) is value corresponding to X(I).
C           PCHIM is designed for monotonic data, but it will work for
C           any F-array.  It will force extrema at points where mono-
C           tonicity switches direction.  If some other treatment of
C           switch points is desired, PCHIC should be used instead.
C                                     -----
C     D -- (output) real array of derivative values at the data points.
C           If the data are monotonic, these values will determine a
C           a monotone cubic Hermite function.
C           The value corresponding to X(I) is stored in
C                D(1+(I-1)*INCFD),  I=1(1)N.
C           No other entries in D are changed.
C
C     INCFD -- (input) increment between successive values in F and D.
C           This argument is provided primarily for 2-D applications.
C           (Error return if  INCFD.LT.1 .)
C
C     IERR -- (output) error flag.
C           Normal return:
C              IERR = 0  (no errors).
C           Warning error:
C              IERR.GT.0  means that IERR switches in the direction
C                 of monotonicity were detected.
C           "Recoverable" errors:
C              IERR = -1  if N.LT.2 .
C              IERR = -2  if INCFD.LT.1 .
C              IERR = -3  if the X-array is not strictly increasing.
C             (The D-array has not been changed in any of these cases.)
C               NOTE:  The above errors are checked in the order listed,
C                   and following arguments have **NOT** been validated.
C
C***REFERENCES  1. F.N.FRITSCH AND R.E.CARLSON, 'MONOTONE PIECEWISE
C                 CUBIC INTERPOLATION,' SIAM J.NUMER.ANAL. 17, 2 (APRIL
C                 1980), 238-246.
C               2. F.N.FRITSCH AND J.BUTLAND, 'A METHOD FOR CONSTRUCTING
C                 LOCAL MONOTONE PIECEWISE CUBIC INTERPOLANTS,' LLNL
C                 PREPRINT UCRL-87559 (APRIL 1982).
C***ROUTINES CALLED  PCHST,XERROR
C***END PROLOGUE  PCHIM
C
C ----------------------------------------------------------------------
C
C  Change record:
C     82-02-01   1. Introduced  PCHST  to reduce possible over/under-
C                   flow problems.
C                2. Rearranged derivative formula for same reason.
C     82-06-02   1. Modified end conditions to be continuous functions
C                   of data when monotonicity switches in next interval.
C                2. Modified formulas so end conditions are less prone
C                   of over/underflow problems.
C     82-08-03   Minor cosmetic changes for release 1.
C
C ----------------------------------------------------------------------
C
C  Programming notes:
C
C     1. The function  PCHST(ARG1,ARG2)  is assumed to return zero if
C        either argument is zero, +1 if they are of the same sign, and
C        -1 if they are of opposite sign.
C     2. To produce a double precision version, simply:
C        a. Change PCHIM to DPCHIM wherever it occurs,
C        b. Change PCHST to DPCHST wherever it occurs,
C        c. Change all references to the Fortran intrinsics to their
C           double precision equivalents,
C        d. Change the real declarations to double precision, and
C        e. Change the constants ZERO and THREE to double precision.
C
C  DECLARE ARGUMENTS.
C
      INTEGER  N, INCFD, IERR
      REAL  X(N), F(INCFD,N), D(INCFD,N)
C
C  DECLARE LOCAL VARIABLES.
C
      INTEGER  I, NLESS1
      REAL  DEL1, DEL2, DMAX, DMIN, DRAT1, DRAT2, DSAVE,
     *      H1, H2, HSUM, HSUMT3, THREE, W1, W2, ZERO
      REAL  PCHST
      DATA  ZERO /0./,  THREE /3./
C
C  VALIDITY-CHECK ARGUMENTS.
C
C***FIRST EXECUTABLE STATEMENT  PCHIM
      IF ( N.LT.2 )  GO TO 5001
      IF ( INCFD.LT.1 )  GO TO 5002
      DO 1  I = 2, N
         IF ( X(I).LE.X(I-1) )  GO TO 5003
    1 CONTINUE
C
C  FUNCTION DEFINITION IS OK, GO ON.
C
      IERR = 0
      NLESS1 = N - 1
      H1 = X(2) - X(1)
      DEL1 = (F(1,2) - F(1,1))/H1
      DSAVE = DEL1
C
C  SPECIAL CASE N=2 -- USE LINEAR INTERPOLATION.
C
      IF (NLESS1 .GT. 1)  GO TO 10
      D(1,1) = DEL1
      D(1,N) = DEL1
      GO TO 5000
C
C  NORMAL CASE  (N .GE. 3).
C
   10 CONTINUE
      H2 = X(3) - X(2)
      DEL2 = (F(1,3) - F(1,2))/H2
C
C  SET D(1) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE
C     SHAPE-PRESERVING.
C
      HSUM = H1 + H2
      W1 = (H1 + HSUM)/HSUM
      W2 = -H1/HSUM
      D(1,1) = W1*DEL1 + W2*DEL2
      IF ( PCHST(D(1,1),DEL1) .LE. ZERO)  THEN
         D(1,1) = ZERO
      ELSE IF ( PCHST(DEL1,DEL2) .LT. ZERO)  THEN
C        NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES.
         DMAX = THREE*DEL1
         IF (ABS(D(1,1)) .GT. ABS(DMAX))  D(1,1) = DMAX
      ENDIF
C
C  LOOP THROUGH INTERIOR POINTS.
C
      DO 50  I = 2, NLESS1
         IF (I .EQ. 2)  GO TO 40
C
         H1 = H2
         H2 = X(I+1) - X(I)
         HSUM = H1 + H2
         DEL1 = DEL2
         DEL2 = (F(1,I+1) - F(1,I))/H2
   40    CONTINUE
C
C        SET D(I)=0 UNLESS DATA ARE STRICTLY MONOTONIC.
C
         D(1,I) = ZERO
         IF ( PCHST(DEL1,DEL2) )  42, 41, 45
C
C        COUNT NUMBER OF CHANGES IN DIRECTION OF MONOTONICITY.
C
   41    CONTINUE
         IF (DEL2 .EQ. ZERO)  GO TO 50
         IF ( PCHST(DSAVE,DEL2) .LT. ZERO)  IERR = IERR + 1
         DSAVE = DEL2
         GO TO 50
C
   42    CONTINUE
         IERR = IERR + 1
         DSAVE = DEL2
         GO TO 50
C
C        USE BRODLIE MODIFICATION OF BUTLAND FORMULA.
C
   45    CONTINUE
         HSUMT3 = HSUM+HSUM+HSUM
         W1 = (HSUM + H1)/HSUMT3
         W2 = (HSUM + H2)/HSUMT3
         DMAX = AMAX1( ABS(DEL1), ABS(DEL2) )
         DMIN = AMIN1( ABS(DEL1), ABS(DEL2) )
         DRAT1 = DEL1/DMAX
         DRAT2 = DEL2/DMAX
         D(1,I) = DMIN/(W1*DRAT1 + W2*DRAT2)
C
   50 CONTINUE
C
C  SET D(N) VIA NON-CENTERED THREE-POINT FORMULA, ADJUSTED TO BE
C     SHAPE-PRESERVING.
C
      W1 = -H2/HSUM
      W2 = (H2 + HSUM)/HSUM
      D(1,N) = W1*DEL1 + W2*DEL2
      IF ( PCHST(D(1,N),DEL2) .LE. ZERO)  THEN
         D(1,N) = ZERO
      ELSE IF ( PCHST(DEL1,DEL2) .LT. ZERO)  THEN
C        NEED DO THIS CHECK ONLY IF MONOTONICITY SWITCHES.
         DMAX = THREE*DEL2
         IF (ABS(D(1,N)) .GT. ABS(DMAX))  D(1,N) = DMAX
      ENDIF
C
C  NORMAL RETURN.
C
 5000 CONTINUE
      RETURN
C
C  ERROR RETURNS.
C
 5001 CONTINUE
C     N.LT.2 RETURN.
      IERR = -1
      CALL XERROR ('PCHIM -- NUMBER OF DATA POINTS LESS THAN TWO'
     *           , 44, IERR, 1)
      RETURN
C
 5002 CONTINUE
C     INCFD.LT.1 RETURN.
      IERR = -2
      CALL XERROR ('PCHIM -- INCREMENT LESS THAN ONE'
     *           , 32, IERR, 1)
      RETURN
C
 5003 CONTINUE
C     X-ARRAY NOT STRICTLY INCREASING.
      IERR = -3
      CALL XERROR ('PCHIM -- X-ARRAY NOT STRICTLY INCREASING'
     *           , 40, IERR, 1)
      RETURN
C------------- LAST LINE OF PCHIM FOLLOWS ------------------------------
      END
      SUBROUTINE PCHMC(N,X,F,D,INCFD,SKIP,ISMON,IERR)
C***BEGIN PROLOGUE  PCHMC
C***DATE WRITTEN   820518   (YYMMDD)
C***REVISION DATE  870707   (YYMMDD)
C***CATEGORY NO.  E1B
C***KEYWORDS  LIBRARY=SLATEC(PCHIP),
C             TYPE=SINGLE PRECISION(PCHMC-S DPCHMC-D),
C             CUBIC HERMITE INTERPOLATION,MONOTONE INTERPOLATION,
C             PIECEWISE CUBIC INTERPOLATION
C***AUTHOR  FRITSCH, F. N., (LLNL)
C             MATHEMATICS AND STATISTICS DIVISION
C             LAWRENCE LIVERMORE NATIONAL LABORATORY
C             P.O. BOX 808  (L-316)
C             LIVERMORE, CA  94550
C             FTS 532-4275, (415) 422-4275
C***PURPOSE  Check a cubic Hermite function for monotonicity.
C***DESCRIPTION
C
C          PCHMC:  Piecewise Cubic Hermite Monotonicity Checker.
C
C     Checks the cubic Hermite function defined by  N, X, F, D  for
C     monotonicity.
C
C     To provide compatibility with PCHIM and PCHIC, includes an
C     increment between successive values of the F- and D-arrays.
C
C ----------------------------------------------------------------------
C
C  Calling sequence:
C
C        PARAMETER  (INCFD = ...)
C        INTEGER  N, ISMON(N), IERR
C        REAL  X(N), F(INCFD,N), D(INCFD,N)
C        LOGICAL  SKIP
C
C        CALL  PCHMC (N, X, F, D, INCFD, SKIP, ISMON, IERR)
C
C   Parameters:
C
C     N -- (input) number of data points.  (Error return if N.LT.2 .)
C
C     X -- (input) real array of independent variable values.  The
C           elements of X must be strictly increasing:
C                X(I-1) .LT. X(I),  I = 2(1)N.
C           (Error return if not.)
C
C     F -- (input) real array of function values.  F(1+(I-1)*INCFD) is
C           the value corresponding to X(I).
C
C     D -- (input) real array of derivative values.  D(1+(I-1)*INCFD) is
C           the value corresponding to X(I).
C
C     INCFD -- (input) increment between successive values in F and D.
C           (Error return if  INCFD.LT.1 .)
C
C     SKIP -- (input/output) logical variable which should be set to
C           .TRUE. if the user wishes to skip checks for validity of
C           preceding parameters, or to .FALSE. otherwise.
C           This will save time in case these checks have already
C           been performed.
C           SKIP will be set to .TRUE. on normal return.
C
C     ISMON -- (output) integer array indicating on which intervals the
C           PCH function defined by  N, X, F, D  is monotonic.
C           For data interval [X(I),X(I+1)],
C             ISMON(I) = -1  if function is strictly decreasing;
C             ISMON(I) =  0  if function is constant;
C             ISMON(I) =  1  if function is strictly increasing;
C             ISMON(I) =  2  if function is non-monotonic;
C             ISMON(I) =  3  if unable to determine.  (This means that
C                            the D-values are near the boundary of the
C                            monotonicity region.  A small increase pro-
C                            duces non-monotonicity; decrease, strict
C                            monotonicity.)
C           The above applies to I=1(1)N-1.  ISMON(N) indicates whether
C              the entire function is monotonic on [X(1),X(N)].
C
C     IERR -- (output) error flag.
C           Normal return:
C              IERR = 0  (no errors).
C           "Recoverable" errors:
C              IERR = -1  if N.LT.2 .
C              IERR = -2  if INCFD.LT.1 .
C              IERR = -3  if the X-array is not strictly increasing.
C          (The ISMON-array has not been changed in any of these cases.)
C               NOTE:  The above errors are checked in the order listed,
C                   and following arguments have **NOT** been validated.
C
C***REFERENCES  F.N.FRITSCH AND R.E.CARLSON, 'MONOTONE PIECEWISE CUBIC
C                 INTERPOLATION,' SIAM J.NUMER.ANAL. 17, 2 (APRIL 1980),
C                 238-246.
C***ROUTINES CALLED  CHFMC,XERROR
C***END PROLOGUE  PCHMC
C
C ----------------------------------------------------------------------
C
C  Change record:
C     82-08-04   Converted to SLATEC library version.
C     83-12-01   Reversed order of subscripts of F and D, so that the
C                routine will work properly when INCFD.GT.1  .
C     87-07-07   Minor cosmetic changes to prologue.
C
C ----------------------------------------------------------------------
C
C  Programming notes:
C
C     To produce a double precision version, simply:
C        a. Change PCHMC to DPCHMC wherever it occurs,
C        b. Change CHFMC to DCHFMC wherever it occurs, and
C        c. Change the real declarations to double precision.
C
C  DECLARE ARGUMENTS.
C
      INTEGER  N, INCFD, ISMON(N), IERR
      REAL  X(N), F(INCFD,N), D(INCFD,N)
      LOGICAL  SKIP
C
C  DECLARE LOCAL VARIABLES.
C
      INTEGER  I, NSEG
      REAL  DELTA
      INTEGER  CHFMC
C
C  VALIDITY-CHECK ARGUMENTS.
C
C***FIRST EXECUTABLE STATEMENT  PCHMC
      IF (SKIP)  GO TO 5
C
      IF ( N.LT.2 )  GO TO 5001
      IF ( INCFD.LT.1 )  GO TO 5002
      DO 1  I = 2, N
         IF ( X(I).LE.X(I-1) )  GO TO 5003
    1 CONTINUE
      SKIP = .TRUE.
C
C  FUNCTION DEFINITION IS OK -- GO ON.
C
    5 CONTINUE
      NSEG = N - 1
      DO 90  I = 1, NSEG
         DELTA = (F(1,I+1)-F(1,I))/(X(I+1)-X(I))
C                   -------------------------------
         ISMON(I) = CHFMC (D(1,I), D(1,I+1), DELTA)
C                   -------------------------------
         IF (I .EQ. 1)  THEN
            ISMON(N) = ISMON(1)
         ELSE
C           NEED TO FIGURE OUT CUMULATIVE MONOTONICITY FROM FOLLOWING
C           'MULTIPLICATION TABLE'--
C
C                    *      I S M O N (I)
C                     *  -1   0   1   2   3
C               I      *--------------------*
C               S   -1 I -1  -1   2   2   3 I
C               M    0 I -1   0   1   2   3 I
C               O    1 I  2   1   1   2   3 I
C               N    2 I  2   2   2   2   2 I
C              (N)   3 I  3   3   3   2   3 I
C                      *--------------------*
C
C           IF EQUAL OR ALREADY DECLARED NONMONOTONIC, NO CHANGE NEEDED.
            IF ((ISMON(I).NE.ISMON(N)) .AND. (ISMON(N).NE.2))  THEN
               IF ( MAX0(ISMON(I), ISMON(N)) .GT. 1)  THEN
C                 AT LEAST ONE IS EITHER 'NO' OR 'MAYBE'.
                  IF (ISMON(I) .EQ. 2)  THEN
                     ISMON(N) = 2
                  ELSE
                     ISMON(N) = 3
                  ENDIF
               ELSE IF (ISMON(I)*ISMON(N) .LT. 0)  THEN
C                 BOTH MONOTONIC, BUT IN OPPOSITE SENSES.
                  ISMON(N) = 2
               ELSE
C                 AT THIS POINT, ONE IS ZERO, THE OTHER IS +-1.
                  ISMON(N) = ISMON(N) + ISMON(I)
               ENDIF
            ENDIF
         ENDIF
   90 CONTINUE
C
C  NORMAL RETURN.
C
      IERR = 0
      RETURN
C
C  ERROR RETURNS.
C
 5001 CONTINUE
C     N.LT.2 RETURN.
      IERR = -1
      CALL XERROR ('PCHMC -- NUMBER OF DATA POINTS LESS THAN TWO'
     *           , 44, IERR, 1)
      RETURN
C
 5002 CONTINUE
C     INCFD.LT.1 RETURN.
      IERR = -2
      CALL XERROR ('PCHMC -- INCREMENT LESS THAN ONE'
     *           , 32, IERR, 1)
      RETURN
C
 5003 CONTINUE
C     X-ARRAY NOT STRICTLY INCREASING.
      IERR = -3
      CALL XERROR ('PCHMC -- X-ARRAY NOT STRICTLY INCREASING'
     *           , 40, IERR, 1)
      RETURN
C------------- LAST LINE OF PCHMC FOLLOWS ------------------------------
      END
      REAL FUNCTION PCHQA(N,X,F,D,A,B,IERR)
C***BEGIN PROLOGUE  PCHQA
C***DATE WRITTEN   870829   (YYMMDD)
C***REVISION DATE  870829   (YYMMDD)
C***CATEGORY NO.  E3,H2A2
C***KEYWORDS  EASY TO USE CUBIC HERMITE OR SPLINE INTEGRATION
C             NUMERICAL INTEGRATION, QUADRATURE
C***AUTHOR  KAHANER, D.K., (NBS)
C             SCIENTIFIC COMPUTING DIVISION
C             NATIONAL BUREAU OF STANDARDS
C             ROOM A161, TECHNOLOGY BUILDING
C             GAITHERSBURG, MARYLAND 20899
C             (301) 975-3808
C***PURPOSE  Evaluates the definite integral of a piecewise cubic Hermite
C            or spline function over an arbitrary interval, easy to use.
C***DESCRIPTION
C
C          PCHQA:  Piecewise Cubic Hermite or Spline Integrator,
C                  Arbitrary limits, Easy to Use.
C
C          From the book "Numerical Methods and Software"
C                  by  D. Kahaner, C. Moler, S. Nash
C                          Prentice Hall 1988
C
C     Evaluates the definite integral of the cubic Hermite or spline
C     function defined by  N, X, F, D  over the interval [A, B].  This
C     is an easy to use driver for the routine PCHIA by F.N. Fritsch
C     described in reference (2) below. That routine also has other
C     capabilities.
C ----------------------------------------------------------------------
C
C  Calling sequence:
C
C           VALUE = PCHQA (N, X, F, D, A, B, IERR)
C
C     INTEGER  N, IERR
C     REAL  X(N), F(N), D(N), A, B
C
C   Parameters:
C
C     VALUE -- (output) VALUE of the requested integral.
C
C     N -- (input) number of data points.  (Error return if N.LT.2 .)
C
C     X -- (input) real array of independent variable values.  The
C           elements of X must be strictly increasing:
C                X(I-1) .LT. X(I),  I = 2(1)N.
C           (Error return if not.)
C
C     F -- (input) real array of function values.  F(I) is
C           the value corresponding to X(I).
C
C     D -- (input) real array of derivative values.  D(I) is
C           the value corresponding to X(I).
C
C     A,B -- (input) the limits of integration.
C           NOTE:  There is no requirement that [A,B] be contained in
C                  [X(1),X(N)].  However, the resulting integral value
C                  will be highly suspect, if not.
C
C     IERR -- (output) error flag.
C           Normal return:
C              IERR = 0  (no errors).
C           Warning errors:
C              IERR = 1  if  A  is outside the interval [X(1),X(N)].
C              IERR = 2  if  B  is outside the interval [X(1),X(N)].
C              IERR = 3  if both of the above are true.  (Note that this
C                        means that either [A,B] contains data interval
C                        or the intervals do not intersect at all.)
C           "Recoverable" errors:
C              IERR = -1  if N.LT.2 .
C              IERR = -3  if the X-array is not strictly increasing.
C                (Value has not been computed in any of these cases.)
C               NOTE:  The above errors are checked in the order listed,
C                   and following arguments have **NOT** been validated.
C
C***REFERENCES  1. F.N.FRITSCH AND R.E.CARLSON, 'MONOTONE PIECEWISE
C                 CUBIC INTERPOLATION,' SIAM J.NUMER.ANAL. 17, 2 (APRIL
C                 1980), 238-246.
C               2. F.N.FRITSCH, 'PIECEWISE CUBIC HERMITE INTERPOLATION
C                 PACKAGE, FINAL SPECIFICATIONS', LAWRENCE LIVERMORE
C                 NATIONAL LABORATORY, COMPUTER DOCUMENTATION UCID-30194,
C                 AUGUST 1982.
C***ROUTINES CALLED  PCHIA
C***END PROLOGUE  PCHQA
      INTEGER  N, IERR
      REAL  X(N), F(N), D(N), A, B
C
C  DECLARE LOCAL VARIABLES.
C
      INTEGER  INCFD
      REAL  PCHIA
      LOGICAL SKIP
C
C  INITIALIZE.
C
      DATA  INCFD /1/
      DATA  SKIP /.TRUE./
C
C
C***FIRST EXECUTABLE STATEMENT  PCHQA

      PCHQA  =  PCHIA( N, X, F, D, INCFD, SKIP, A, B, IERR )
C
C ERROR MESSAGES ARE FROM LOWER LEVEL ROUTINES
      RETURN
C
C------------- LAST LINE OF PCHQA FOLLOWS ------------------------------
      END
      SUBROUTINE PCHSP(IC,VC,N,X,F,D,INCFD,WK,NWK,IERR)
C***BEGIN PROLOGUE  PCHSP
C***DATE WRITTEN   820503   (YYMMDD)
C***REVISION DATE  870707   (YYMMDD)
C***CATEGORY NO.  E1B
C***KEYWORDS  LIBRARY=SLATEC(PCHIP),
C             TYPE=SINGLE PRECISION(PCHSP-S DPCHSP-D),
C             CUBIC HERMITE INTERPOLATION,PIECEWISE CUBIC INTERPOLATION,
C             SPLINE INTERPOLATION
C***AUTHOR  FRITSCH, F. N., (LLNL)
C             MATHEMATICS AND STATISTICS DIVISION
C             LAWRENCE LIVERMORE NATIONAL LABORATORY
C             P.O. BOX 808  (L-316)
C             LIVERMORE, CA  94550
C             FTS 532-4275, (415) 422-4275
C***PURPOSE  Set derivatives needed to determine the Hermite represen-
C            tation of the cubic spline interpolant to given data, with
C            specified boundary conditions.
C***DESCRIPTION
C
C          PCHSP:   Piecewise Cubic Hermite Spline
C
C     Computes the Hermite representation of the cubic spline inter-
C     polant to the data given in X and F satisfying the boundary
C     conditions specified by IC and VC.
C
C     To facilitate two-dimensional applications, includes an increment
C     between successive values of the F- and D-arrays.
C
C     The resulting piecewise cubic Hermite function may be evaluated
C     by PCHFE or PCHFD.
C
C     NOTE:  This is a modified version of C. de Boor'S cubic spline
C            routine CUBSPL.
C
C ----------------------------------------------------------------------
C
C  Calling sequence:
C
C        PARAMETER  (INCFD = ...)
C        INTEGER  IC(2), N, NWK, IERR
C        REAL  VC(2), X(N), F(INCFD,N), D(INCFD,N), WK(NWK)
C
C        CALL  PCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR)
C
C   Parameters:
C
C     IC -- (input) integer array of length 2 specifying desired
C           boundary conditions:
C           IC(1) = IBEG, desired condition at beginning of data.
C           IC(2) = IEND, desired condition at end of data.
C
C           IBEG = 0  to set D(1) so that the third derivative is con-
C              tinuous at X(2).  This is the "not a knot" condition
C              provided by de Boor'S cubic spline routine CUBSPL.
C              < This is the default boundary condition. >
C           IBEG = 1  if first derivative at X(1) is given in VC(1).
C           IBEG = 2  if second derivative at X(1) is given in VC(1).
C           IBEG = 3  to use the 3-point difference formula for D(1).
C                     (Reverts to the default b.c. if N.LT.3 .)
C           IBEG = 4  to use the 4-point difference formula for D(1).
C                     (Reverts to the default b.c. if N.LT.4 .)
C          NOTES:
C           1. An error return is taken if IBEG is out of range.
C           2. For the "natural" boundary condition, use IBEG=2 and
C              VC(1)=0.
C
C           IEND may take on the same values as IBEG, but applied to
C           derivative at X(N).  In case IEND = 1 or 2, the value is
C           given in VC(2).
C
C          NOTES:
C           1. An error return is taken if IEND is out of range.
C           2. For the "natural" boundary condition, use IEND=2 and
C              VC(2)=0.
C
C     VC -- (input) real array of length 2 specifying desired boundary
C           values, as indicated above.
C           VC(1) need be set only if IC(1) = 1 or 2 .
C           VC(2) need be set only if IC(2) = 1 or 2 .
C
C     N -- (input) number of data points.  (Error return if N.LT.2 .)
C
C     X -- (input) real array of independent variable values.  The
C           elements of X must be strictly increasing:
C                X(I-1) .LT. X(I),  I = 2(1)N.
C           (Error return if not.)
C
C     F -- (input) real array of dependent variable values to be inter-
C           polated.  F(1+(I-1)*INCFD) is value corresponding to X(I).
C
C     D -- (output) real array of derivative values at the data points.
C           These values will determine the cubic spline interpolant
C           with the requested boundary conditions.
C           The value corresponding to X(I) is stored in
C                D(1+(I-1)*INCFD),  I=1(1)N.
C           No other entries in D are changed.
C
C     INCFD -- (input) increment between successive values in F and D.
C           This argument is provided primarily for 2-D applications.
C           (Error return if  INCFD.LT.1 .)
C
C     WK -- (scratch) real array of working storage.
C
C     NWK -- (input) length of work array.
C           (Error return if NWK.LT.2*N .)
C
C     IERR -- (output) error flag.
C           Normal return:
C              IERR = 0  (no errors).
C           "Recoverable" errors:
C              IERR = -1  if N.LT.2 .
C              IERR = -2  if INCFD.LT.1 .
C              IERR = -3  if the X-array is not strictly increasing.
C              IERR = -4  if IBEG.LT.0 or IBEG.GT.4 .
C              IERR = -5  if IEND.LT.0 of IEND.GT.4 .
C              IERR = -6  if both of the above are true.
C              IERR = -7  if NWK is too small.
C               NOTE:  The above errors are checked in the order listed,
C                   and following arguments have **NOT** been validated.
C             (The D-array has not been changed in any of these cases.)
C              IERR = -8  in case of trouble solving the linear system
C                         for the interior derivative values.
C             (The D-array may have been changed in this case.)
C             (             Do **NOT** use it!                )
C
C***REFERENCES  CARL DE BOOR, A PRACTICAL GUIDE TO SPLINES, SPRINGER-
C                 VERLAG (NEW YORK, 1978), PP. 53-59.
C***ROUTINES CALLED  PCHDF,XERROR
C***END PROLOGUE  PCHSP
C
C ----------------------------------------------------------------------
C
C  Change record:
C     82-08-04   Converted to SLATEC library version.
C     87-07-07   Minor cosmetic changes to prologue.
C
C ----------------------------------------------------------------------
C
C  Programming notes:
C
C     To produce a double precision version, simply:
C        a. Change PCHSP to DPCHSP wherever it occurs,
C        b. Change the real declarations to double precision, and
C        c. Change the constants ZERO, HALF, ... to double precision.
C
C  DECLARE ARGUMENTS.
C
      INTEGER  IC(2), N, INCFD, NWK, IERR
      REAL  VC(2), X(N), F(INCFD,N), D(INCFD,N), WK(2,N)
C
C  DECLARE LOCAL VARIABLES.
C
      INTEGER  IBEG, IEND, INDEX, J, NM1
      REAL  G, HALF, ONE, STEMP(3), THREE, TWO, XTEMP(4), ZERO
      REAL  PCHDF
C
      DATA  ZERO /0./,  HALF /0.5/,  ONE /1./,  TWO /2./,  THREE /3./
C
C  VALIDITY-CHECK ARGUMENTS.
C
C***FIRST EXECUTABLE STATEMENT  PCHSP
      IF ( N.LT.2 )  GO TO 5001
      IF ( INCFD.LT.1 )  GO TO 5002
      DO 1  J = 2, N
         IF ( X(J).LE.X(J-1) )  GO TO 5003
    1 CONTINUE
C
      IBEG = IC(1)
      IEND = IC(2)
      IERR = 0
      IF ( (IBEG.LT.0).OR.(IBEG.GT.4) )  IERR = IERR - 1
      IF ( (IEND.LT.0).OR.(IEND.GT.4) )  IERR = IERR - 2
      IF ( IERR.LT.0 )  GO TO 5004
C
C  FUNCTION DEFINITION IS OK -- GO ON.
C
      IF ( NWK .LT. 2*N )  GO TO 5007
C
C  COMPUTE FIRST DIFFERENCES OF X SEQUENCE AND STORE IN WK(1,.). ALSO,
C  COMPUTE FIRST DIVIDED DIFFERENCE OF DATA AND STORE IN WK(2,.).
      DO 5  J=2,N
         WK(1,J) = X(J) - X(J-1)
         WK(2,J) = (F(1,J) - F(1,J-1))/WK(1,J)
    5 CONTINUE
C
C  SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL.
C
      IF ( IBEG.GT.N )  IBEG = 0
      IF ( IEND.GT.N )  IEND = 0
C
C  SET UP FOR BOUNDARY CONDITIONS.
C
      IF ( (IBEG.EQ.1).OR.(IBEG.EQ.2) )  THEN
         D(1,1) = VC(1)
      ELSE IF (IBEG .GT. 2)  THEN
C        PICK UP FIRST IBEG POINTS, IN REVERSE ORDER.
         DO 10  J = 1, IBEG
            INDEX = IBEG-J+1
C           INDEX RUNS FROM IBEG DOWN TO 1.
            XTEMP(J) = X(INDEX)
            IF (J .LT. IBEG)  STEMP(J) = WK(2,INDEX)
   10    CONTINUE
C                 --------------------------------
         D(1,1) = PCHDF (IBEG, XTEMP, STEMP, IERR)
C                 --------------------------------
         IF (IERR .NE. 0)  GO TO 5009
         IBEG = 1
      ENDIF
C
      IF ( (IEND.EQ.1).OR.(IEND.EQ.2) )  THEN
         D(1,N) = VC(2)
      ELSE IF (IEND .GT. 2)  THEN
C        PICK UP LAST IEND POINTS.
         DO 15  J = 1, IEND
            INDEX = N-IEND+J
C           INDEX RUNS FROM N+1-IEND UP TO N.
            XTEMP(J) = X(INDEX)
            IF (J .LT. IEND)  STEMP(J) = WK(2,INDEX+1)
   15    CONTINUE
C                 --------------------------------
         D(1,N) = PCHDF (IEND, XTEMP, STEMP, IERR)
C                 --------------------------------
         IF (IERR .NE. 0)  GO TO 5009
         IEND = 1
      ENDIF
C
C --------------------( BEGIN CODING FROM CUBSPL )--------------------
C
C  **** A TRIDIAGONAL LINEAR SYSTEM FOR THE UNKNOWN SLOPES S(J) OF
C  F  AT X(J), J=1,...,N, IS GENERATED AND THEN SOLVED BY GAUSS ELIM-
C  INATION, WITH S(J) ENDING UP IN D(1,J), ALL J.
C     WK(1,.) AND WK(2,.) ARE USED FOR TEMPORARY STORAGE.
C
C  CONSTRUCT FIRST EQUATION FROM FIRST BOUNDARY CONDITION, OF THE FORM
C             WK(2,1)*S(1) + WK(1,1)*S(2) = D(1,1)
C
      IF (IBEG .EQ. 0)  THEN
         IF (N .EQ. 2)  THEN
C           NO CONDITION AT LEFT END AND N = 2.
            WK(2,1) = ONE
            WK(1,1) = ONE
            D(1,1) = TWO*WK(2,2)
         ELSE
C           NOT-A-KNOT CONDITION AT LEFT END AND N .GT. 2.
            WK(2,1) = WK(1,3)
            WK(1,1) = WK(1,2) + WK(1,3)
            D(1,1) =((WK(1,2) + TWO*WK(1,1))*WK(2,2)*WK(1,3)
     *                        + WK(1,2)**2*WK(2,3)) / WK(1,1)
         ENDIF
      ELSE IF (IBEG .EQ. 1)  THEN
C        SLOPE PRESCRIBED AT LEFT END.
         WK(2,1) = ONE
         WK(1,1) = ZERO
      ELSE
C        SECOND DERIVATIVE PRESCRIBED AT LEFT END.
         WK(2,1) = TWO
         WK(1,1) = ONE
         D(1,1) = THREE*WK(2,2) - HALF*WK(1,2)*D(1,1)
      ENDIF
C
C  IF THERE ARE INTERIOR KNOTS, GENERATE THE CORRESPONDING EQUATIONS AND
C  CARRY OUT THE FORWARD PASS OF GAUSS ELIMINATION, AFTER WHICH THE J-TH
C  EQUATION READS    WK(2,J)*S(J) + WK(1,J)*S(J+1) = D(1,J).
C
      NM1 = N-1
      IF (NM1 .GT. 1)  THEN
         DO 20 J=2,NM1
            IF (WK(2,J-1) .EQ. ZERO)  GO TO 5008
            G = -WK(1,J+1)/WK(2,J-1)
            D(1,J) = G*D(1,J-1)
     *                  + THREE*(WK(1,J)*WK(2,J+1) + WK(1,J+1)*WK(2,J))
            WK(2,J) = G*WK(1,J-1) + TWO*(WK(1,J) + WK(1,J+1))
   20    CONTINUE
      ENDIF
C
C  CONSTRUCT LAST EQUATION FROM SECOND BOUNDARY CONDITION, OF THE FORM
C           (-G*WK(2,N-1))*S(N-1) + WK(2,N)*S(N) = D(1,N)
C
C     IF SLOPE IS PRESCRIBED AT RIGHT END, ONE CAN GO DIRECTLY TO BACK-
C     SUBSTITUTION, SINCE ARRAYS HAPPEN TO BE SET UP JUST RIGHT FOR IT
C     AT THIS POINT.
      IF (IEND .EQ. 1)  GO TO 30
C
      IF (IEND .EQ. 0)  THEN
         IF (N.EQ.2 .AND. IBEG.EQ.0)  THEN
C           NOT-A-KNOT AT RIGHT ENDPOINT AND AT LEFT ENDPOINT AND N = 2.
            D(1,2) = WK(2,2)
            GO TO 30
         ELSE IF ((N.EQ.2) .OR. (N.EQ.3 .AND. IBEG.EQ.0))  THEN
C           EITHER (N=3 AND NOT-A-KNOT ALSO AT LEFT) OR (N=2 AND *NOT*
C           NOT-A-KNOT AT LEFT END POINT).
            D(1,N) = TWO*WK(2,N)
            WK(2,N) = ONE
            IF (WK(2,N-1) .EQ. ZERO)  GO TO 5008
            G = -ONE/WK(2,N-1)
         ELSE
C           NOT-A-KNOT AND N .GE. 3, AND EITHER N.GT.3 OR  ALSO NOT-A-
C           KNOT AT LEFT END POINT.
            G = WK(1,N-1) + WK(1,N)
C           DO NOT NEED TO CHECK FOLLOWING DENOMINATORS (X-DIFFERENCES).
            D(1,N) = ((WK(1,N)+TWO*G)*WK(2,N)*WK(1,N-1)
     *                  + WK(1,N)**2*(F(1,N-1)-F(1,N-2))/WK(1,N-1))/G
            IF (WK(2,N-1) .EQ. ZERO)  GO TO 5008
            G = -G/WK(2,N-1)
            WK(2,N) = WK(1,N-1)
         ENDIF
      ELSE
C        SECOND DERIVATIVE PRESCRIBED AT RIGHT ENDPOINT.
         D(1,N) = THREE*WK(2,N) + HALF*WK(1,N)*D(1,N)
         WK(2,N) = TWO
         IF (WK(2,N-1) .EQ. ZERO)  GO TO 5008
         G = -ONE/WK(2,N-1)
      ENDIF
C
C  COMPLETE FORWARD PASS OF GAUSS ELIMINATION.
C
      WK(2,N) = G*WK(1,N-1) + WK(2,N)
      IF (WK(2,N) .EQ. ZERO)   GO TO 5008
      D(1,N) = (G*D(1,N-1) + D(1,N))/WK(2,N)
C
C  CARRY OUT BACK SUBSTITUTION
C
   30 CONTINUE
      DO 40 J=NM1,1,-1
         IF (WK(2,J) .EQ. ZERO)  GO TO 5008
         D(1,J) = (D(1,J) - WK(1,J)*D(1,J+1))/WK(2,J)
   40 CONTINUE
C --------------------(  END  CODING FROM CUBSPL )--------------------
C
C  NORMAL RETURN.
C
      RETURN
C
C  ERROR RETURNS.
C
 5001 CONTINUE
C     N.LT.2 RETURN.
      IERR = -1
      CALL XERROR ('PCHSP -- NUMBER OF DATA POINTS LESS THAN TWO'
     *           , 44, IERR, 1)
      RETURN
C
 5002 CONTINUE
C     INCFD.LT.1 RETURN.
      IERR = -2
      CALL XERROR ('PCHSP -- INCREMENT LESS THAN ONE'
     *           , 32, IERR, 1)
      RETURN
C
 5003 CONTINUE
C     X-ARRAY NOT STRICTLY INCREASING.
      IERR = -3
      CALL XERROR ('PCHSP -- X-ARRAY NOT STRICTLY INCREASING'
     *           , 40, IERR, 1)
      RETURN
C
 5004 CONTINUE
C     IC OUT OF RANGE RETURN.
      IERR = IERR - 3
      CALL XERROR ('PCHSP -- IC OUT OF RANGE'
     *           , 24, IERR, 1)
      RETURN
C
 5007 CONTINUE
C     NWK TOO SMALL RETURN.
      IERR = -7
      CALL XERROR ('PCHSP -- WORK ARRAY TOO SMALL'
     *           , 29, IERR, 1)
      RETURN
C
 5008 CONTINUE
C     SINGULAR SYSTEM.
C   *** THEORETICALLY, THIS CAN ONLY OCCUR IF SUCCESSIVE X-VALUES   ***
C   *** ARE EQUAL, WHICH SHOULD ALREADY HAVE BEEN CAUGHT (IERR=-3). ***
      IERR = -8
      CALL XERROR ('PCHSP -- SINGULAR LINEAR SYSTEM'
     *           , 31, IERR, 1)
      RETURN
C
 5009 CONTINUE
C     ERROR RETURN FROM PCHDF.
C   *** THIS CASE SHOULD NEVER OCCUR ***
      IERR = -9
      CALL XERROR ('PCHSP -- ERROR RETURN FROM PCHDF'
     *           , 32, IERR, 1)
      RETURN
C------------- LAST LINE OF PCHSP FOLLOWS ------------------------------
      END
      REAL FUNCTION PCHST(ARG1,ARG2)
C***BEGIN PROLOGUE  PCHST
C***REFER TO  PCHCE,PCHCI,PCHCS,PCHIM
C***ROUTINES CALLED  (NONE)
C***DESCRIPTION
C
C         PCHST:  PCHIP Sign-Testing Routine.
C
C
C     Returns:
C        -1. if ARG1 and ARG2 are of opposite sign.
C         0. if either argument is zero.
C        +1. if ARG1 and ARG2 are of the same sign.
C
C     The object is to do this without multiplying ARG1*ARG2, to avoid
C     possible over/underflow problems.
C
C  Fortran intrinsics used:  SIGN.
C
C ----------------------------------------------------------------------
C
C  Programmed by:  Fred N. Fritsch,  FTS 532-4275, (415) 422-4275,
C                  Mathematics and Statistics Division,
C                  Lawrence Livermore National Laboratory.
C
C  Change record:
C     82-08-05   Converted to SLATEC library version.
C
C ----------------------------------------------------------------------
C
C  Programming notes:
C
C     To produce a double precision version, simply:
C        a. Change PCHST to DPCHST wherever it occurs,
C        b. Change all references to the Fortran intrinsics to their
C           double presision equivalents,
C        c. Change the real declarations to double precision, and
C        d. Change the constants  ZERO  and  ONE  to double precision.
C***END PROLOGUE  PCHST
      REAL  ARG1, ARG2
C
C  DECLARE LOCAL VARIABLES.
C
      REAL  ONE, ZERO
      DATA  ZERO /0./,  ONE /1./
C
C  PERFORM THE TEST.
C
C***FIRST EXECUTABLE STATEMENT  PCHST
      PCHST = SIGN(ONE,ARG1) * SIGN(ONE,ARG2)
      IF ((ARG1.EQ.ZERO) .OR. (ARG2.EQ.ZERO))  PCHST = ZERO
C
      RETURN
C------------- LAST LINE OF PCHST FOLLOWS ------------------------------
      END
      SUBROUTINE PCHSW(DFMAX,IEXTRM,D1,D2,H,SLOPE,IERR)
C***BEGIN PROLOGUE  PCHSW
C***REFER TO  PCHCS
C***ROUTINES CALLED  R1MACH,XERROR
C***REVISION DATE  870707   (YYMMDD)
C***DESCRIPTION
C
C         PCHSW:  PCHCS Switch Excursion Limiter.
C
C     Called by  PCHCS  to adjust D1 and D2 if necessary to insure that
C     the extremum on this interval is not further than DFMAX from the
C     extreme data value.
C
C ----------------------------------------------------------------------
C
C  Calling sequence:
C
C        INTEGER  IEXTRM, IERR
C        REAL  DFMAX, D1, D2, H, SLOPE
C
C        CALL  PCHSW (DFMAX, IEXTRM, D1, D2, H, SLOPE, IERR)
C
C   Parameters:
C
C     DFMAX -- (input) maximum allowed difference between F(IEXTRM) and
C           the cubic determined by derivative values D1,D2.  (assumes
C           DFMAX.GT.0.)
C
C     IEXTRM -- (input) index of the extreme data value.  (assumes
C           IEXTRM = 1 or 2 .  Any value .NE.1 is treated as 2.)
C
C     D1,D2 -- (input) derivative values at the ends of the interval.
C           (Assumes D1*D2 .LE. 0.)
C          (output) may be modified if necessary to meet the restriction
C           imposed by DFMAX.
C
C     H -- (input) interval length.  (Assumes  H.GT.0.)
C
C     SLOPE -- (input) data SLOPE on the interval.
C
C     IERR -- (output) error flag.  should be zero.
C           If IERR=-1, assumption on D1 and D2 is not satisfied.
C           If IERR=-2, quadratic equation locating extremum has
C                       negative descriminant (should never occur).
C
C    -------
C    WARNING:  This routine does no validity-checking of arguments.
C    -------
C
C
C  Fortran intrinsics used:  ABS, SIGN, SQRT.
C
C***END PROLOGUE  PCHSW
C
C ----------------------------------------------------------------------
C
C  Programmed by:  Fred N. Fritsch,  FTS 532-4275, (415) 422-4275,
C                  Mathematics and Statistics Division,
C                  Lawrence Livermore National Laboratory.
C
C  Change record:
C     82-08-05   Converted to SLATEC library version.
C     87-07-07   Replaced DATA statement for SMALL with a use of R1MACH.
C
C ----------------------------------------------------------------------
C
C  Programming notes:
C
C     To produce a double precision version, simply:
C        a. Change PCHSW to DPCHSW wherever it occurs,
C        b. Change PCHCS to DPCHCS wherever it occurs,
C        c. Change R1MACH to D1MACH wherever it occurs,
C        d. Change all references to the Fortran intrinsics to their
C           double precision equivalents,
C        e. Change the real declarations to double precision, and
C        f. Change constants ZERO, ONE, TWO, ... to double precision.
C
C  DECLARE ARGUMENTS.
C
      INTEGER  IEXTRM, IERR
      REAL  DFMAX, D1, D2, H, SLOPE
C
C  DECLARE LOCAL VARIABLES.
C
      REAL  CP, DMAX, FACT, LAMBDA, NU, ONE, PHI, RADCAL, RHO, SIGMA,
     *      SMALL, THAT, THIRD, THREE, TWO, ZERO
      DATA  ZERO /0./,  ONE /1./,  TWO /2./,  THREE /3./, FACT /100./
C        THIRD SHOULD BE SLIGHTLY LESS THAN 1/3.
      DATA  THIRD /0.33333/
C
C  NOTATION AND GENERAL REMARKS.
C
C     RHO IS THE RATIO OF THE DATA SLOPE TO THE DERIVATIVE BEING TESTED.
C     LAMBDA IS THE RATIO OF D2 TO D1.
C     THAT = T-HAT(RHO) IS THE NORMALIZED LOCATION OF THE EXTREMUM.
C     PHI IS THE NORMALIZED VALUE OF P(X)-F1 AT X = XHAT = X-HAT(RHO),
C           WHERE  THAT = (XHAT - X1)/H .
C        THAT IS, P(XHAT)-F1 = D*H*PHI,  WHERE D=D1 OR D2.
C     SIMILARLY,  P(XHAT)-F2 = D*H*(PHI-RHO) .
C
C      SMALL SHOULD BE A FEW ORDERS OF MAGNITUDE GREATER THAN MACHEPS.
C***FIRST EXECUTABLE STATEMENT  PCHSW
      SMALL = FACT*R1MACH(4)
C
C  DO MAIN CALCULATION.
C
      IF (D1 .EQ. ZERO)  THEN
C
C        SPECIAL CASE -- D1.EQ.ZERO .
C
C          IF D2 IS ALSO ZERO, THIS ROUTINE SHOULD NOT HAVE BEEN CALLED.
         IF (D2 .EQ. ZERO)  GO TO 5001
C
         RHO = SLOPE/D2
C          EXTREMUM IS OUTSIDE INTERVAL WHEN RHO .GE. 1/3 .
         IF (RHO .GE. THIRD)  GO TO 5000
         THAT = (TWO*(THREE*RHO-ONE)) / (THREE*(TWO*RHO-ONE))
         PHI = THAT**2 * ((THREE*RHO-ONE)/THREE)
C
C          CONVERT TO DISTANCE FROM F2 IF IEXTRM.NE.1 .
         IF (IEXTRM .NE. 1)  PHI = PHI - RHO
C
C          TEST FOR EXCEEDING LIMIT, AND ADJUST ACCORDINGLY.
         DMAX = DFMAX / (H*ABS(PHI))
         IF (ABS(D2) .GT. DMAX)  D2 = SIGN (DMAX, D2)
      ELSE
C
         RHO = SLOPE/D1
         LAMBDA = -D2/D1
         IF (D2 .EQ. ZERO)  THEN
C
C           SPECIAL CASE -- D2.EQ.ZERO .
C
C             EXTREMUM IS OUTSIDE INTERVAL WHEN RHO .GE. 1/3 .
            IF (RHO .GE. THIRD)  GO TO 5000
            CP = TWO - THREE*RHO
            NU = ONE - TWO*RHO
            THAT = ONE / (THREE*NU)
         ELSE
            IF (LAMBDA .LE. ZERO)  GO TO 5001
C
C           NORMAL CASE -- D1 AND D2 BOTH NONZERO, OPPOSITE SIGNS.
C
            NU = ONE - LAMBDA - TWO*RHO
            SIGMA = ONE - RHO
            CP = NU + SIGMA
            IF (ABS(NU) .GT. SMALL)  THEN
               RADCAL = (NU - (TWO*RHO+ONE))*NU + SIGMA**2
               IF (RADCAL .LT. ZERO)  GO TO 5002
               THAT = (CP - SQRT(RADCAL)) / (THREE*NU)
            ELSE
               THAT = ONE/(TWO*SIGMA)
            ENDIF
         ENDIF
         PHI = THAT*((NU*THAT - CP)*THAT + ONE)
C
C          CONVERT TO DISTANCE FROM F2 IF IEXTRM.NE.1 .
         IF (IEXTRM .NE. 1)  PHI = PHI - RHO
C
C          TEST FOR EXCEEDING LIMIT, AND ADJUST ACCORDINGLY.
         DMAX = DFMAX / (H*ABS(PHI))
         IF (ABS(D1) .GT. DMAX)  THEN
            D1 = SIGN (DMAX, D1)
            D2 = -LAMBDA*D1
         ENDIF
      ENDIF
C
C  NORMAL RETURN.
C
 5000 CONTINUE
      IERR = 0
      RETURN
C
C  ERROR RETURNS.
C
 5001 CONTINUE
C     D1 AND D2 BOTH ZERO, OR BOTH NONZERO AND SAME SIGN.
      IERR = -1
      CALL XERROR ('PCHSW -- D1 AND/OR D2 INVALID'
     *           , 29, IERR, 1)
      RETURN
C
 5002 CONTINUE
C     NEGATIVE VALUE OF RADICAL (SHOULD NEVER OCCUR).
      IERR = -2
      CALL XERROR ('PCHSW -- NEGATIVE RADICAL'
     *           , 25, IERR, 1)
      RETURN
C------------- LAST LINE OF PCHSW FOLLOWS ------------------------------
      END
