      SUBROUTINE DE1(F,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,
     1  YY,WT,P,YP,YPOUT,PHI,ALPHA,BETA,SIG,V,W,G,PHASE1,PSI,X,H,HOLD,
     2  START,TOLD,DELSGN,NS,NORND,K,KOLD,ISNOLD)
C
C   ODE  MERELY ALLOCATES STORAGE FOR  DE1 TO RELIEVE THE USER OF THE
C   INCONVENIENCE OF A LONG CALL LIST.  CONSEQUENTLY  DE1 IS USED AS
C   DESCRIBED IN THE COMMENTS FOR  ODE .
C
C   THIS CODE IS COMPLETELY EXPLAINED AND DOCUMENTED IN THE TEXT,
C   COMPUTER SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS0  THE INITIAL
C   VALUE PROBLEM  BY L. F. SHAMPINE AND M. K. GORDON.
C
      LOGICAL STIFF,CRASH,START,PHASE1,NORND
      DIMENSION Y(NEQN),YY(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),
     1  YPOUT(NEQN),PSI(12),ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13)
      EXTERNAL F
C
C   THE CONSTANT  MAXNUM  IS THE MAXIMUM NUMBER OF STEPS ALLOWED IN ONE
C   CALL TO  DE1.  THE USER MAY CHANGE THIS LIMIT BY ALTERING THE
C   FOLLOWING STATEMENT
      DATA MAXNUM/500/
C
C   ****** U IS A MACHINE DEPENDENT PARAMETER. IT IS THE SMALLEST
C          POSITIVE NUMBER FOR WHICH 1.0 + U .GT. 1.0.
C
      U = SPMPAR(1)
      FOURU = 4.0*U
C
C            ***            ***            ***
C   TEST FOR IMPROPER PARAMETERS
C
      IF(NEQN .LT. 1) GO TO 10
      IF(T .EQ. TOUT) GO TO 10
      IF(RELERR .LT. 0.0  .OR.  ABSERR .LT. 0.0) GO TO 10
      EPS = AMAX1(RELERR,ABSERR)
      IF(EPS .LE. 0.0) GO TO 10
      IF(IFLAG .EQ. 0) GO TO 10
      ISN = ISIGN(1,IFLAG)
      IFLAG = IABS(IFLAG)
      IF(IFLAG .EQ. 1) GO TO 20
      IF(T .NE. TOLD) GO TO 10
      IF(IFLAG .GE. 2  .AND.  IFLAG .LE. 5) GO TO 20
      IF(IFLAG .EQ. 6  .AND.  ABSERR .GT. 0.0) GO TO 20
 10   IFLAG = 7
      RETURN
C
C   ON EACH CALL SET INTERVAL OF INTEGRATION AND COUNTER FOR NUMBER OF
C   STEPS.  ADJUST INPUT ERROR TOLERANCES TO DEFINE WEIGHT VECTOR FOR
C   SUBROUTINE  STEP1
C
 20   DEL = TOUT - T
      ABSDEL = ABS(DEL)
      TEND = T + 10.0*DEL
      IF(ISN .LT. 0) TEND = TOUT
      NOSTEP = 0
      KLE4 = 0
      STIFF = .FALSE.
      RELEPS = RELERR/EPS
      ABSEPS = ABSERR/EPS
      IF(IFLAG .EQ. 1) GO TO 30
      IF(ISNOLD .LT. 0) GO TO 30
      IF(DELSGN*DEL .GT. 0.0) GO TO 50
C
C   ON START AND RESTART ALSO SET WORK VARIABLES X AND YY(*), STORE THE
C   DIRECTION OF INTEGRATION AND INITIALIZE THE STEP SIZE
C
 30   START = .TRUE.
      X = T
      DO 40 L = 1,NEQN
 40     YY(L) = Y(L)
      DELSGN = SIGN(1.0,DEL)
      H = SIGN(AMAX1(FOURU*ABS(X),ABS(TOUT-X)),TOUT-X)
C
C   IF ALREADY PAST OUTPUT POINT, INTERPOLATE AND RETURN
C
 50   IF(ABS(X-T) .LT. ABSDEL) GO TO 60
      CALL INTRP(X,YY,TOUT,Y,YPOUT,NEQN,KOLD,PHI,PSI)
      IFLAG = 2
      T = TOUT
      TOLD = T
      ISNOLD = ISN
      RETURN
C
C   IF CANNOT GO PAST OUTPUT POINT AND SUFFICIENTLY CLOSE,
C   EXTRAPOLATE AND RETURN
C
 60   IF(ISN .GT. 0  .OR.  ABS(TOUT-X) .GE. FOURU*ABS(X)) GO TO 80
      H = TOUT - X
      CALL F(X,YY,YP)
      DO 70 L = 1,NEQN
 70     Y(L) = YY(L) + H*YP(L)
      IFLAG = 2
      T = TOUT
      TOLD = T
      ISNOLD = ISN
      RETURN
C
C   TEST FOR TOO MANY STEPS
C
 80   IF(NOSTEP .LT. MAXNUM) GO TO 100
      IFLAG = ISN*4
      IF(STIFF) IFLAG = ISN*5
      DO 90 L = 1,NEQN
 90     Y(L) = YY(L)
      T = X
      TOLD = T
      ISNOLD = 1
      RETURN
C
C   LIMIT STEP SIZE, SET WEIGHT VECTOR AND TAKE A STEP
C
 100  H = SIGN(AMIN1(ABS(H),ABS(TEND-X)),H)
      DO 110 L = 1,NEQN
        WT(L) = RELEPS*ABS(YY(L)) + ABSEPS
        IF(WT(L) .LE. 0.0) GO TO 140
 110    CONTINUE
      CALL STEP1(F,NEQN,YY,X,H,EPS,WT,START,
     1  HOLD,K,KOLD,CRASH,PHI,P,YP,PSI,
     2  ALPHA,BETA,SIG,V,W,G,PHASE1,NS,NORND)
C
C   TEST FOR TOLERANCES TOO SMALL
C
      IF(.NOT.CRASH) GO TO 130
      IFLAG = ISN*3
      RELERR = EPS*RELEPS
      ABSERR = EPS*ABSEPS
      DO 120 L = 1,NEQN
 120    Y(L) = YY(L)
      T = X
      TOLD = T
      ISNOLD = 1
      RETURN
C
C   AUGMENT COUNTER ON NUMBER OF STEPS AND TEST FOR STIFFNESS
C
 130  NOSTEP = NOSTEP + 1
      KLE4 = KLE4 + 1
      IF(KOLD .GT. 4) KLE4 = 0
      IF(KLE4 .GE. 50) STIFF = .TRUE.
      GO TO 50
C
C   RELATIVE ERROR CRITERION INAPPROPRIATE
C
 140  IFLAG = ISN*6
      DO 150 L = 1,NEQN
 150    Y(L) = YY(L)
      T = X
      TOLD = T
      ISNOLD = 1
      RETURN
      END
