      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
