      SUBROUTINE DRIVE(N,T0,H0,Y0,TOUT,AERR,RERR,NUM,MF,INDEX,
     1     Y,YMAX,ERROR,SAVE1,SAVE2,PW,IPIV,DIFFUN,PEDERV)
C-----------------------------------------------------------------------
C THIS IS A MODIFICATION OF THE JUNE 24, 1975 VERSION OF EPISODE..
C A PACKAGE FOR SOLVING SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS
C    DY/DT = F(Y,T),  Y = (Y(1),Y(2),...,Y(N)) TRANSPOSE,
C GIVEN THE INITIAL VALUE OF Y.
C DEVELOPED BY G.D.BYRNE (UNIV OF PITTSBURGH) AND A.C.HINDMARSH
C (LAWRENCE LIVERMORE LAB). MODIFIED BY A.H.MORRIS (NSWC/DL).
C                 REFERENCES
C 1.  G. D. BYRNE AND A. C. HINDMARSH, A POLYALGORITHM FOR THE
C       NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS,
C       UCRL-75652, LAWRENCE LIVERMORE LABORATORY, P. O. BOX 808,
C       LIVERMORE, CA 94550, APRIL 1974. ALSO IN ACM TRANSACTIONS
C       ON MATHEMATICAL SOFTWARE, 1 (1975), PP. 71-96.
C
C 2.  A. C. HINDMARSH AND G. D. BYRNE, EPISODE.. AN EXPERIMENTAL
C       PACKAGE FOR THE INTEGRATION OF SYSTEMS OF ORDINARY
C       DIFFERENTIAL EQUATIONS, UCID-30112, L.L.L., MAY, 1975.
C
C 3.  A. C. HINDMARSH, GEAR.. ORDINARY DIFFERENTIAL EQUATION
C       SYSTEM SOLVER, UCID-30001, REV. 3, L.L.L., DECEMBER, 1974.
C
C-----------------------------------------------------------------------
C DRIVE IS A DRIVER SUBROUTINE FOR THE EPISODE PACKAGE.
C DRIVE IS TO BE CALLED ONCE FOR EACH OUTPUT VALUE OF T.
C IT THEN MAKES REPEATED CALLS TO THE CORE INTEGRATOR
C SUBROUTINE, TSTEP.
C
C THE INPUT PARAMETERS ARE AS FOLLOWS.
C   N     =  THE NUMBER OF DIFFERENTIAL EQUATIONS. N MAY BE MODIFIED
C              WHEN INDEX=-1. N MAY BE DECREASED, BUT NEVER INCREASED
C              DURING A GIVEN PROBLEM.
C   T0    =  THE INITIAL VALUE OF T, THE INDEPENDENT VARIABLE
C              (USED FOR INPUT ONLY ON FIRST CALL).
C   H0    =  THE STEP SIZE IN T (USED FOR INPUT ONLY ON THE
C              FIRST CALL).
C   Y0    =  A VECTOR OF LENGTH N CONTAINING THE INITIAL VALUES OF
C              Y (USED FOR INPUT ONLY ON FIRST CALL).
C   TOUT  =  THE VALUE OF T AT WHICH OUTPUT IS DESIRED NEXT.
C              INTEGRATION WILL NORMALLY GO BEYOND TOUT AND
C              INTERPOLATE TO T = TOUT.  (USED ONLY FOR INPUT.)
C   AERR  =  THE ABSOLUTE ERROR TOLERANCE TO BE SATISFIED
C   RERR  =  THE RELATIVE ERROR TOLERANCE TO BE SATISFIED
C   NUM   =  THE MAXIMUM NUMBER OF STEPS TO BE PERFORMED. IF
C              NUM=1 THEN INTEGRATION WILL NOT GO BEYOND TOUT.
C   MF    =  THE METHOD FLAG (USED ONLY ON FIRST CALL, UNLESS
C              INDEX = -1).  ALLOWED VALUES ARE 10, 11, 12, 13,
C              20, 21, 22, 23.  MF IS AN INTEGER WITH TWO DECIMAL
C              DIGITS, METH AND MITER (MF = 10*METH + MITER).  (MF
C              CAN BE THOUGHT OF AS THE ORDERED PAIR (METH,MITER).)
C              METH IS THE BASIC METHOD INDICATOR.
C                METH = 1 INDICATES VARIABLE-STEP SIZE, VARIABLE-
C                         ORDER ADAMS METHOD, SUITABLE FOR NON-
C                         STIFF PROBLEMS.
C                METH = 2 INDICATES VARIABLE-STEP SIZE, VARIABLE-
C                         ORDER BACKWARD DIFFERENTIATION METHOD,
C                         SUITABLE FOR STIFF PROBLEMS.
C              MITER INDICATES THE METHOD OF ITERATIVE CORRECTION
C                (NONLINEAR SYSTEM SOLUTION).
C                MITER = 0 INDICATES FUNCTIONAL ITERATION (NO
C                          PARTIAL DERIVATIVES NEEDED).
C                MITER = 1 INDICATES A CHORD OR SEMI-STATIONARY
C                          NEWTON METHOD WITH CLOSED FORM (EXACT)
C                          JACOBIAN, WHICH IS COMPUTED IN THE
C                          USER SUPPLIED SUBROUTINE
C                          PEDERV(N,T,Y,PD,N0) DESCRIBED BELOW.
C                MITER = 2 INDICATES A CHORD  OR SEMI-STATIONARY
C                          NEWTON METHOD WITH AN INTERNALLY
C                          COMPUTED FINITE DIFFERENCE APPROXIMATION
C                          TO THE JACOBIAN.
C                MITER = 3 INDICATES A CHORD OR SEMI-STATIONARY
C                          NEWTON METHOD WITH AN INTERNALLY
C                          COMPUTED DIAGONAL MATRIX APPROXIMATION
C                          TO THE JACOBIAN, BASED ON A DIRECTIONAL
C                          DERIVATIVE.
C   INDEX =  INTEGER USED ON INPUT TO INDICATE TYPE OF CALL,
C              WITH THE FOLLOWING VALUES AND MEANINGS..
C         1    THIS IS THE FIRST CALL FOR THIS PROBLEM.
C         0    THIS IS NOT THE FIRST CALL FOR THIS PROBLEM,
C              AND INTEGRATION IS TO CONTINUE.
C        -1    THIS IS NOT THE FIRST CALL FOR THE PROBLEM.
C              THE USER HAS RESET N, AERR, RERR, AND/OR MF.
C         2    SAME AS 0 EXCEPT THAT TOUT IS TO BE HIT
C              EXACTLY (NO INTERPOLATION IS DONE). THIS
C              ASSUMES TOUT IS BEYOND THE CURRENT T.
C            SINCE THE NORMAL OUTPUT VALUE OF INDEX IS 0,
C            IT NEED NOT BE RESET FOR NORMAL CONTINUATION.
C THE REMAINING INPUT PARAMETERS ARE ARRAYS THAT ARE USED AS WORK
C SPACES. THESE ARRAYS HAVE THE FOLLOWING DIMENSIONS..
C     DIMENSION Y(N,13)
C     DIMENSION YMAX(N),ERROR(N),SAVE1(N),SAVE2(N)
C     DIMENSION PW(N*N),IPIV(N)
C THE ROW LENGTH OF Y CAN BE REDUCED FROM 13 TO 6 IF METH = 2.
C THE ARRAY IPIV IS USED IF AND ONLY IF MITER = 1 OR MITER = 2.
C THE SIZE OF THE PW ARRAY CAN BE REDUCED TO 1 IF MITER = 0 OR
C TO N IF MITER = 3.
C
C AFTER THE INITIAL CALL, IF A NORMAL RETURN OCCURRED AND A NORMAL
C CONTINUATION IS DESIRED, SIMPLY RESET TOUT AND CALL AGAIN.
C ALL OTHER PARAMETERS WILL BE READY FOR THE NEXT CALL.
C A CHANGE OF PARAMETERS WITH INDEX = -1 CAN BE MADE AFTER
C EITHER A SUCCESSFUL OR AN UNSUCCESSFUL RETURN.
C
C THE OUTPUT PARAMETERS ARE AS FOLLOWS.
C   T0    =  THE FURTHEST POINT REACHED IN THE INTEGRATION PROCESS.
C   H0    =  THE STEP SIZE H USED LAST, WHETHER SUCCESSFULLY OR NOT.
C   Y0    =  THE COMPUTED VALUE OF Y AT TOUT IF THE INTEGRATION WAS
C            SUCCESSFUL. OTHERWISE Y0 IS THE COMPUTED VALUE OF Y AT
C            THE (OUTPUT) POINT T0.
C   INDEX =  INTEGER USED ON OUTPUT TO INDICATE RESULTS,
C              WITH THE FOLLOWING VALUES AND MEANINGS..
C         0    INTEGRATION WAS COMPLETED TO TOUT OR BEYOND.
C        -1    THE INTEGRATION WAS HALTED AFTER FAILING TO PASS THE
C              ERROR TEST EVEN AFTER REDUCING H BY A FACTOR OF
C              1.E10 FROM ITS INITIAL VALUE.
C        -2    AFTER SOME INITIAL SUCCESS, THE INTEGRATION WAS
C              HALTED EITHER BY REPEATED ERROR TEST FAILURES OR
C              BY A TEST ON EPS.  POSSIBLY TOO MUCH ACCURACY HAS
C              BEEN REQUESTED, OR A BAD CHOICE OF MF WAS MADE.
C        -3    THE INTEGRATION WAS HALTED AFTER FAILING TO ACHIEVE
C              CORRECTOR CONVERGENCE EVEN AFTER REDUCING H BY A
C              FACTOR OF 1.E10 FROM ITS INITIAL VALUE.
C        -4    IMMEDIATE HALT BECAUSE OF ILLEGAL VALUES OF INPUT
C              PARAMETERS.
C        -5    INDEX WAS -1 ON INPUT, BUT THE DESIRED CHANGES OF
C              PARAMETERS WERE NOT IMPLEMENTED BECAUSE TOUT
C              WAS NOT BEYOND T.  INTERPOLATION TO TOUT WAS
C              PERFORMED AS ON A NORMAL RETURN.  TO CONTINUE,
C              SIMPLY CALL AGAIN WITH INDEX = -1 AND A NEW TOUT.
C        -6    INDEX WAS 2 OR NUM WAS 1 ON INPUT, BUT TOUT WAS
C              NOT BEYOND T. NO ACTION WAS TAKEN.
C        -7    IMMEDIATE HALT SINCE AERR = 0 AND SOME INPUT
C              Y0(I) = 0. AERR MUST BE MADE POSITIVE.
C        -8    THE INTEGRATION WAS HALTED BECAUSE AERR = 0
C              AND SOME Y0(I) = 0. AERR MUST BE MADE POSITIVE.
C        -9    NUM STEPS HAVE BEEN PERFORMED.
C
C IN ADDITION TO DRIVE, THE FOLLOWING SUBROUTINES ARE USED BY AND
C PROVIDED IN THIS PACKAGE
C   INTERP(TOUT,Y,N0,Y0)  INTERPOLATES TO GIVE OUTPUT VALUES AT
C                         T = TOUT  BY USING DATA IN THE Y ARRAY.
C   TSTEP(----)  IS THE CORE INTEGRATION SUBROUTINE, WHICH INTEGRATES
C                OVER A SINGLE STEP AND DOES ASSOCIATED ERROR
C                CONTROL.
C   COSET  SETS COEFFICIENTS FOR USE IN TSTEP.
C   ADJUST(Y,N0)  ADJUSTS THE HISTORY ARRAY Y ON REDUCTION OF ORDER.
C   PSET(----)   COMPUTES AND PROCESSES THE JACOBIAN MATRIX DF/DY.
C   DEC(N,N0,A,IP,IER)  PERFORMS THE LU DECOMPOSITION OF A MATRIX.
C   SOL(N,N0,A,B,IP)  SOLVES A LINEAR SYSTEM A*X = B, AFTER DEC
C                       HAS BEEN CALLED FOR THE MATRIX A.
C NOTE   PSET, DEC, AND SOL ARE CALLED IF AND ONLY IF MITER = 1
C        OR MITER = 2.
C
C THE USER MUST FURNISH THE FOLLOWING SUBROUTINES
C   DIFFUN(N,T,Y,YDOT)  COMPUTES THE FUNCTION  YDOT = F(Y,T),
C                       THE RIGHT HAND SIDE OF THE ORDINARY
C                       DIFFERENTIAL EQUATION SYSTEM, WHERE Y
C                       AND YDOT ARE VECTORS OF LENGTH N.
C   PEDERV(N,T,Y,PD,N0)  COMPUTES THE N BY N JACOBIAN MATRIX OF
C                        PARTIAL DERIVATIVES AND STORES IT IN PD AS
C                        AN N0 BY N0 ARRAY.   PD(I,J) IS TO BE SET
C                        TO THE PARTIAL DERIVATIVE OF YDOT(I) WITH
C                        RESPECT TO Y(J).  PEDERV IS CALLED IF AND
C                        ONLY IF MITER = 1. FOR OTHER VALUES OF
C                        MITER, PEDERV CAN BE A DUMMY SUBROUTINE.
C
C THE COMMON BLOCK EPCOM9 CAN BE ACCESSED EXTERNALLY BY THE USER,
C IF HE DESIRES.  IT CONTAINS THE STEP SIZE LAST USED SUCCESSFULLY
C (HUSED), THE ORDER LAST USED SUCCESSFULLY (NQUSED), THE
C NUMBER OF STEPS TAKEN SO FAR (NSTEP), THE NUMBER OF FUNCTION
C EVALUATIONS (DIFFUN CALLS) SO FAR (NFE), AND THE NUMBER OF
C JACOBIAN EVALUATIONS SO FAR (NJE).
C
C-----------------------------------------------------------------------
C
      INTEGER NUM, INDEX, MF, N
      INTEGER IPIV, JSTART, KFLAG, MFC, NC, NFE, NJE,
     1        NQUSED, NSQ, NSTEP
      INTEGER M, I, KGO, NHCUT, N0
      INTEGER IDUMMY, IWK
      REAL AERR, RERR, H0, TOUT, T0, Y0
      REAL EPS, AEPS, REPS, TPH
      REAL EPSC, EPSJ, ERROR, HMAX, H, HMIN, HUSED,
     1     PW, SAVE1, SAVE2, SS, T, UROUND, YMAX
      REAL AYI, D, Y
      REAL DUMMY, PRL1, RC, TAU, EL, TQ
      REAL HCUT
      REAL FOUR, HUNDRD, ONE, TEN, ZERO
      DIMENSION Y0(N),Y(N,13)
      DIMENSION YMAX(1),ERROR(1),SAVE1(1),SAVE2(1),PW(1),IPIV(1)
      REAL SPMPAR
      EXTERNAL DIFFUN,PEDERV
C
      COMMON /XXXDRV/ AEPS,REPS,NHCUT
      COMMON /XXXSTP/ DUMMY(7),IDUMMY(6),PRL1,RC
      COMMON /EPCOM1/T,H,HMIN,HMAX,EPSC,SS,UROUND,NC,MFC,KFLAG,JSTART
      COMMON /EPCOM8/ EPSJ,NSQ
      COMMON /EPCOM9/ HUSED,NQUSED,NSTEP,NFE,NJE
      COMMON /EPCM10/ TAU(13),EL(13),TQ(5),IWK(5)
C
      DATA HCUT /0.1E0/
      DATA FOUR /4.0E0/, HUNDRD /1.0E2/, ONE  /1.0E0/,
     1     TEN  /1.0E1/, ZERO   /0.0E0/
C
C ****** UROUND IS A MACHINE DEPENDENT PARAMETER. ASSIGN UROUND
C        THE VALUE U WHERE U IS THE SMALLEST POSITIVE FLOATING
C        POINT NUMBER FOR WHICH 1.0 + U .GT. 1.0.
C
      UROUND = SPMPAR(1)
C
C
      M = 0
      IF (NUM .LE. 0) GO TO 400
      IF (INDEX .EQ. 0) GO TO 20
      IF (INDEX .EQ. 2) GO TO 25
      IF (AERR .LT. ZERO .OR. RERR .LT. ZERO) GO TO 400
      EPS = AMAX1(AERR,RERR)
      IF (EPS .LE. ZERO .OR. N .LE. 0) GO TO 400
      AEPS = AERR/EPS
      REPS = RERR/EPS
      IF (INDEX .EQ. -1) GO TO 30
      IF (INDEX .NE. 1 .OR. (T0-TOUT)*H0 .GE. ZERO) GO TO 400
C-----------------------------------------------------------------------
C IF INITIAL VALUES FOR YMAX OTHER THAN THOSE BELOW ARE DESIRED,
C THEY SHOULD BE SET HERE.  ALL YMAX(I) MUST BE POSITIVE.  IF
C VALUES FOR HMIN OR HMAX, THE BOUNDS ON THE ABSOLUTE VALUE OF H,
C OTHER THAN THOSE BELOW, ARE DESIRED, THEY ALSO SHOULD BE SET HERE.
C IF ERROR PER SS UNITS OF T IS TO BE CONTROLLED, SS SHOULD BE SET
C TO A POSITIVE VALUE BELOW.  ERROR PER UNIT STEP IS CONTROLLED
C WHEN SS = 1.  THE DEFAULT VALUE FOR SS IS 0 AND YIELDS CONTROL
C OF ERROR PER STEP.
C-----------------------------------------------------------------------
      DO 10 I = 1,N
        YMAX(I) = AEPS+ABS(Y0(I))*REPS
        IF (YMAX(I) .EQ. ZERO) GO TO 300
 10     Y(I,1) = Y0(I)
      NC = N
      T = T0
      HMAX = ABS(T0 - TOUT)
      IF (NUM .NE. 1) HMAX = HMAX*TEN
      HMIN = AMIN1(ABS(H0),HMAX)
      H = SIGN(HMIN,H0)
      EPSC = EPS
      MFC = MF
      JSTART = 0
      SS = ZERO
      N0 = N
      NSQ = N0*N0
      EPSJ = SQRT(UROUND)
      NHCUT = 0
      GO TO 50
C
 20   IF (NUM .EQ. 1) GO TO 25
      HMAX = ABS(T - TOUT)*TEN
      HMIN = AMIN1(HMIN,HMAX)
      GO TO 80
 25   IF ((T-TOUT)*H .GE. ZERO) GO TO 460
      HMAX = ABS(T - TOUT)
      HMIN = AMIN1(HMIN,HMAX)
      GO TO 85
C
 30   IF ((T-TOUT)*H .GE. ZERO) GO TO 450
      IF (MF .NE. MFC) JSTART = -1
      NC = N
      EPSC = EPS
      MFC = MF
      IF (NUM .NE. 1) GO TO 35
      HMAX = ABS(T - TOUT)
      HMIN = AMIN1(HMIN,HMAX)
      GO TO 85
 35   HMAX = ABS(T - TOUT)*TEN
      HMIN = AMIN1(HMIN,HMAX)
C
 50   CALL TSTEP(Y,N0,YMAX,ERROR,SAVE1,SAVE2,PW,IPIV,DIFFUN,PEDERV)
      KGO = 1 - KFLAG
      GO TO (60, 100, 500, 100), KGO
C KFLAG  =   0,  -1,  -2,  -3  -----------------------------------------
C
 60   CONTINUE
C-----------------------------------------------------------------------
C NORMAL RETURN FROM TSTEP.
C
C THE WEIGHTS YMAX(I) ARE UPDATED.  IF DIFFERENT VALUES ARE DESIRED,
C THEY SHOULD BE SET HERE.  IF SS IS TO BE UPDATED FOR CONTROL OF
C ERROR PER SS UNITS OF T, IT SHOULD ALSO BE DONE HERE.  A TEST IS
C MADE TO DETERMINE IF EPS IS TOO SMALL FOR MACHINE PRECISION.
C
C ANY OTHER TESTS OR CALCULATIONS THAT ARE REQUIRED AFTER EACH STEP
C SHOULD BE INSERTED HERE.
C
C IF INDEX = 2, H IS CONTROLLED TO HIT TOUT (WITHIN ROUNDOFF
C ERROR), AND THEN THE CURRENT Y VALUES ARE PUT IN Y0 ON
C RETURN.  FOR ANY OTHER VALUE OF INDEX, CONTROL RETURNS TO
C THE INTEGRATOR UNLESS TOUT HAS BEEN REACHED.  THEN
C INTERPOLATED VALUES OF Y ARE COMPUTED AND STORED IN Y0 ON
C RETURN.
C IF INTERPOLATION IS NOT DESIRED, THE CALL TO INTERP SHOULD
C BE DELETED AND CONTROL TRANSFERRED TO STATEMENT 500 INSTEAD
C OF 520.
C-----------------------------------------------------------------------
      D = ZERO
      DO 70 I = 1,N
        AYI = ABS(Y(I,1))
        YMAX(I) = AEPS+AYI*REPS
        IF (YMAX(I) .EQ. ZERO) GO TO 350
 70     D = D + (AYI/YMAX(I))**2
      D = D*(UROUND/EPS)**2
      IF (D .GT. FLOAT(N)) GO TO 250
      M = M+1
      IF (M .GE. NUM) GO TO 200
      IF (INDEX .EQ. 2) GO TO 85
 80   IF ((T-TOUT)*H .LT. ZERO) GO TO 50
      CALL INTERP (TOUT, Y, N0, Y0)
      T0 = T
      GO TO 520
 85   TPH = T+H
      IF ((TPH-TOUT)*H .LE. ZERO) GO TO 50
      IF (ABS(T-TOUT) .LE. HUNDRD*UROUND*HMAX) GO TO 500
      IF ((T-TOUT)*H .GE. ZERO) GO TO 500
      H = (TOUT - T)*(ONE - FOUR*UROUND)
      JSTART = -1
      GO TO 50
C-----------------------------------------------------------------------
C ON AN ERROR RETURN FROM TSTEP, AN IMMEDIATE RETURN OCCURS IF
C KFLAG = -2, AND RECOVERY ATTEMPTS ARE MADE OTHERWISE.
C TO RECOVER, H AND HMIN ARE REDUCED BY A FACTOR OF .1 UP TO 10
C TIMES BEFORE GIVING UP.
C-----------------------------------------------------------------------
 100  IF (NHCUT .EQ. 10) GO TO 500
      NHCUT = NHCUT + 1
      HMIN = HCUT*HMIN
      H = HCUT*H
      JSTART = -1
      GO TO 50
C
 200  KFLAG = -9
      GO TO 500
C
 250  KFLAG = -2
      GO TO 500
C
 300  INDEX = -7
      RETURN
C
 350  INDEX = -8
      GO TO 500
C
 400  INDEX = -4
      RETURN
C
 450  CALL INTERP (TOUT, Y, N0, Y0)
      T0 = T
      INDEX = -5
      RETURN
C
 460  INDEX = -6
      RETURN
C
 500  T0 = T
      DO 510 I = 1,N
 510    Y0(I) = Y(I,1)
 520  INDEX = KFLAG
      H0 = HUSED
      IF (KFLAG .NE. 0) H0 = H
      RETURN
C-----------------------  END OF SUBROUTINE DRIVE  ---------------------
      END
