      SUBROUTINE ODE(F,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,WORK,IWORK)
C
C     SANDIA MATHEMATICAL PROGRAM LIBRARY
C     APPLIED MATHEMATICS DIVISION 2642
C     SANDIA LABORATORIES
C     ALBUQUERQUE, NEW MEXICO  87115
C     CONTROL DATA 6600 VERSION 6.1  JANUARY 1976
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C   *                 ISSUED BY SANDIA LABORATORIES,
C   *                   A PRIME CONTRACTOR TO THE
C   *   UNITED STATES ENERGY RESEARCH AND DEVELOPMENT ADMINISTRATION
C   * * * * * * * * * * * * * *  NOTICE   * * * * * * * * * * * * * * *
C   * THIS REPORT WAS PREPARED AS AN ACCOUNT OF WORK SPONSORED BY THE
C   * UNITED STATES GOVERNMENT.  NEITHER THE UNITED STATES NOR THE
C   * UNITED STATES ENERGY RESEARCH AND DEVELOPMENT ADMINISTRATION,
C   * NOR ANY OF THEIR EMPLOYEES, NOR ANY OF THEIR CONTRACTORS,
C   * SUBCONTRACTORS, OR THEIR EMPLOYEES, MAKES ANY WARRANTY, EXPRESS
C   * OR IMPLIED, OR ASSUMES ANY LEGAL LIABILITY OR RESPONSIBILITY
C   * FOR THE ACCURACY, COMPLETENESS OR USEFULNESS OF ANY INFORMATION,
C   * APPARATUS, PRODUCT OR PROCESS DISCLOSED, OR REPRESENTS THAT ITS
C   * USE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS.
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C   * THE PRIMARY DOCUMENT FOR THE LIBRARY OF WHICH THIS ROUTINE IS
C   * IS A PART IS SAND75-0545, PRINTED FEB 1976.
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C   WRITTEN BY L. F. SHAMPINE AND M. K. GORDON
C
C   ABSTRACT
C
C
C   SUBROUTINE  ODE  INTEGRATES A SYSTEM OF  NEQN  FIRST ORDER
C   ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM
C             DY(I)/DT = F(T,Y(1),Y(2),...,Y(NEQN))
C             Y(I) GIVEN AT  T .
C   THE SUBROUTINE INTEGRATES FROM  T  TO  TOUT .  ON RETURN THE
C   PARAMETERS IN THE CALL LIST ARE SET FOR CONTINUING THE INTEGRATION.
C   THE USER HAS ONLY TO DEFINE A NEW VALUE  TOUT  AND CALL  ODE  AGAIN.
C
C   THE DIFFERENTIAL EQUATIONS ARE ACTUALLY SOLVED BY A SUITE OF CODES
C   DE1,  STEP1, AND  INTRP .  ODE  ALLOCATES VIRTUAL STORAGE IN THE
C   ARRAYS  WORK  AND  IWORK  AND CALLS  DE1.  DE1 IS A SUPERVISOR WHICH
C   DIRECTS THE SOLUTION.  IT CALLS ON THE ROUTINES  STEP1 AND  INTRP
C   TO ADVANCE THE INTEGRATION AND TO INTERPOLATE AT OUTPUT POINTS.
C   STEP1 USES A MODIFIED DIVIDED DIFFERENCE FORM OF THE ADAMS PECE
C   FORMULAS AND LOCAL EXTRAPOLATION.  IT ADJUSTS THE ORDER AND STEP
C   SIZE TO CONTROL THE LOCAL ERROR PER UNIT STEP IN A GENERALIZED
C   SENSE.  NORMALLY EACH CALL TO  STEP1 ADVANCES THE SOLUTION ONE STEP
C   IN THE DIRECTION OF  TOUT .  FOR REASONS OF EFFICIENCY  DE1
C   INTEGRATES BEYOND  TOUT  INTERNALLY, THOUGH NEVER BEYOND
C   T+10*(TOUT-T), AND CALLS  INTRP  TO INTERPOLATE THE SOLUTION AT
C   TOUT .  AN OPTION IS PROVIDED TO STOP THE INTEGRATION AT  TOUT  BUT
C   IT SHOULD BE USED ONLY IF IT IS IMPOSSIBLE TO CONTINUE THE
C   INTEGRATION BEYOND  TOUT .
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
C   THE PARAMETERS REPRESENT0
C      F -- SUBROUTINE F(T,Y,YP) TO EVALUATE DERIVATIVES YP(I)=DY(I)/DT
C      NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED
C      Y(*) -- SOLUTION VECTOR AT  T
C      T -- INDEPENDENT VARIABLE
C      TOUT -- POINT AT WHICH SOLUTION IS DESIRED
C      RELERR,ABSERR -- RELATIVE AND ABSOLUTE ERROR TOLERANCES FOR LOCAL
C           ERROR TEST.  AT EACH STEP THE CODE REQUIRES
C             ABS(LOCAL ERROR) .LE. ABS(Y)*RELERR + ABSERR
C           FOR EACH COMPONENT OF THE LOCAL ERROR AND SOLUTION VECTORS
C      IFLAG -- INDICATES STATUS OF INTEGRATION
C      WORK(*),IWORK(*) -- ARRAYS TO HOLD INFORMATION INTERNAL TO CODE
C           WHICH IS NECESSARY FOR SUBSEQUENT CALLS
C
C   FIRST CALL TO ODE --
C
C   THE USER MUST PROVIDE STORAGE IN HIS CALLING PROGRAM FOR THE ARRAYS
C   IN THE CALL LIST,
C      Y(NEQN), WORK(100+21*NEQN), IWORK(5),
C   DECLARE  F  IN AN EXTERNAL STATEMENT, SUPPLY THE SUBROUTINE
C   F(T,Y,YP) TO EVALUATE
C      DY(I)/DT = YP(I) = F(T,Y(1),Y(2),...,Y(NEQN))
C   AND INITIALIZE THE PARAMETERS0
C      NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED
C      Y(*) -- VECTOR OF INITIAL CONDITIONS
C      T -- STARTING POINT OF INTEGRATION
C      TOUT -- POINT AT WHICH SOLUTION IS DESIRED
C      RELERR,ABSERR -- RELATIVE AND ABSOLUTE LOCAL ERROR TOLERANCES
C      IFLAG -- +1,-1.  INDICATOR TO INITIALIZE THE CODE.  NORMAL INPUT
C           IS +1.  THE USER SHOULD SET IFLAG=-1 ONLY IF IT IS
C           IMPOSSIBLE TO CONTINUE THE INTEGRATION BEYOND  TOUT .
C   ALL PARAMETERS EXCEPT  F ,  NEQN  AND  TOUT  MAY BE ALTERED BY THE
C   CODE ON OUTPUT SO MUST BE VARIABLES IN THE CALLING PROGRAM.
C
C   OUTPUT FROM  ODE  --
C
C      NEQN -- UNCHANGED
C      Y(*) -- SOLUTION AT  T
C      T -- LAST POINT REACHED IN INTEGRATION.  NORMAL RETURN HAS
C           T = TOUT .
C      TOUT -- UNCHANGED
C      RELERR,ABSERR -- NORMAL RETURN HAS TOLERANCES UNCHANGED.  IFLAG=3
C           SIGNALS TOLERANCES INCREASED
C      IFLAG = 2 -- NORMAL RETURN.  INTEGRATION REACHED  TOUT
C            = 3 -- INTEGRATION DID NOT REACH  TOUT  BECAUSE ERROR
C                   TOLERANCES TOO SMALL.  RELERR ,  ABSERR  INCREASED
C                   APPROPRIATELY FOR CONTINUING
C            = 4 -- INTEGRATION DID NOT REACH  TOUT  BECAUSE MORE THAN
C                   MAXNUM STEPS NEEDED
C            = 5 -- INTEGRATION DID NOT REACH  TOUT  BECAUSE EQUATIONS
C                   APPEAR TO BE STIFF
C            = 6 -- INTEGRATION DID NOT REACH  TOUT  BECAUSE SOLUTION
C                   VANISHED MAKING PURE RELATIVE ERROR IMPOSSIBLE.
C                   MUST USE NON-ZERO  ABSERR  TO CONTINUE.
C            = 7 -- INVALID INPUT PARAMETERS (FATAL ERROR)
C           THE VALUE OF  IFLAG  IS RETURNED NEGATIVE WHEN THE INPUT
C           VALUE IS NEGATIVE AND THE INTEGRATION DOES NOT REACH  TOUT ,
C           I.E., -3, -4, -5, -6.
C      WORK(*),IWORK(*) -- INFORMATION GENERALLY OF NO INTEREST TO THE
C           USER BUT NECESSARY FOR SUBSEQUENT CALLS.
C
C   SUBSEQUENT CALLS TO  ODE --
C
C   SUBROUTINE  ODE  RETURNS WITH ALL INFORMATION NEEDED TO CONTINUE
C   THE INTEGRATION.  IF THE INTEGRATION REACHED  TOUT , THE USER NEED
C   ONLY DEFINE A NEW  TOUT  AND CALL AGAIN.  IF THE INTEGRATION DID NOT
C   REACH  TOUT  AND THE USER WANTS TO CONTINUE, HE JUST CALLS AGAIN.
C   IN THE CASE  IFLAG=6 , THE USER MUST ALSO ALTER THE ERROR CRITERION.
C   THE OUTPUT VALUE OF  IFLAG  IS THE APPROPRIATE INPUT VALUE FOR
C   SUBSEQUENT CALLS.  THE ONLY SITUATION IN WHICH IT SHOULD BE ALTERED
C   IS TO STOP THE INTEGRATION INTERNALLY AT THE NEW  TOUT , I.E.,
C   CHANGE OUTPUT  IFLAG=2  TO INPUT  IFLAG=-2 .  ERROR TOLERANCES MAY
C   BE CHANGED BY THE USER BEFORE CONTINUING.  ALL OTHER PARAMETERS MUST
C   REMAIN UNCHANGED.
C
      LOGICAL START,PHASE1,NORND
      DIMENSION Y(NEQN),WORK(1),IWORK(5)
      EXTERNAL F
      DATA IALPHA,IBETA,ISIG,IV,IW,IG,IPHASE,IPSI,IX,IH,IHOLD,ISTART,
     1  ITOLD,IDELSN/1,13,25,38,50,62,75,76,88,89,90,91,92,93/
      IYY = 100
      IWT = IYY + NEQN
      IP = IWT + NEQN
      IYP = IP + NEQN
      IYPOUT = IYP + NEQN
      IPHI = IYPOUT + NEQN
      IF(IABS(IFLAG) .LT. 2  .OR.  IABS(IFLAG) .GT. 6) GO TO 1
      START = WORK(ISTART) .GT. 0.0
      PHASE1 = WORK(IPHASE) .GT. 0.0
      NORND = IWORK(2) .NE. -1
 1    CALL DE1(F,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,WORK(IYY),
     1  WORK(IWT),WORK(IP),WORK(IYP),WORK(IYPOUT),WORK(IPHI),
     2  WORK(IALPHA),WORK(IBETA),WORK(ISIG),WORK(IV),WORK(IW),WORK(IG),
     3  PHASE1,WORK(IPSI),WORK(IX),WORK(IH),WORK(IHOLD),START,
     4  WORK(ITOLD),WORK(IDELSN),IWORK(1),NORND,IWORK(3),IWORK(4),
     5  IWORK(5))
      WORK(ISTART) = -1.0
      IF(START) WORK(ISTART) = 1.0
      WORK(IPHASE) = -1.0
      IF(PHASE1) WORK(IPHASE) = 1.0
      IWORK(2) = -1
      IF(NORND) IWORK(2) = 1
      RETURN
      END
