      SUBROUTINE SPLNCF(C,NPNTS,MCFLAG)
C
C*******************************************************************
C   NAME- SPLNCF
C
C   LANGUAGE- VS FORTRAN    TYPE- SUBROUTINE   COMPUTER- HDS9000
C
C   VERSION- 1.0     DATE- 12/17/91    PROGRAMMER- CARLA NEIGH(SMSRC)
C   VERSION- 3.0     DATE- 06/04/96    PROGRAMMER- TIM GARDNER (HSTX)
C
C   FUNCTIONS- THIS SUBROUTINE COMPUTES SPLINE COEFFICIENTS FOR A
C              PIECEWISE CUBIC POLYNOMIAL FUNCTIONS REPRESENTATION
C              OF A GIVEN NUMBER OF DATA POINTS.  THE 4 COEFFICIENTS
C              (C(1,I) - C(4,I)) FOR EACH DATA POINT ARE RETURNED SO
C              THAT THE ITH POLYNOMIAL PIECE P(X,I) MAY BE WRITTEN
C              AS:
C                P(X,I) = C(1,I) + C(2,I) (X - X(I))
C                         + (1/2) C(3,I) (X - X(I))**2
C                         + (1/6) C(4,I) (X - X(I))**3
C
C              THE FIRST DERIVATIVES AT EACH DATA POINT ARE CALCULATED
C              USING GAUSSIAN ELIMINATION, AND THEN A CUBIC HERMITE
C              POLYNOMIAL REPRESENTATION IS OBTAINED USING THE ROUTINE
C              HRM3CF.
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, AND THE DIFFERENCES BETWEEN
C                       ADJACENT ABCISSAE VALUES IN THE THIRD ROW
C                       (I.E. C(3,I) = X(I) - X(I-1)) ON INPUT; ON
C                       OUTPUT THE MATRIX CONTAINS VALUES OF THE
C                       VARIOUS 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, ELMINIATING
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   G       REAL  VALUE USED IN GAUSSIAN ELIMINATION
C
C   SUBPROGRAMS CALLED-  HRM3CF
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), G
      LOGICAL MCFLAG
C
C  COMPUTE THE FIRST DIVIDED DIFFERENCE OF DATA AND STORE IN C(4,.)
C
      DO 10 I = 2, NPNTS
        C(4,I) = (C(1,I) - C(1,I-1)) / C(3,I)
   10 CONTINUE
C
C  NOT-A-KNOT CONDITION, I.E., FIT A CUBIC POLYNOMIAL TO THE FIRST
C  THREE POINTS
C
      C(4,1) = C(3,3)
      C(3,1) = C(3,2) + C(3,3)
      C(2,1) = ((C(3,2) + 2.0 * C(3,1)) * C(4,2) * C(3,3) + C(3,2) *
     1         C(3,2) * C(4,3)) / C(3,1)
C
C  GENERATE EQUATIONS FOR INTERIOR POINTS AND CARRY OUT THE FORWARD
C  PASS OF THE GAUSS ELIMINATION
C
      DO 20 I = 2, NPNTS-1
        G =  - C(3,I+1) / C(4,I-1)
        C(2,I) = G * C(2,I-1) + 3.0 * (C(3,I) * C(4,I+1) +
     1           C(3,I+1) * C(4,I))
        C(4,I) = G * C(3,I-1) + 2.0 *(C(3,I) + C(3,I+1))
   20 CONTINUE
C
C  NOT-A-KNOT CONDITION, I.E., FIT A CUBIC POLYNOMIAL TO THE LAST
C  THREE POINTS
C
      G = C(3,NPNTS-1) + C(3,NPNTS)
      C(2,NPNTS) = ((C(3,NPNTS) + 2.0 * G) * C(4,NPNTS) * C(3,NPNTS-1)
     1        + C(3,NPNTS) * C(3,NPNTS) * (C(1,NPNTS-1) - C(1,NPNTS-2))
     2        / C(3,NPNTS-1)) / G
      G = - G/C(4,NPNTS-1)
      C(4,NPNTS) = C(3,NPNTS-1)
C
C  COMPLETE THE FORWARD PASS OF THE GAUSS ELIMINATION
C
      C(4,NPNTS) = G * C(3,NPNTS-1) + C(4,NPNTS)
      C(2,NPNTS) = (G * C(2,NPNTS-1) + C(2,NPNTS)) / C(4,NPNTS)
C
C  CARRY OUT BACK SUBSTITUTION
C
      DO 30 I = NPNTS-1, 1, -1
        C(2,I) = (C(2,I) - C(3,I) * C(2,I+1)) / C(4,I)
   30 CONTINUE
C
C  GENERATE THE CUBIC SPLINE COEFFICIENTS IN EACH INTERVAL
C
      CALL HRM3CF(C,NPNTS,MCFLAG)
C
      RETURN
      END
