      SUBROUTINE CURV1 (N,X,Y,SLP1,SLPN,ISLPSW,YP,TEMP,
     *                  SIGMA,IERR)
C
      INTEGER N,ISLPSW,IERR
      REAL X(N),Y(N),SLP1,SLPN,YP(N),TEMP(N),SIGMA
C
C                      FROM THE SPLINE UNDER TENSION PACKAGE
C                       CODED BY A. K. CLINE AND R. J. RENKA
C                            DEPARTMENT OF COMPUTER SCIENCES
C                              UNIVERSITY OF TEXAS AT AUSTIN
C
C THIS SUBROUTINE DETERMINES THE PARAMETERS NECESSARY TO
C COMPUTE AN INTERPOLATORY SPLINE UNDER TENSION THROUGH
C A SEQUENCE OF FUNCTIONAL VALUES. THE SLOPES AT THE TWO
C ENDS OF THE CURVE MAY BE SPECIFIED OR OMITTED.  FOR ACTUAL
C COMPUTATION OF POINTS ON THE CURVE IT IS NECESSARY TO CALL
C THE FUNCTION CURV2.
C
C ON INPUT--
C
C   N IS THE NUMBER OF VALUES TO BE INTERPOLATED (N.GE.2).
C
C   X IS AN ARRAY OF THE N INCREASING ABSCISSAE OF THE
C   FUNCTIONAL VALUES.
C
C   Y IS AN ARRAY OF THE N ORDINATES OF THE VALUES, (I. E.
C   Y(K) IS THE FUNCTIONAL VALUE CORRESPONDING TO X(K) ).
C
C   SLP1 AND SLPN CONTAIN THE DESIRED VALUES FOR THE FIRST
C   DERIVATIVE OF THE CURVE AT X(1) AND X(N), RESPECTIVELY.
C   THE USER MAY OMIT VALUES FOR EITHER OR BOTH OF THESE
C   PARAMETERS AND SIGNAL THIS WITH ISLPSW.
C
C   ISLPSW CONTAINS A SWITCH INDICATING WHICH SLOPE DATA
C   SHOULD BE USED AND WHICH SHOULD BE ESTIMATED BY THIS
C   SUBROUTINE,
C          = 0 IF SLP1 AND SLPN ARE TO BE USED,
C          = 1 IF SLP1 IS TO BE USED BUT NOT SLPN,
C          = 2 IF SLPN IS TO BE USED BUT NOT SLP1,
C          = 3 IF BOTH SLP1 AND SLPN ARE TO BE ESTIMATED
C              INTERNALLY.
C
C   YP IS AN ARRAY OF LENGTH AT LEAST N.
C
C   TEMP IS AN ARRAY OF LENGTH AT LEAST N WHICH IS USED FOR
C   SCRATCH STORAGE.
C
C AND
C
C   SIGMA CONTAINS THE TENSION FACTOR. THIS VALUE INDICATES
C   THE CURVINESS DESIRED. IF ABS(SIGMA) IS NEARLY ZERO
C   (E.G. .001) THE RESULTING CURVE IS APPROXIMATELY A
C   CUBIC SPLINE. IF ABS(SIGMA) IS LARGE (E.G. 50.) THE
C   RESULTING CURVE IS NEARLY A POLYGONAL LINE. IF SIGMA
C   EQUALS ZERO A CUBIC SPLINE RESULTS.  A STANDARD VALUE
C   FOR SIGMA IS APPROXIMATELY 1. IN ABSOLUTE VALUE.
C
C ON OUTPUT--
C
C   YP CONTAINS THE VALUES OF THE SECOND DERIVATIVE OF THE
C   CURVE AT THE GIVEN NODES.
C
C   IERR CONTAINS AN ERROR FLAG,
C        = 0 FOR NORMAL RETURN,
C        = 1 IF N IS LESS THAN 2,
C        = 2 IF X-VALUES ARE NOT STRICTLY INCREASING.
C
C AND
C
C   N, X, Y, SLP1, SLPN, ISLPSW AND SIGMA ARE UNALTERED.
C
C THIS SUBROUTINE REFERENCES PACKAGE MODULES CEEZ, TERMS,
C AND SNHCSH.
C
C-----------------------------------------------------------
C
      NM1 = N-1
      NP1 = N+1
      IERR = 0
      IF (N .LE. 1) GO TO 8
      IF (X(N) .LE. X(1)) GO TO 9
C
C DENORMALIZE TENSION FACTOR
C
      SIGMAP = ABS(SIGMA)*FLOAT(N-1)/(X(N)-X(1))
C
C APPROXIMATE END SLOPES
C
      IF (ISLPSW .GE. 2) GO TO 1
      SLPP1 = SLP1
      GO TO 2
    1 DELX1 = X(2)-X(1)
      DELX2 = DELX1+DELX1
      IF (N .GT. 2) DELX2 = X(3)-X(1)
      IF (DELX1 .LE. 0. .OR. DELX2 .LE. DELX1) GO TO 9
      CALL CEEZ (DELX1,DELX2,SIGMAP,C1,C2,C3,N)
      SLPP1 = C1*Y(1)+C2*Y(2)
      IF (N .GT. 2) SLPP1 = SLPP1+C3*Y(3)
    2 IF (ISLPSW .EQ. 1 .OR. ISLPSW .EQ. 3) GO TO 3
      SLPPN = SLPN
      GO TO 4
    3 DELXN = X(N)-X(NM1)
      DELXNM = DELXN+DELXN
      IF (N .GT. 2) DELXNM = X(N)-X(N-2)
      IF (DELXN .LE. 0. .OR. DELXNM .LE. DELXN) GO TO 9
      CALL CEEZ (-DELXN,-DELXNM,SIGMAP,C1,C2,C3,N)
      SLPPN = C1*Y(N)+C2*Y(NM1)
      IF (N .GT. 2) SLPPN = SLPPN+C3*Y(N-2)
C
C SET UP RIGHT HAND SIDE AND TRIDIAGONAL SYSTEM FOR YP AND
C PERFORM FORWARD ELIMINATION
C
    4 DELX1 = X(2)-X(1)
      IF (DELX1 .LE. 0.) GO TO 9
      DX1 = (Y(2)-Y(1))/DELX1
      CALL TERMS (DIAG1,SDIAG1,SIGMAP,DELX1)
      YP(1) = (DX1-SLPP1)/DIAG1
      TEMP(1) = SDIAG1/DIAG1
      IF (N .EQ. 2) GO TO 6
      DO 5 I = 2,NM1
        DELX2 = X(I+1)-X(I)
        IF (DELX2 .LE. 0.) GO TO 9
        DX2 = (Y(I+1)-Y(I))/DELX2
        CALL TERMS (DIAG2,SDIAG2,SIGMAP,DELX2)
        DIAG = DIAG1+DIAG2-SDIAG1*TEMP(I-1)
        YP(I) = (DX2-DX1-SDIAG1*YP(I-1))/DIAG
        TEMP(I) = SDIAG2/DIAG
        DX1 = DX2
        DIAG1 = DIAG2
    5   SDIAG1 = SDIAG2
    6 DIAG = DIAG1-SDIAG1*TEMP(NM1)
      YP(N) = (SLPPN-DX1-SDIAG1*YP(NM1))/DIAG
C
C PERFORM BACK SUBSTITUTION
C
      DO 7 I = 2,N
        IBAK = NP1-I
    7   YP(IBAK) = YP(IBAK)-TEMP(IBAK)*YP(IBAK+1)
      RETURN
C
C TOO FEW POINTS
C
    8 IERR = 1
      RETURN
C
C X-VALUES NOT STRICTLY INCREASING
C
    9 IERR = 2
      RETURN
      END
