      SUBROUTINE HRM3CF(C,NPNTS,MCFLAG)
C
C*********************************************************************
C
C   NAME- HRM3CF
C
C   LANGUAGE- VS FORTRAN       TYPE- SUBROUTINE    COMPUTER- HDS9000
C
C   VERSION- 1.0       DATE- 12/12/91   PROGRAMMER- CARLA NEIGH, SMSRC
C   VERSION- 3.0       DATE- 06/04/96   PROGRAMMER- TIM GARDNER (HSTX)
C
C   FUNCTION- THIS SUBROUTINE COMPUTES THE SPLINE COEFFICIENTS FOR A
C             PIECEWISE CUBIC HERMITE POLYNOMIAL FUNCTIONAL
C             REPRESENTATION OF A GIVEN NUMBER OF DATA POINTS.
C
C   VERSION 3.0- DOCUMENTATION AND SYNTAX CHANGES FOR ATOVS BASELINE
C                ON THE CRAY.
C
C   CALLING SEQUENCE PARAMETERS-
C   NAME   TYPE   I/O   DESCRIPTION
C   ----   ----   ---   -----------
C   C      REAL   I/O   MATRIX (DIM 4 X NPNTS) CONTAINING ORDINATES
C                       IN THE FIRST ROW, THE FIRST DERIVATIVES
C                       IN THE SECOND ROW, AND THE DIFFERENCES BETWEEN
C                       ADJACENT ABCISSA VALUES IN THE THIRD ROW (I.E.
C                       C(3,I) = X(I) - X(I-1)) ON INPUT; ON OUTPUT
C                       MATRIX CONTAINS VALUES OF THE VARIOUS
C                       ORDER DERIVATIVES AT THE DATA POINTS
C   NPNTS  INT     I    NUMBER OF DATA POINTS
C   MCFLAG LOG     I    MONOTONICITY FLAG-
C                       IF SET TO .TRUE. THEN THE CALCULATED SLOPES
C                       ARE CONSTRAINED PRIOR TO THE CALCULATION OF
C                       THE HIGHER ORDER DERIVATIVES AT THE SUPPLIED
C                       DATA POINTS.  THIS STEP PRESERVES LOCAL
C                       MONOTONICITY PROPERTIES (E.G. SHAPE) OF THE
C                       INTERPOLANT IN REGIONS WHERE THE GRID SPACING
C                       IS COURSE COMPARED WITH UNCONSTRAINED
C                       PIECEWISE POLYNOMIAL VARIATION, ELIMNIATING
C                       UNDESIRABLE OSCILLATIONS.  THE CONSTRAINING
C                       ALGORITHM IS TAKEN FROM THE HYMAN, JAMES M.,
C                       1983:  ACCURATE MONOTONICITY PRESERVING CUBIC
C                       INTERPOLATION, SIAM J. SCI. STAT. COMPUT., 4,
C                       645-654.  IT IS RECOMMENDED THAT MCFLAG BE
C                       SET TO .TRUE. FOR MOST APPLICATIONS, AND SET
C                       TO .FALSE. ONLY IF ANY RESULTING ISOLATED
C                       DISCONTINUITIES IN THE SECOND DERIVATIVE ARE
C                       NOT DESIRED
C
C   LOCAL VARIABLES-
C   NAME    TYPE   DESCRIPTION
C   ----    ----   -----------
C   I       INT    DO LOOP VARIABLE
C   DX      REAL   DIFFERENCE BETWEEN ADJACENT X VALUES
C   DIVDF1  REAL   DIFFERENCE BETWEEN ADJACENT ORDINATES DIVIDED BY DX
C   DIVDF3  REAL   DIFFERENCE BETWEEN ADJACENT FIRST DERIVATIVES
C                  MINUS 2 TIMES DIVDF1
C   CNSTRN  REAL   CONSTRAINT APPLIED WHEN PRESERVING MONOTONICITY
C
C   FUNCTIONS USED- AMIN1, AMAX1, ABS
C
C   EXIT STATES- NONE
C
C   COMMENTS-
C     SEE THE FOLLOWING REFERENCES:
C       DE BOOR, CARL,1978: A PRACTICAL GUIDE TO SPLINES, SPRINGER-
C          VERLAG, 392 PP.
C       HYMAN, JAMES M.,1983: ACCURATE MONOTONICITY PRESERVING CUBIC
C          INTERPOLATION, SIAM J.SCI.STAT.COMPUT., 4, 645-654.
C
C*********************************************************************
C
      INTEGER NPNTS, I
      REAL C(4,NPNTS), DX, DIVDF1, DIVDF3, CNSTRN
      LOGICAL MCFLAG
C
C  SET INITIAL VALUES FOR THE 4TH ROW OF THE COEFFICIENT MATRIX.
C
      DO 10 I = 2, NPNTS
        C(4,I) = (C(1,I) - C(1,I-1)) / C(3,I)
   10 CONTINUE
C
C  PRESERVE MONOTONICITY, IF DESIRED
C
      IF (MCFLAG) THEN
        DO 20 I = 2, NPNTS-1
          CNSTRN = 3.0 * AMIN1(ABS(C(4,I)),ABS(C(4,I+1)))
          IF (C(2,I) .GE. 0.0) THEN
             C(2,I) = AMIN1(C(2,I), CNSTRN)
          ELSE
             C(2,I) = AMAX1(C(2,I), -CNSTRN)
          ENDIF
   20   CONTINUE
        C(2,1) = ((2.0 * C(3,2) + C(3,3)) * C(4,2) - C(3,2) * C(4,3))
     1            / (C(3,2) + C(3,3))
        C(2,NPNTS) = ((2.0 * C(3,NPNTS) + C(3,NPNTS-1)) * C(4,NPNTS)
     1               - C(3,NPNTS) * C(4,NPNTS-1)) / (C(3,NPNTS)
     2               + C(3,NPNTS-1))
      ENDIF
C
C GENERATE THE CUBIC HERMITE PIECEWISE-POLYNOMIALS
C
      DO 30 I = 2, NPNTS
        DX = C(3,I)
        DIVDF1 = C(4,I)
        DIVDF3 = C(2,I-1) + C(2,I) - 2.0 * DIVDF1
        C(3,I-1) = 2.0 * (DIVDF1 - C(2,I-1) - DIVDF3) / DX
        C(4,I-1) = (DIVDF3 / DX) * (6.0 / DX)
   30 CONTINUE
      DO 40 I = 2,4
        C(I,NPNTS) = 0.0
   40 CONTINUE
C
      RETURN
      END
