changeset 4000:1760b2ce8ef6

[project @ 2002-07-24 18:26:59 by jwe]
author jwe
date Wed, 24 Jul 2002 18:27:00 +0000
parents b4b4515af951
children c3b8a9ff236b
files libcruft/dasrt/Makefile.in libcruft/dasrt/ddasrt.f libcruft/dasrt/drchek.f libcruft/dasrt/droots.f libcruft/dasrt/xerrwv.f
diffstat 5 files changed, 2035 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
new file mode 100644
--- /dev/null
+++ b/libcruft/dasrt/Makefile.in
@@ -0,0 +1,19 @@
+#
+# Makefile for octave's libcruft/dasrt directory
+#
+# John W. Eaton
+# jwe@bevo.che.wisc.edu
+# University of Wisconsin-Madison
+# Department of Chemical Engineering
+
+TOPDIR = ../..
+
+srcdir = @srcdir@
+top_srcdir = @top_srcdir@
+VPATH = @srcdir@
+
+EXTERNAL_DISTFILES = $(DISTFILES)
+
+include $(TOPDIR)/Makeconf
+
+include ../Makerules
new file mode 100644
--- /dev/null
+++ b/libcruft/dasrt/ddasrt.f
@@ -0,0 +1,1559 @@
+      SUBROUTINE DDASRT (RES,NEQ,T,Y,YPRIME,TOUT,
+     *  INFO,RTOL,ATOL,IDID,RWORK,LRW,IWORK,LIW,RPAR,IPAR,JAC,
+     *  G,NG,JROOT)
+C
+C***BEGIN PROLOGUE  DDASRT
+C***DATE WRITTEN   821001   (YYMMDD)
+C***REVISION DATE  910624   (YYMMDD)
+C***KEYWORDS  DIFFERENTIAL/ALGEBRAIC,BACKWARD DIFFERENTIATION FORMULAS
+C             IMPLICIT DIFFERENTIAL SYSTEMS
+C***AUTHOR  PETZOLD,LINDA R.,COMPUTING AND MATHEMATICS RESEARCH DIVISION
+C             LAWRENCE LIVERMORE NATIONAL LABORATORY
+C             L - 316, P.O. Box 808,
+C             LIVERMORE, CA.    94550
+C***PURPOSE  This code solves a system of differential/algebraic
+C            equations of the form F(T,Y,YPRIME) = 0.
+C***DESCRIPTION
+C
+C *Usage:
+C
+C      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+C      EXTERNAL RES, JAC, G
+C      INTEGER NEQ, INFO(N), IDID, LRW, LIW, IWORK(LIW), IPAR, NG,
+C     *   JROOT(NG)
+C      DOUBLE PRECISION T, Y(NEQ), YPRIME(NEQ), TOUT, RTOL, ATOL,
+C     *   RWORK(LRW), RPAR
+C
+C      CALL DDASRT (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
+C     *   IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC)
+C
+C
+C
+C *Arguments:
+C
+C  RES:EXT  This is a subroutine which you provide to define the
+C           differential/algebraic system.
+C
+C  NEQ:IN  This is the number of equations to be solved.
+C
+C  T:INOUT  This is the current value of the independent variable.
+C
+C  Y(*):INOUT  This array contains the solution components at T.
+C
+C  YPRIME(*):INOUT  This array contains the derivatives of the solution
+C                   components at T.
+C
+C  TOUT:IN  This is a point at which a solution is desired.
+C
+C  INFO(N):IN  The basic task of the code is to solve the system from T
+C              to TOUT and return an answer at TOUT.  INFO is an integer
+C              array which is used to communicate exactly how you want
+C              this task to be carried out.  N must be greater than or
+C              equal to 15.
+C
+C  RTOL,ATOL:INOUT  These quantities represent absolute and relative
+C                   error tolerances which you provide to indicate how
+C                   accurately you wish the solution to be computed.
+C                   You may choose them to be both scalars or else
+C                   both vectors.
+C
+C  IDID:OUT  This scalar quantity is an indicator reporting what the
+C            code did.  You must monitor this integer variable to decide
+C            what action to take next.
+C
+C  RWORK:WORK  A real work array of length LRW which provides the
+C               code with needed storage space.
+C
+C  LRW:IN  The length of RWORK.
+C
+C  IWORK:WORK  An integer work array of length LIW which probides the
+C               code with needed storage space.
+C
+C  LIW:IN  The length of IWORK.
+C
+C  RPAR,IPAR:IN  These are real and integer parameter arrays which
+C                you can use for communication between your calling
+C                program and the RES subroutine (and the JAC subroutine)
+C
+C  JAC:EXT  This is the name of a subroutine which you may choose to
+C           provide for defining a matrix of partial derivatives
+C           described below.
+C
+C  G  This is the name of the subroutine for defining
+C     constraint functions, G(T,Y), whose roots are desired
+C     during the integration.  This name must be declared
+C     external in the calling program.
+C
+C  NG  This is the number of constraint functions G(I).
+C      If there are none, set NG=0, and pass a dummy name
+C      for G.
+C
+C  JROOT  This is an integer array of length NG for output
+C         of root information.
+C
+C
+C *Description
+C
+C  QUANTITIES WHICH MAY BE ALTERED BY THE CODE ARE
+C     T,Y(*),YPRIME(*),INFO(1),RTOL,ATOL,
+C     IDID,RWORK(*) AND IWORK(*).
+C
+C  Subroutine DDASRT uses the backward differentiation formulas of
+C  orders one through five to solve a system of the above form for Y and
+C  YPRIME.  Values for Y and YPRIME at the initial time must be given as
+C  input.  These values must be consistent, (that is, if T,Y,YPRIME are
+C  the given initial values, they must satisfy F(T,Y,YPRIME) = 0.).  The
+C  subroutine solves the system from T to TOUT.
+C  It is easy to continue the solution to get results at additional
+C  TOUT.  This is the interval mode of operation.  Intermediate results
+C  can also be obtained easily by using the intermediate-output
+C  capability.  If DDASRT detects a sign-change in G(T,Y), then
+C  it will return the intermediate value of T and Y for which
+C  G(T,Y) = 0.
+C
+C  ---------INPUT-WHAT TO DO ON THE FIRST CALL TO DDASRT---------------
+C
+C
+C  The first call of the code is defined to be the start of each new
+C  problem. Read through the descriptions of all the following items,
+C  provide sufficient storage space for designated arrays, set
+C  appropriate variables for the initialization of the problem, and
+C  give information about how you want the problem to be solved.
+C
+C
+C  RES -- Provide a subroutine of the form
+C             SUBROUTINE RES(T,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
+C         to define the system of differential/algebraic
+C         equations which is to be solved. For the given values
+C         of T,Y and YPRIME, the subroutine should
+C         return the residual of the defferential/algebraic
+C         system
+C             DELTA = F(T,Y,YPRIME)
+C         (DELTA(*) is a vector of length NEQ which is
+C         output for RES.)
+C
+C         Subroutine RES must not alter T,Y or YPRIME.
+C         You must declare the name RES in an external
+C         statement in your program that calls DDASRT.
+C         You must dimension Y,YPRIME and DELTA in RES.
+C
+C         IRES is an integer flag which is always equal to
+C         zero on input. Subroutine RES should alter IRES
+C         only if it encounters an illegal value of Y or
+C         a stop condition. Set IRES = -1 if an input value
+C         is illegal, and DDASRT will try to solve the problem
+C         without getting IRES = -1. If IRES = -2, DDASRT
+C         will return control to the calling program
+C         with IDID = -11.
+C
+C         RPAR and IPAR are real and integer parameter arrays which
+C         you can use for communication between your calling program
+C         and subroutine RES. They are not altered by DDASRT. If you
+C         do not need RPAR or IPAR, ignore these parameters by treat-
+C         ing them as dummy arguments. If you do choose to use them,
+C         dimension them in your calling program and in RES as arrays
+C         of appropriate length.
+C
+C  NEQ -- Set it to the number of differential equations.
+C         (NEQ .GE. 1)
+C
+C  T -- Set it to the initial point of the integration.
+C       T must be defined as a variable.
+C
+C  Y(*) -- Set this vector to the initial values of the NEQ solution
+C          components at the initial point. You must dimension Y of
+C          length at least NEQ in your calling program.
+C
+C  YPRIME(*) -- Set this vector to the initial values of
+C               the NEQ first derivatives of the solution
+C               components at the initial point. You
+C               must dimension YPRIME at least NEQ
+C               in your calling program. If you do not
+C               know initial values of some of the solution
+C               components, see the explanation of INFO(11).
+C
+C  TOUT - Set it to the first point at which a solution
+C         is desired. You can not take TOUT = T.
+C         integration either forward in T (TOUT .GT. T) or
+C         backward in T (TOUT .LT. T) is permitted.
+C
+C         The code advances the solution from T to TOUT using
+C         step sizes which are automatically selected so as to
+C         achieve the desired accuracy. If you wish, the code will
+C         return with the solution and its derivative at
+C         intermediate steps (intermediate-output mode) so that
+C         you can monitor them, but you still must provide TOUT in
+C         accord with the basic aim of the code.
+C
+C         the first step taken by the code is a critical one
+C         because it must reflect how fast the solution changes near
+C         the initial point. The code automatically selects an
+C         initial step size which is practically always suitable for
+C         the problem. By using the fact that the code will not step
+C         past TOUT in the first step, you could, if necessary,
+C         restrict the length of the initial step size.
+C
+C         For some problems it may not be permissable to integrate
+C         past a point TSTOP because a discontinuity occurs there
+C         or the solution or its derivative is not defined beyond
+C         TSTOP. When you have declared a TSTOP point (SEE INFO(4)
+C         and RWORK(1)), you have told the code not to integrate
+C         past TSTOP. In this case any TOUT beyond TSTOP is invalid
+C         input.
+C
+C  INFO(*) - Use the INFO array to give the code more details about
+C            how you want your problem solved. This array should be
+C            dimensioned of length 15, though DDASRT uses
+C            only the first twelve entries. You must respond to all of
+C            the following items which are arranged as questions. The
+C            simplest use of the code corresponds to answering all
+C            questions as yes, i.e. setting all entries of INFO to 0.
+C
+C       INFO(1) - This parameter enables the code to initialize
+C              itself. You must set it to indicate the start of every
+C              new problem.
+C
+C          **** Is this the first call for this problem ...
+C                Yes - Set INFO(1) = 0
+C                 No - Not applicable here.
+C                      See below for continuation calls.  ****
+C
+C       INFO(2) - How much accuracy you want of your solution
+C              is specified by the error tolerances RTOL and ATOL.
+C              The simplest use is to take them both to be scalars.
+C              To obtain more flexibility, they can both be vectors.
+C              The code must be told your choice.
+C
+C          **** Are both error tolerances RTOL, ATOL scalars ...
+C                Yes - Set INFO(2) = 0
+C                      and input scalars for both RTOL and ATOL
+C                 No - Set INFO(2) = 1
+C                      and input arrays for both RTOL and ATOL ****
+C
+C       INFO(3) - The code integrates from T in the direction
+C              of TOUT by steps. If you wish, it will return the
+C              computed solution and derivative at the next
+C              intermediate step (the intermediate-output mode) or
+C              TOUT, whichever comes first. This is a good way to
+C              proceed if you want to see the behavior of the solution.
+C              If you must have solutions at a great many specific
+C              TOUT points, this code will compute them efficiently.
+C
+C          **** Do you want the solution only at
+C                TOUT (and not at the next intermediate step) ...
+C                 Yes - Set INFO(3) = 0
+C                  No - Set INFO(3) = 1 ****
+C
+C       INFO(4) - To handle solutions at a great many specific
+C              values TOUT efficiently, this code may integrate past
+C              TOUT and interpolate to obtain the result at TOUT.
+C              Sometimes it is not possible to integrate beyond some
+C              point TSTOP because the equation changes there or it is
+C              not defined past TSTOP. Then you must tell the code
+C              not to go past.
+C
+C           **** Can the integration be carried out without any
+C                restrictions on the independent variable T ...
+C                 Yes - Set INFO(4)=0
+C                  No - Set INFO(4)=1
+C                       and define the stopping point TSTOP by
+C                       setting RWORK(1)=TSTOP ****
+C
+C       INFO(5) - To solve differential/algebraic problems it is
+C              necessary to use a matrix of partial derivatives of the
+C              system of differential equations. If you do not
+C              provide a subroutine to evaluate it analytically (see
+C              description of the item JAC in the call list), it will
+C              be approximated by numerical differencing in this code.
+C              although it is less trouble for you to have the code
+C              compute partial derivatives by numerical differencing,
+C              the solution will be more reliable if you provide the
+C              derivatives via JAC. Sometimes numerical differencing
+C              is cheaper than evaluating derivatives in JAC and
+C              sometimes it is not - this depends on your problem.
+C
+C           **** Do you want the code to evaluate the partial
+C                derivatives automatically by numerical differences ...
+C                   Yes - Set INFO(5)=0
+C                    No - Set INFO(5)=1
+C                  and provide subroutine JAC for evaluating the
+C                  matrix of partial derivatives ****
+C
+C       INFO(6) - DDASRT will perform much better if the matrix of
+C              partial derivatives, DG/DY + CJ*DG/DYPRIME,
+C              (here CJ is a scalar determined by DDASRT)
+C              is banded and the code is told this. In this
+C              case, the storage needed will be greatly reduced,
+C              numerical differencing will be performed much cheaper,
+C              and a number of important algorithms will execute much
+C              faster. The differential equation is said to have
+C              half-bandwidths ML (lower) and MU (upper) if equation i
+C              involves only unknowns Y(J) with
+C                             I-ML .LE. J .LE. I+MU
+C              for all I=1,2,...,NEQ. Thus, ML and MU are the widths
+C              of the lower and upper parts of the band, respectively,
+C              with the main diagonal being excluded. If you do not
+C              indicate that the equation has a banded matrix of partial
+C              derivatives, the code works with a full matrix of NEQ**2
+C              elements (stored in the conventional way). Computations
+C              with banded matrices cost less time and storage than with
+C              full matrices if 2*ML+MU .LT. NEQ. If you tell the
+C              code that the matrix of partial derivatives has a banded
+C              structure and you want to provide subroutine JAC to
+C              compute the partial derivatives, then you must be careful
+C              to store the elements of the matrix in the special form
+C              indicated in the description of JAC.
+C
+C          **** Do you want to solve the problem using a full
+C               (dense) matrix (and not a special banded
+C               structure) ...
+C                Yes - Set INFO(6)=0
+C                 No - Set INFO(6)=1
+C                       and provide the lower (ML) and upper (MU)
+C                       bandwidths by setting
+C                       IWORK(1)=ML
+C                       IWORK(2)=MU ****
+C
+C
+C        INFO(7) -- You can specify a maximum (absolute value of)
+C              stepsize, so that the code
+C              will avoid passing over very
+C              large regions.
+C
+C          ****  Do you want the code to decide
+C                on its own maximum stepsize?
+C                Yes - Set INFO(7)=0
+C                 No - Set INFO(7)=1
+C                      and define HMAX by setting
+C                      RWORK(2)=HMAX ****
+C
+C        INFO(8) -- Differential/algebraic problems
+C              may occaisionally suffer from
+C              severe scaling difficulties on the
+C              first step. If you know a great deal
+C              about the scaling of your problem, you can
+C              help to alleviate this problem by
+C              specifying an initial stepsize H0.
+C
+C          ****  Do you want the code to define
+C                its own initial stepsize?
+C                Yes - Set INFO(8)=0
+C                 No - Set INFO(8)=1
+C                      and define H0 by setting
+C                      RWORK(3)=H0 ****
+C
+C        INFO(9) -- If storage is a severe problem,
+C              you can save some locations by
+C              restricting the maximum order MAXORD.
+C              the default value is 5. for each
+C              order decrease below 5, the code
+C              requires NEQ fewer locations, however
+C              it is likely to be slower. In any
+C              case, you must have 1 .LE. MAXORD .LE. 5
+C          ****  Do you want the maximum order to
+C                default to 5?
+C                Yes - Set INFO(9)=0
+C                 No - Set INFO(9)=1
+C                      and define MAXORD by setting
+C                      IWORK(3)=MAXORD ****
+C
+C        INFO(10) --If you know that the solutions to your equations
+C               will always be nonnegative, it may help to set this
+C               parameter. However, it is probably best to
+C               try the code without using this option first,
+C               and only to use this option if that doesn't
+C               work very well.
+C           ****  Do you want the code to solve the problem without
+C                 invoking any special nonnegativity constraints?
+C                  Yes - Set INFO(10)=0
+C                   No - Set INFO(10)=1
+C
+C        INFO(11) --DDASRT normally requires the initial T,
+C               Y, and YPRIME to be consistent. That is,
+C               you must have F(T,Y,YPRIME) = 0 at the initial
+C               time. If you do not know the initial
+C               derivative precisely, you can let DDASRT try
+C               to compute it.
+C          ****   Are the initial T, Y, YPRIME consistent?
+C                 Yes - Set INFO(11) = 0
+C                  No - Set INFO(11) = 1,
+C                       and set YPRIME to an initial approximation
+C                       to YPRIME.  (If you have no idea what
+C                       YPRIME should be, set it to zero. Note
+C                       that the initial Y should be such
+C                       that there must exist a YPRIME so that
+C                       F(T,Y,YPRIME) = 0.)
+C
+C        INFO(12) --Maximum number of steps.
+C          ****   Do you want to let DDASRT use the default limit for
+C                 the number of steps?
+C                 Yes - Set INFO(11) = 0
+C                  No - Set INFO(11) = 1,
+C                       and define the maximum number of steps
+C                       by setting IWORK(19)=MXSTEP
+C
+C   RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOL
+C               error tolerances to tell the code how accurately you
+C               want the solution to be computed. They must be defined
+C               as variables because the code may change them. You
+C               have two choices --
+C                     Both RTOL and ATOL are scalars. (INFO(2)=0)
+C                     Both RTOL and ATOL are vectors. (INFO(2)=1)
+C               in either case all components must be non-negative.
+C
+C               The tolerances are used by the code in a local error
+C               test at each step which requires roughly that
+C                     ABS(LOCAL ERROR) .LE. RTOL*ABS(Y)+ATOL
+C               for each vector component.
+C               (More specifically, a root-mean-square norm is used to
+C               measure the size of vectors, and the error test uses the
+C               magnitude of the solution at the beginning of the step.)
+C
+C               The true (global) error is the difference between the
+C               true solution of the initial value problem and the
+C               computed approximation. Practically all present day
+C               codes, including this one, control the local error at
+C               each step and do not even attempt to control the global
+C               error directly.
+C               Usually, but not always, the true accuracy of the
+C               computed Y is comparable to the error tolerances. This
+C               code will usually, but not always, deliver a more
+C               accurate solution if you reduce the tolerances and
+C               integrate again. By comparing two such solutions you
+C               can get a fairly reliable idea of the true error in the
+C               solution at the bigger tolerances.
+C
+C               Setting ATOL=0. results in a pure relative error test on
+C               that component. Setting RTOL=0. results in a pure
+C               absolute error test on that component. A mixed test
+C               with non-zero RTOL and ATOL corresponds roughly to a
+C               relative error test when the solution component is much
+C               bigger than ATOL and to an absolute error test when the
+C               solution component is smaller than the threshhold ATOL.
+C
+C               The code will not attempt to compute a solution at an
+C               accuracy unreasonable for the machine being used. It
+C               will advise you if you ask for too much accuracy and
+C               inform you as to the maximum accuracy it believes
+C               possible.
+C
+C  RWORK(*) --  Dimension this real work array of length LRW in your
+C               calling program.
+C
+C  LRW -- Set it to the declared length of the RWORK array.
+C               You must have
+C                    LRW .GE. 50+(MAXORD+4)*NEQ+NEQ**2
+C               for the full (dense) JACOBIAN case (when INFO(6)=0), or
+C                    LRW .GE. 50+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQ
+C               for the banded user-defined JACOBIAN case
+C               (when INFO(5)=1 and INFO(6)=1), or
+C                     LRW .GE. 50+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQ
+C                           +2*(NEQ/(ML+MU+1)+1)
+C               for the banded finite-difference-generated JACOBIAN case
+C               (when INFO(5)=0 and INFO(6)=1)
+C
+C  IWORK(*) --  Dimension this integer work array of length LIW in
+C               your calling program.
+C
+C  LIW -- Set it to the declared length of the IWORK array.
+C               you must have LIW .GE. 20+NEQ
+C
+C  RPAR, IPAR -- These are parameter arrays, of real and integer
+C               type, respectively. You can use them for communication
+C               between your program that calls DDASRT and the
+C               RES subroutine (and the JAC subroutine). They are not
+C               altered by DDASRT. If you do not need RPAR or IPAR,
+C               ignore these parameters by treating them as dummy
+C               arguments. If you do choose to use them, dimension
+C               them in your calling program and in RES (and in JAC)
+C               as arrays of appropriate length.
+C
+C  JAC -- If you have set INFO(5)=0, you can ignore this parameter
+C               by treating it as a dummy argument. Otherwise, you must
+C               provide a subroutine of the form
+C               JAC(T,Y,YPRIME,PD,CJ,RPAR,IPAR)
+C               to define the matrix of partial derivatives
+C               PD=DG/DY+CJ*DG/DYPRIME
+C               CJ is a scalar which is input to JAC.
+C               For the given values of T,Y,YPRIME, the
+C               subroutine must evaluate the non-zero partial
+C               derivatives for each equation and each solution
+C               component, and store these values in the
+C               matrix PD. The elements of PD are set to zero
+C               before each call to JAC so only non-zero elements
+C               need to be defined.
+C
+C               Subroutine JAC must not alter T,Y,(*),YPRIME(*), or CJ.
+C               You must declare the name JAC in an
+C               EXTERNAL STATEMENT in your program that calls
+C               DDASRT. You must dimension Y, YPRIME and PD
+C               in JAC.
+C
+C               The way you must store the elements into the PD matrix
+C               depends on the structure of the matrix which you
+C               indicated by INFO(6).
+C               *** INFO(6)=0 -- Full (dense) matrix ***
+C                   Give PD a first dimension of NEQ.
+C                   When you evaluate the (non-zero) partial derivative
+C                   of equation I with respect to variable J, you must
+C                   store it in PD according to
+C                   PD(I,J) = * DF(I)/DY(J)+CJ*DF(I)/DYPRIME(J)*
+C               *** INFO(6)=1 -- Banded JACOBIAN with ML lower and MU
+C                   upper diagonal bands (refer to INFO(6) description
+C                   of ML and MU) ***
+C                   Give PD a first dimension of 2*ML+MU+1.
+C                   when you evaluate the (non-zero) partial derivative
+C                   of equation I with respect to variable J, you must
+C                   store it in PD according to
+C                   IROW = I - J + ML + MU + 1
+C                   PD(IROW,J) = *DF(I)/DY(J)+CJ*DF(I)/DYPRIME(J)*
+C               RPAR and IPAR are real and integer parameter arrays
+C               which you can use for communication between your calling
+C               program and your JACOBIAN subroutine JAC. They are not
+C               altered by DDASRT. If you do not need RPAR or IPAR,
+C               ignore these parameters by treating them as dummy
+C               arguments. If you do choose to use them, dimension
+C               them in your calling program and in JAC as arrays of
+C               appropriate length.
+C
+C  G -- This is the name of the subroutine for defining constraint
+C               functions, whose roots are desired during the
+C               integration.  It is to have the form
+C                   SUBROUTINE G(NEQ,T,Y,NG,GOUT,RPAR,IPAR)
+C                   DIMENSION Y(NEQ),GOUT(NG),
+C               where NEQ, T, Y and NG are INPUT, and the array GOUT is
+C               output.  NEQ, T, and Y have the same meaning as in the
+C               RES routine, and GOUT is an array of length NG.
+C               For I=1,...,NG, this routine is to load into GOUT(I)
+C               the value at (T,Y) of the I-th constraint function G(I).
+C               DDASRT will find roots of the G(I) of odd multiplicity
+C               (that is, sign changes) as they occur during
+C               the integration.  G must be declared EXTERNAL in the
+C               calling program.
+C
+C               CAUTION..because of numerical errors in the functions
+C               G(I) due to roundoff and integration error, DDASRT
+C               may return false roots, or return the same root at two
+C               or more nearly equal values of T.  If such false roots
+C               are suspected, the user should consider smaller error
+C               tolerances and/or higher precision in the evaluation of
+C               the G(I).
+C
+C               If a root of some G(I) defines the end of the problem,
+C               the input to DDASRT should nevertheless allow
+C               integration to a point slightly past that ROOT, so
+C               that DDASRT can locate the root by interpolation.
+C
+C  NG -- The number of constraint functions G(I).  If there are none,
+C               set NG = 0, and pass a dummy name for G.
+C
+C JROOT -- This is an integer array of length NG.  It is used only for
+C               output.  On a return where one or more roots have been
+C               found, JROOT(I)=1 If G(I) has a root at T,
+C               or JROOT(I)=0 if not.
+C
+C
+C
+C  OPTIONALLY REPLACEABLE NORM ROUTINE:
+C  DDASRT uses a weighted norm DDANRM to measure the size
+C  of vectors such as the estimated error in each step.
+C  A FUNCTION subprogram
+C    DOUBLE PRECISION FUNCTION DDANRM(NEQ,V,WT,RPAR,IPAR)
+C    DIMENSION V(NEQ),WT(NEQ)
+C  is used to define this norm. Here, V is the vector
+C  whose norm is to be computed, and WT is a vector of
+C  weights.  A DDANRM routine has been included with DDASRT
+C  which computes the weighted root-mean-square norm
+C  given by
+C    DDANRM=SQRT((1/NEQ)*SUM(V(I)/WT(I))**2)
+C  this norm is suitable for most problems. In some
+C  special cases, it may be more convenient and/or
+C  efficient to define your own norm by writing a function
+C  subprogram to be called instead of DDANRM. This should
+C  ,however, be attempted only after careful thought and
+C  consideration.
+C
+C
+C------OUTPUT-AFTER ANY RETURN FROM DDASRT----
+C
+C  The principal aim of the code is to return a computed solution at
+C  TOUT, although it is also possible to obtain intermediate results
+C  along the way. To find out whether the code achieved its goal
+C  or if the integration process was interrupted before the task was
+C  completed, you must check the IDID parameter.
+C
+C
+C   T -- The solution was successfully advanced to the
+C               output value of T.
+C
+C   Y(*) -- Contains the computed solution approximation at T.
+C
+C   YPRIME(*) -- Contains the computed derivative
+C               approximation at T.
+C
+C   IDID -- Reports what the code did.
+C
+C                     *** Task completed ***
+C                Reported by positive values of IDID
+C
+C           IDID = 1 -- A step was successfully taken in the
+C                   intermediate-output mode. The code has not
+C                   yet reached TOUT.
+C
+C           IDID = 2 -- The integration to TSTOP was successfully
+C                   completed (T=TSTOP) by stepping exactly to TSTOP.
+C
+C           IDID = 3 -- The integration to TOUT was successfully
+C                   completed (T=TOUT) by stepping past TOUT.
+C                   Y(*) is obtained by interpolation.
+C                   YPRIME(*) is obtained by interpolation.
+C
+C           IDID = 4 -- The integration was successfully completed
+C                   by finding one or more roots of G at T.
+C
+C                    *** Task interrupted ***
+C                Reported by negative values of IDID
+C
+C           IDID = -1 -- A large amount of work has been expended.
+C                   (About INFO(12) steps)
+C
+C           IDID = -2 -- The error tolerances are too stringent.
+C
+C           IDID = -3 -- The local error test cannot be satisfied
+C                   because you specified a zero component in ATOL
+C                   and the corresponding computed solution
+C                   component is zero. Thus, a pure relative error
+C                   test is impossible for this component.
+C
+C           IDID = -6 -- DDASRT had repeated error test
+C                   failures on the last attempted step.
+C
+C           IDID = -7 -- The corrector could not converge.
+C
+C           IDID = -8 -- The matrix of partial derivatives
+C                   is singular.
+C
+C           IDID = -9 -- The corrector could not converge.
+C                   there were repeated error test failures
+C                   in this step.
+C
+C           IDID =-10 -- The corrector could not converge
+C                   because IRES was equal to minus one.
+C
+C           IDID =-11 -- IRES equal to -2 was encountered
+C                   and control is being returned to the
+C                   calling program.
+C
+C           IDID =-12 -- DDASRT failed to compute the initial
+C                   YPRIME.
+C
+C
+C
+C           IDID = -13,..,-32 -- Not applicable for this code
+C
+C                    *** Task terminated ***
+C                Reported by the value of IDID=-33
+C
+C           IDID = -33 -- The code has encountered trouble from which
+C                   it cannot recover. A message is printed
+C                   explaining the trouble and control is returned
+C                   to the calling program. For example, this occurs
+C                   when invalid input is detected.
+C
+C   RTOL, ATOL -- These quantities remain unchanged except when
+C               IDID = -2. In this case, the error tolerances have been
+C               increased by the code to values which are estimated to
+C               be appropriate for continuing the integration. However,
+C               the reported solution at T was obtained using the input
+C               values of RTOL and ATOL.
+C
+C   RWORK, IWORK -- Contain information which is usually of no
+C               interest to the user but necessary for subsequent calls.
+C               However, you may find use for
+C
+C               RWORK(3)--Which contains the step size H to be
+C                       attempted on the next step.
+C
+C               RWORK(4)--Which contains the current value of the
+C                       independent variable, i.e., the farthest point
+C                       integration has reached. This will be different
+C                       from T only when interpolation has been
+C                       performed (IDID=3).
+C
+C               RWORK(7)--Which contains the stepsize used
+C                       on the last successful step.
+C
+C               IWORK(7)--Which contains the order of the method to
+C                       be attempted on the next step.
+C
+C               IWORK(8)--Which contains the order of the method used
+C                       on the last step.
+C
+C               IWORK(11)--Which contains the number of steps taken so
+C                        far.
+C
+C               IWORK(12)--Which contains the number of calls to RES
+C                        so far.
+C
+C               IWORK(13)--Which contains the number of evaluations of
+C                        the matrix of partial derivatives needed so
+C                        far.
+C
+C               IWORK(14)--Which contains the total number
+C                        of error test failures so far.
+C
+C               IWORK(15)--Which contains the total number
+C                        of convergence test failures so far.
+C                        (includes singular iteration matrix
+C                        failures.)
+C
+C               IWORK(16)--Which contains the total number of calls
+C                        to the constraint function g so far
+C
+C
+C
+C   INPUT -- What to do to continue the integration
+C            (calls after the first)                **
+C
+C     This code is organized so that subsequent calls to continue the
+C     integration involve little (if any) additional effort on your
+C     part. You must monitor the IDID parameter in order to determine
+C     what to do next.
+C
+C     Recalling that the principal task of the code is to integrate
+C     from T to TOUT (the interval mode), usually all you will need
+C     to do is specify a new TOUT upon reaching the current TOUT.
+C
+C     Do not alter any quantity not specifically permitted below,
+C     in particular do not alter NEQ,T,Y(*),YPRIME(*),RWORK(*),IWORK(*)
+C     or the differential equation in subroutine RES. Any such
+C     alteration constitutes a new problem and must be treated as such,
+C     i.e., you must start afresh.
+C
+C     You cannot change from vector to scalar error control or vice
+C     versa (INFO(2)), but you can change the size of the entries of
+C     RTOL, ATOL. Increasing a tolerance makes the equation easier
+C     to integrate. Decreasing a tolerance will make the equation
+C     harder to integrate and should generally be avoided.
+C
+C     You can switch from the intermediate-output mode to the
+C     interval mode (INFO(3)) or vice versa at any time.
+C
+C     If it has been necessary to prevent the integration from going
+C     past a point TSTOP (INFO(4), RWORK(1)), keep in mind that the
+C     code will not integrate to any TOUT beyond the currently
+C     specified TSTOP. Once TSTOP has been reached you must change
+C     the value of TSTOP or set INFO(4)=0. You may change INFO(4)
+C     or TSTOP at any time but you must supply the value of TSTOP in
+C     RWORK(1) whenever you set INFO(4)=1.
+C
+C     Do not change INFO(5), INFO(6), IWORK(1), or IWORK(2)
+C     unless you are going to restart the code.
+C
+C                    *** Following a completed task ***
+C     If
+C     IDID = 1, call the code again to continue the integration
+C                  another step in the direction of TOUT.
+C
+C     IDID = 2 or 3, define a new TOUT and call the code again.
+C                  TOUT must be different from T. You cannot change
+C                  the direction of integration without restarting.
+C
+C     IDID = 4, call the code again to continue the integration
+C                  another step in the direction of TOUT.  You may
+C                  change the functions in G after a return with IDID=4,
+C                  but the number of constraint functions NG must remain
+C                  the same.  If you wish to change
+C                  the functions in RES or in G, then you
+C                  must restart the code.
+C
+C                    *** Following an interrupted task ***
+C                  To show the code that you realize the task was
+C                  interrupted and that you want to continue, you
+C                  must take appropriate action and set INFO(1) = 1
+C     If
+C     IDID = -1, The code has reached the step iteration.
+C                  If you want to continue, set INFO(1) = 1 and
+C                  call the code again.  See also INFO(12).
+C
+C     IDID = -2, The error tolerances RTOL, ATOL have been
+C                  increased to values the code estimates appropriate
+C                  for continuing. You may want to change them
+C                  yourself. If you are sure you want to continue
+C                  with relaxed error tolerances, set INFO(1)=1 and
+C                  call the code again.
+C
+C     IDID = -3, A solution component is zero and you set the
+C                  corresponding component of ATOL to zero. If you
+C                  are sure you want to continue, you must first
+C                  alter the error criterion to use positive values
+C                  for those components of ATOL corresponding to zero
+C                  solution components, then set INFO(1)=1 and call
+C                  the code again.
+C
+C     IDID = -4,-5  --- Cannot occur with this code.
+C
+C     IDID = -6, Repeated error test failures occurred on the
+C                  last attempted step in DDASRT. A singularity in the
+C                  solution may be present. If you are absolutely
+C                  certain you want to continue, you should restart
+C                  the integration. (Provide initial values of Y and
+C                  YPRIME which are consistent)
+C
+C     IDID = -7, Repeated convergence test failures occurred
+C                  on the last attempted step in DDASRT. An inaccurate
+C                  or ill-conditioned JACOBIAN may be the problem. If
+C                  you are absolutely certain you want to continue, you
+C                  should restart the integration.
+C
+C     IDID = -8, The matrix of partial derivatives is singular.
+C                  Some of your equations may be redundant.
+C                  DDASRT cannot solve the problem as stated.
+C                  It is possible that the redundant equations
+C                  could be removed, and then DDASRT could
+C                  solve the problem. It is also possible
+C                  that a solution to your problem either
+C                  does not exist or is not unique.
+C
+C     IDID = -9, DDASRT had multiple convergence test
+C                  failures, preceeded by multiple error
+C                  test failures, on the last attempted step.
+C                  It is possible that your problem
+C                  is ill-posed, and cannot be solved
+C                  using this code. Or, there may be a
+C                  discontinuity or a singularity in the
+C                  solution. If you are absolutely certain
+C                  you want to continue, you should restart
+C                  the integration.
+C
+C    IDID =-10, DDASRT had multiple convergence test failures
+C                  because IRES was equal to minus one.
+C                  If you are absolutely certain you want
+C                  to continue, you should restart the
+C                  integration.
+C
+C    IDID =-11, IRES=-2 was encountered, and control is being
+C                  returned to the calling program.
+C
+C    IDID =-12, DDASRT failed to compute the initial YPRIME.
+C               This could happen because the initial
+C               approximation to YPRIME was not very good, or
+C               if a YPRIME consistent with the initial Y
+C               does not exist. The problem could also be caused
+C               by an inaccurate or singular iteration matrix.
+C
+C
+C
+C     IDID = -13,..,-32 --- Cannot occur with this code.
+C
+C                       *** Following a terminated task ***
+C     If IDID= -33, you cannot continue the solution of this
+C                  problem. An attempt to do so will result in your
+C                  run being terminated.
+C
+C  ---------------------------------------------------------------------
+C
+C***REFERENCE
+C      K. E. Brenan, S. L. Campbell, and L. R. Petzold, Numerical 
+C      Solution of Initial-Value Problems in Differential-Algebraic
+C      Equations, Elsevier, New York, 1989.
+C
+C***ROUTINES CALLED  DDASTP,DDAINI,DDANRM,DDAWTS,DDATRP,DRCHEK,DROOTS,
+C                    XERRWV,D1MACH
+C***END PROLOGUE  DDASRT
+C
+C**End
+C
+      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
+      LOGICAL DONE
+      EXTERNAL RES, JAC, G
+      DIMENSION Y(*),YPRIME(*)
+      DIMENSION INFO(15)
+      DIMENSION RWORK(*),IWORK(*)
+      DIMENSION RTOL(*),ATOL(*)
+      DIMENSION RPAR(*),IPAR(*)
+      CHARACTER MSG*80
+C
+C     SET POINTERS INTO IWORK
+      PARAMETER (LML=1, LMU=2, LMXORD=3, LMTYPE=4, LNST=11,
+     *  LNRE=12, LNJE=13, LETF=14, LCTF=15, LNGE=16, LNPD=17,
+     *  LIRFND=18, LMXSTP=19, LIPVT=21, LJCALC=5, LPHASE=6, LK=7,
+     *  LKOLD=8, LNS=9, LNSTL=10, LIWM=1)
+C
+C     SET RELATIVE OFFSET INTO RWORK
+      PARAMETER (NPD=1)
+C
+C     SET POINTERS INTO RWORK
+      PARAMETER (LTSTOP=1, LHMAX=2, LH=3, LTN=4,
+     *  LCJ=5, LCJOLD=6, LHOLD=7, LS=8, LROUND=9,
+     *  LALPHA=11, LBETA=17, LGAMMA=23,
+     *  LPSI=29, LSIGMA=35, LT0=41, LTLAST=42, LALPHR=43, LX2=44,
+     *  LDELTA=51)
+C
+C***FIRST EXECUTABLE STATEMENT  DDASRT
+      IF(INFO(1).NE.0)GO TO 100
+C
+C-----------------------------------------------------------------------
+C     THIS BLOCK IS EXECUTED FOR THE INITIAL CALL ONLY.
+C     IT CONTAINS CHECKING OF INPUTS AND INITIALIZATIONS.
+C-----------------------------------------------------------------------
+C
+C     FIRST CHECK INFO ARRAY TO MAKE SURE ALL ELEMENTS OF INFO
+C     ARE EITHER ZERO OR ONE.
+      DO 10 I=2,12
+         IF(INFO(I).NE.0.AND.INFO(I).NE.1)GO TO 701
+10       CONTINUE
+C
+      IF(NEQ.LE.0)GO TO 702
+C
+C     CHECK AND COMPUTE MAXIMUM ORDER
+      MXORD=5
+      IF(INFO(9).EQ.0)GO TO 20
+         MXORD=IWORK(LMXORD)
+         IF(MXORD.LT.1.OR.MXORD.GT.5)GO TO 703
+20       IWORK(LMXORD)=MXORD
+C
+C     COMPUTE MTYPE,LENPD,LENRW.CHECK ML AND MU.
+      IF(INFO(6).NE.0)GO TO 40
+         LENPD=NEQ**2
+         LENRW=50+(IWORK(LMXORD)+4)*NEQ+LENPD
+         IF(INFO(5).NE.0)GO TO 30
+            IWORK(LMTYPE)=2
+            GO TO 60
+30          IWORK(LMTYPE)=1
+            GO TO 60
+40    IF(IWORK(LML).LT.0.OR.IWORK(LML).GE.NEQ)GO TO 717
+      IF(IWORK(LMU).LT.0.OR.IWORK(LMU).GE.NEQ)GO TO 718
+      LENPD=(2*IWORK(LML)+IWORK(LMU)+1)*NEQ
+      IF(INFO(5).NE.0)GO TO 50
+         IWORK(LMTYPE)=5
+         MBAND=IWORK(LML)+IWORK(LMU)+1
+         MSAVE=(NEQ/MBAND)+1
+         LENRW=50+(IWORK(LMXORD)+4)*NEQ+LENPD+2*MSAVE
+         GO TO 60
+50       IWORK(LMTYPE)=4
+         LENRW=50+(IWORK(LMXORD)+4)*NEQ+LENPD
+C
+C     CHECK LENGTHS OF RWORK AND IWORK
+60    LENIW=20+NEQ
+      IWORK(LNPD)=LENPD
+      IF(LRW.LT.LENRW)GO TO 704
+      IF(LIW.LT.LENIW)GO TO 705
+C
+C     CHECK TO SEE THAT TOUT IS DIFFERENT FROM T
+C     Also check to see that NG is larger than 0.
+      IF(TOUT .EQ. T)GO TO 719
+      IF(NG .LT. 0) GO TO 730
+C
+C     CHECK HMAX
+      IF(INFO(7).EQ.0)GO TO 70
+         HMAX=RWORK(LHMAX)
+         IF(HMAX.LE.0.0D0)GO TO 710
+70    CONTINUE
+C
+C     CHECK AND COMPUTE MAXIMUM STEPS
+      MXSTP=500
+      IF(INFO(12).EQ.0)GO TO 80
+        MXSTP=IWORK(LMXSTP)
+        IF(MXSTP.LT.0)GO TO 703
+80      IWORK(LMXSTP)=MXSTP
+C
+C     INITIALIZE COUNTERS
+      IWORK(LNST)=0
+      IWORK(LNRE)=0
+      IWORK(LNJE)=0
+      IWORK(LNGE)=0
+C
+      IWORK(LNSTL)=0
+      IDID=1
+      GO TO 200
+C
+C-----------------------------------------------------------------------
+C     THIS BLOCK IS FOR CONTINUATION CALLS
+C     ONLY. HERE WE CHECK INFO(1),AND IF THE
+C     LAST STEP WAS INTERRUPTED WE CHECK WHETHER
+C     APPROPRIATE ACTION WAS TAKEN.
+C-----------------------------------------------------------------------
+C
+100   CONTINUE
+      IF(INFO(1).EQ.1)GO TO 110
+      IF(INFO(1).NE.-1)GO TO 701
+C     IF WE ARE HERE, THE LAST STEP WAS INTERRUPTED
+C     BY AN ERROR CONDITION FROM DDASTP,AND
+C     APPROPRIATE ACTION WAS NOT TAKEN. THIS
+C     IS A FATAL ERROR.
+      MSG = 'DASRT--  THE LAST STEP TERMINATED WITH A NEGATIVE'
+      CALL XERRWV(MSG,49,201,0,0,0,0,0,0.0D0,0.0D0)
+      MSG = 'DASRT--  VALUE (=I1) OF IDID AND NO APPROPRIATE'
+      CALL XERRWV(MSG,47,202,0,1,IDID,0,0,0.0D0,0.0D0)
+      MSG = 'DASRT--  ACTION WAS TAKEN. RUN TERMINATED'
+      CALL XERRWV(MSG,41,203,1,0,0,0,0,0.0D0,0.0D0)
+      RETURN
+110   CONTINUE
+      IWORK(LNSTL)=IWORK(LNST)
+C
+C-----------------------------------------------------------------------
+C     THIS BLOCK IS EXECUTED ON ALL CALLS.
+C     THE ERROR TOLERANCE PARAMETERS ARE
+C     CHECKED, AND THE WORK ARRAY POINTERS
+C     ARE SET.
+C-----------------------------------------------------------------------
+C
+200   CONTINUE
+C     CHECK RTOL,ATOL
+      NZFLG=0
+      RTOLI=RTOL(1)
+      ATOLI=ATOL(1)
+      DO 210 I=1,NEQ
+         IF(INFO(2).EQ.1)RTOLI=RTOL(I)
+         IF(INFO(2).EQ.1)ATOLI=ATOL(I)
+         IF(RTOLI.GT.0.0D0.OR.ATOLI.GT.0.0D0)NZFLG=1
+         IF(RTOLI.LT.0.0D0)GO TO 706
+         IF(ATOLI.LT.0.0D0)GO TO 707
+210      CONTINUE
+      IF(NZFLG.EQ.0)GO TO 708
+C
+C     SET UP RWORK STORAGE.IWORK STORAGE IS FIXED
+C     IN DATA STATEMENT.
+      LG0=LDELTA+NEQ
+      LG1=LG0+NG
+      LGX=LG1+NG
+      LE=LGX+NG
+      LWT=LE+NEQ
+      LPHI=LWT+NEQ
+      LPD=LPHI+(IWORK(LMXORD)+1)*NEQ
+      LWM=LPD
+      NTEMP=NPD+IWORK(LNPD)
+      IF(INFO(1).EQ.1)GO TO 400
+C
+C-----------------------------------------------------------------------
+C     THIS BLOCK IS EXECUTED ON THE INITIAL CALL
+C     ONLY. SET THE INITIAL STEP SIZE, AND
+C     THE ERROR WEIGHT VECTOR, AND PHI.
+C     COMPUTE INITIAL YPRIME, IF NECESSARY.
+C-----------------------------------------------------------------------
+C
+300   CONTINUE
+      TN=T
+      IDID=1
+C
+C     SET ERROR WEIGHT VECTOR WT
+      CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR)
+      DO 305 I = 1,NEQ
+         IF(RWORK(LWT+I-1).LE.0.0D0) GO TO 713
+305      CONTINUE
+C
+C     COMPUTE UNIT ROUNDOFF AND HMIN
+      UROUND = D1MACH(4)
+      RWORK(LROUND) = UROUND
+      HMIN = 4.0D0*UROUND*DMAX1(DABS(T),DABS(TOUT))
+C
+C     CHECK INITIAL INTERVAL TO SEE THAT IT IS LONG ENOUGH
+      TDIST = DABS(TOUT - T)
+      IF(TDIST .LT. HMIN) GO TO 714
+C
+C     CHECK H0, IF THIS WAS INPUT
+      IF (INFO(8) .EQ. 0) GO TO 310
+         HO = RWORK(LH)
+         IF ((TOUT - T)*HO .LT. 0.0D0) GO TO 711
+         IF (HO .EQ. 0.0D0) GO TO 712
+         GO TO 320
+310    CONTINUE
+C
+C     COMPUTE INITIAL STEPSIZE, TO BE USED BY EITHER
+C     DDASTP OR DDAINI, DEPENDING ON INFO(11)
+      HO = 0.001D0*TDIST
+      YPNORM = DDANRM(NEQ,YPRIME,RWORK(LWT),RPAR,IPAR)
+      IF (YPNORM .GT. 0.5D0/HO) HO = 0.5D0/YPNORM
+      HO = DSIGN(HO,TOUT-T)
+C     ADJUST HO IF NECESSARY TO MEET HMAX BOUND
+320   IF (INFO(7) .EQ. 0) GO TO 330
+         RH = DABS(HO)/RWORK(LHMAX)
+         IF (RH .GT. 1.0D0) HO = HO/RH
+C     COMPUTE TSTOP, IF APPLICABLE
+330   IF (INFO(4) .EQ. 0) GO TO 340
+         TSTOP = RWORK(LTSTOP)
+         IF ((TSTOP - T)*HO .LT. 0.0D0) GO TO 715
+         IF ((T + HO - TSTOP)*HO .GT. 0.0D0) HO = TSTOP - T
+         IF ((TSTOP - TOUT)*HO .LT. 0.0D0) GO TO 709
+C
+C     COMPUTE INITIAL DERIVATIVE, UPDATING TN AND Y, IF APPLICABLE
+340   IF (INFO(11) .EQ. 0) GO TO 350
+      CALL DDAINI(TN,Y,YPRIME,NEQ,
+     *  RES,JAC,HO,RWORK(LWT),IDID,RPAR,IPAR,
+     *  RWORK(LPHI),RWORK(LDELTA),RWORK(LE),
+     *  RWORK(LWM),IWORK(LIWM),HMIN,RWORK(LROUND),
+     *  INFO(10),NTEMP)
+      IF (IDID .LT. 0) GO TO 390
+C
+C     LOAD H WITH H0.  STORE H IN RWORK(LH)
+350   H = HO
+      RWORK(LH) = H
+C
+C     LOAD Y AND H*YPRIME INTO PHI(*,1) AND PHI(*,2)
+360   ITEMP = LPHI + NEQ
+      DO 370 I = 1,NEQ
+         RWORK(LPHI + I - 1) = Y(I)
+370      RWORK(ITEMP + I - 1) = H*YPRIME(I)
+C
+C     INITIALIZE T0 IN RWORK AND CHECK FOR A ZERO OF G NEAR THE
+C     INITIAL T.
+C
+      RWORK(LT0) = T
+      IWORK(LIRFND) = 0
+      RWORK(LPSI)=H
+      RWORK(LPSI+1)=2.0D0*H
+      IWORK(LKOLD)=1
+      IF(NG .EQ. 0) GO TO 390
+      CALL DRCHEK(1,G,NG,NEQ,T,TOUT,Y,RWORK(LE),RWORK(LPHI),
+     *  RWORK(LPSI),IWORK(LKOLD),RWORK(LG0),RWORK(LG1),
+     *  RWORK(LGX),JROOT,IRT,RWORK(LROUND),INFO(3),
+     *  RWORK,IWORK,RPAR,IPAR)
+      IF(IRT .NE. 0) GO TO 732
+C
+C     Check for a root in the interval (T0,TN], unless DDASRT
+C     did not have to initialize YPRIME.
+C
+      IF(NG .EQ. 0 .OR. INFO(11) .EQ. 0) GO TO 390
+      CALL DRCHEK(3,G,NG,NEQ,TN,TOUT,Y,RWORK(LE),RWORK(LPHI),
+     *  RWORK(LPSI),IWORK(LKOLD),RWORK(LG0),RWORK(LG1),
+     *  RWORK(LGX),JROOT,IRT,RWORK(LROUND),INFO(3),
+     *  RWORK,IWORK,RPAR,IPAR)
+      IF(IRT .NE. 1) GO TO 390
+      IWORK(LIRFND) = 1
+      IDID = 4
+      T = RWORK(LT0)
+      GO TO 580
+C
+390   GO TO 500
+C
+C-------------------------------------------------------
+C     THIS BLOCK IS FOR CONTINUATION CALLS ONLY. ITS
+C     PURPOSE IS TO CHECK STOP CONDITIONS BEFORE
+C     TAKING A STEP.
+C     ADJUST H IF NECESSARY TO MEET HMAX BOUND
+C-------------------------------------------------------
+C
+400   CONTINUE
+      UROUND=RWORK(LROUND)
+      DONE = .FALSE.
+      TN=RWORK(LTN)
+      H=RWORK(LH)
+      IF(NG .EQ. 0) GO TO 405
+C
+C     Check for a zero of G near TN.
+C
+      CALL DRCHEK(2,G,NG,NEQ,TN,TOUT,Y,RWORK(LE),RWORK(LPHI),
+     *  RWORK(LPSI),IWORK(LKOLD),RWORK(LG0),RWORK(LG1),
+     *  RWORK(LGX),JROOT,IRT,RWORK(LROUND),INFO(3),
+     *  RWORK,IWORK,RPAR,IPAR)
+      IF(IRT .NE. 1) GO TO 405
+      IWORK(LIRFND) = 1
+      IDID = 4
+      T = RWORK(LT0)
+      DONE = .TRUE.
+      GO TO 490
+C
+405   CONTINUE
+      IF(INFO(7) .EQ. 0) GO TO 410
+         RH = DABS(H)/RWORK(LHMAX)
+         IF(RH .GT. 1.0D0) H = H/RH
+410   CONTINUE
+      IF(T .EQ. TOUT) GO TO 719
+      IF((T - TOUT)*H .GT. 0.0D0) GO TO 711
+      IF(INFO(4) .EQ. 1) GO TO 430
+      IF(INFO(3) .EQ. 1) GO TO 420
+      IF((TN-TOUT)*H.LT.0.0D0)GO TO 490
+      CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
+     *  RWORK(LPHI),RWORK(LPSI))
+      T=TOUT
+      IDID = 3
+      DONE = .TRUE.
+      GO TO 490
+420   IF((TN-T)*H .LE. 0.0D0) GO TO 490
+      IF((TN - TOUT)*H .GT. 0.0D0) GO TO 425
+      CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD),
+     *  RWORK(LPHI),RWORK(LPSI))
+      T = TN
+      IDID = 1
+      DONE = .TRUE.
+      GO TO 490
+425   CONTINUE
+      CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
+     *  RWORK(LPHI),RWORK(LPSI))
+      T = TOUT
+      IDID = 3
+      DONE = .TRUE.
+      GO TO 490
+430   IF(INFO(3) .EQ. 1) GO TO 440
+      TSTOP=RWORK(LTSTOP)
+      IF((TN-TSTOP)*H.GT.0.0D0) GO TO 715
+      IF((TSTOP-TOUT)*H.LT.0.0D0)GO TO 709
+      IF((TN-TOUT)*H.LT.0.0D0)GO TO 450
+      CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
+     *   RWORK(LPHI),RWORK(LPSI))
+      T=TOUT
+      IDID = 3
+      DONE = .TRUE.
+      GO TO 490
+440   TSTOP = RWORK(LTSTOP)
+      IF((TN-TSTOP)*H .GT. 0.0D0) GO TO 715
+      IF((TSTOP-TOUT)*H .LT. 0.0D0) GO TO 709
+      IF((TN-T)*H .LE. 0.0D0) GO TO 450
+      IF((TN - TOUT)*H .GT. 0.0D0) GO TO 445
+      CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD),
+     *  RWORK(LPHI),RWORK(LPSI))
+      T = TN
+      IDID = 1
+      DONE = .TRUE.
+      GO TO 490
+445   CONTINUE
+      CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
+     *  RWORK(LPHI),RWORK(LPSI))
+      T = TOUT
+      IDID = 3
+      DONE = .TRUE.
+      GO TO 490
+450   CONTINUE
+C     CHECK WHETHER WE ARE WITH IN ROUNDOFF OF TSTOP
+      IF(DABS(TN-TSTOP).GT.100.0D0*UROUND*
+     *   (DABS(TN)+DABS(H)))GO TO 460
+      CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,IWORK(LKOLD),
+     *  RWORK(LPHI),RWORK(LPSI))
+      IDID=2
+      T=TSTOP
+      DONE = .TRUE.
+      GO TO 490
+460   TNEXT=TN+H
+      IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 490
+      H=TSTOP-TN
+      RWORK(LH)=H
+C
+490   IF (DONE) GO TO 590
+C
+C-------------------------------------------------------
+C     THE NEXT BLOCK CONTAINS THE CALL TO THE
+C     ONE-STEP INTEGRATOR DDASTP.
+C     THIS IS A LOOPING POINT FOR THE INTEGRATION STEPS.
+C     CHECK FOR TOO MANY STEPS.
+C     UPDATE WT.
+C     CHECK FOR TOO MUCH ACCURACY REQUESTED.
+C     COMPUTE MINIMUM STEPSIZE.
+C-------------------------------------------------------
+C
+500   CONTINUE
+C     CHECK FOR FAILURE TO COMPUTE INITIAL YPRIME
+      IF (IDID .EQ. -12) GO TO 527
+C
+C     CHECK FOR TOO MANY STEPS
+      IF((IWORK(LNST)-IWORK(LNSTL)).LT.IWORK(LMXSTP))
+     *   GO TO 510
+           IDID=-1
+           GO TO 527
+C
+C     UPDATE WT
+510   CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,RWORK(LPHI),
+     *  RWORK(LWT),RPAR,IPAR)
+      DO 520 I=1,NEQ
+         IF(RWORK(I+LWT-1).GT.0.0D0)GO TO 520
+           IDID=-3
+           GO TO 527
+520   CONTINUE
+C
+C     TEST FOR TOO MUCH ACCURACY REQUESTED.
+      R=DDANRM(NEQ,RWORK(LPHI),RWORK(LWT),RPAR,IPAR)*
+     *   100.0D0*UROUND
+      IF(R.LE.1.0D0)GO TO 525
+C     MULTIPLY RTOL AND ATOL BY R AND RETURN
+      IF(INFO(2).EQ.1)GO TO 523
+           RTOL(1)=R*RTOL(1)
+           ATOL(1)=R*ATOL(1)
+           IDID=-2
+           GO TO 527
+523   DO 524 I=1,NEQ
+           RTOL(I)=R*RTOL(I)
+524        ATOL(I)=R*ATOL(I)
+      IDID=-2
+      GO TO 527
+525   CONTINUE
+C
+C     COMPUTE MINIMUM STEPSIZE
+      HMIN=4.0D0*UROUND*DMAX1(DABS(TN),DABS(TOUT))
+C
+C     TEST H VS. HMAX
+      IF (INFO(7) .EQ. 0) GO TO 526
+         RH = ABS(H)/RWORK(LHMAX)
+         IF (RH .GT. 1.0D0) H = H/RH
+526   CONTINUE     
+C
+      CALL DDASTP(TN,Y,YPRIME,NEQ,
+     *   RES,JAC,H,RWORK(LWT),INFO(1),IDID,RPAR,IPAR,
+     *   RWORK(LPHI),RWORK(LDELTA),RWORK(LE),
+     *   RWORK(LWM),IWORK(LIWM),
+     *   RWORK(LALPHA),RWORK(LBETA),RWORK(LGAMMA),
+     *   RWORK(LPSI),RWORK(LSIGMA),
+     *   RWORK(LCJ),RWORK(LCJOLD),RWORK(LHOLD),
+     *   RWORK(LS),HMIN,RWORK(LROUND),
+     *   IWORK(LPHASE),IWORK(LJCALC),IWORK(LK),
+     *   IWORK(LKOLD),IWORK(LNS),INFO(10),NTEMP)
+527   IF(IDID.LT.0)GO TO 600
+C
+C--------------------------------------------------------
+C     THIS BLOCK HANDLES THE CASE OF A SUCCESSFUL RETURN
+C     FROM DDASTP (IDID=1).  TEST FOR STOP CONDITIONS.
+C--------------------------------------------------------
+C
+      IF(NG .EQ. 0) GO TO 529
+C
+C     Check for a zero of G near TN.
+C
+      CALL DRCHEK(3,G,NG,NEQ,TN,TOUT,Y,RWORK(LE),RWORK(LPHI),
+     *  RWORK(LPSI),IWORK(LKOLD),RWORK(LG0),RWORK(LG1),
+     *  RWORK(LGX),JROOT,IRT,RWORK(LROUND),INFO(3),
+     *  RWORK,IWORK,RPAR,IPAR)
+      IF(IRT .NE. 1) GO TO 529
+      IWORK(LIRFND) = 1
+      IDID = 4
+      T = RWORK(LT0)
+      GO TO 580
+C
+529   CONTINUE
+      IF(INFO(4).NE.0)GO TO 540
+           IF(INFO(3).NE.0)GO TO 530
+             IF((TN-TOUT)*H.LT.0.0D0)GO TO 500
+             CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
+     *         IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
+             IDID=3
+             T=TOUT
+             GO TO 580
+530          IF((TN-TOUT)*H.GE.0.0D0)GO TO 535
+             T=TN
+             IDID=1
+             GO TO 580
+535          CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
+     *         IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
+             IDID=3
+             T=TOUT
+             GO TO 580
+540   IF(INFO(3).NE.0)GO TO 550
+      IF((TN-TOUT)*H.LT.0.0D0)GO TO 542
+         CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
+     *     IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
+         T=TOUT
+         IDID=3
+         GO TO 580
+542   IF(DABS(TN-TSTOP).LE.100.0D0*UROUND*
+     *   (DABS(TN)+DABS(H)))GO TO 545
+      TNEXT=TN+H
+      IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 500
+      H=TSTOP-TN
+      GO TO 500
+545   CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,
+     *  IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
+      IDID=2
+      T=TSTOP
+      GO TO 580
+550   IF((TN-TOUT)*H.GE.0.0D0)GO TO 555
+      IF(DABS(TN-TSTOP).LE.100.0D0*UROUND*(DABS(TN)+DABS(H)))GO TO 552
+      T=TN
+      IDID=1
+      GO TO 580
+552   CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,
+     *  IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
+      IDID=2
+      T=TSTOP
+      GO TO 580
+555   CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
+     *   IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
+      T=TOUT
+      IDID=3
+580   CONTINUE
+C
+C--------------------------------------------------------
+C     ALL SUCCESSFUL RETURNS FROM DDASRT ARE MADE FROM
+C     THIS BLOCK.
+C--------------------------------------------------------
+C
+590   CONTINUE
+      RWORK(LTN)=TN
+      RWORK(LH)=H
+      RWORK(LTLAST) = T
+      RETURN
+C
+C-----------------------------------------------------------------------
+C     THIS BLOCK HANDLES ALL UNSUCCESSFUL
+C     RETURNS OTHER THAN FOR ILLEGAL INPUT.
+C-----------------------------------------------------------------------
+C
+600   CONTINUE
+      ITEMP=-IDID
+      GO TO (610,620,630,690,690,640,650,660,670,675,
+     *  680,685), ITEMP
+C
+C     THE MAXIMUM NUMBER OF STEPS WAS TAKEN BEFORE
+C     REACHING TOUT
+610   MSG = 'DASRT--  AT CURRENT T (=R1)  500 STEPS'
+      CALL XERRWV(MSG,38,610,0,0,0,0,1,TN,0.0D0)
+      MSG = 'DASRT--  TAKEN ON THIS CALL BEFORE REACHING TOUT'
+      CALL XERRWV(MSG,48,611,0,0,0,0,0,0.0D0,0.0D0)
+      GO TO 690
+C
+C     TOO MUCH ACCURACY FOR MACHINE PRECISION
+620   MSG = 'DASRT--  AT T (=R1) TOO MUCH ACCURACY REQUESTED'
+      CALL XERRWV(MSG,47,620,0,0,0,0,1,TN,0.0D0)
+      MSG = 'DASRT--  FOR PRECISION OF MACHINE. RTOL AND ATOL'
+      CALL XERRWV(MSG,48,621,0,0,0,0,0,0.0D0,0.0D0)
+      MSG = 'DASRT--  WERE INCREASED TO APPROPRIATE VALUES'
+      CALL XERRWV(MSG,45,622,0,0,0,0,0,0.0D0,0.0D0)
+C
+      GO TO 690
+C     WT(I) .LE. 0.0D0 FOR SOME I (NOT AT START OF PROBLEM)
+630   MSG = 'DASRT--  AT T (=R1) SOME ELEMENT OF WT'
+      CALL XERRWV(MSG,38,630,0,0,0,0,1,TN,0.0D0)
+      MSG = 'DASRT--  HAS BECOME .LE. 0.0'
+      CALL XERRWV(MSG,28,631,0,0,0,0,0,0.0D0,0.0D0)
+      GO TO 690
+C
+C     ERROR TEST FAILED REPEATEDLY OR WITH H=HMIN
+640   MSG = 'DASRT--  AT T (=R1) AND STEPSIZE H (=R2) THE'
+      CALL XERRWV(MSG,44,640,0,0,0,0,2,TN,H)
+      MSG='DASRT--  ERROR TEST FAILED REPEATEDLY OR WITH ABS(H)=HMIN'
+      CALL XERRWV(MSG,57,641,0,0,0,0,0,0.0D0,0.0D0)
+      GO TO 690
+C
+C     CORRECTOR CONVERGENCE FAILED REPEATEDLY OR WITH H=HMIN
+650   MSG = 'DASRT--  AT T (=R1) AND STEPSIZE H (=R2) THE'
+      CALL XERRWV(MSG,44,650,0,0,0,0,2,TN,H)
+      MSG = 'DASRT--  CORRECTOR FAILED TO CONVERGE REPEATEDLY'
+      CALL XERRWV(MSG,48,651,0,0,0,0,0,0.0D0,0.0D0)
+      MSG = 'DASRT--  OR WITH ABS(H)=HMIN'
+      CALL XERRWV(MSG,28,652,0,0,0,0,0,0.0D0,0.0D0)
+      GO TO 690
+C
+C     THE ITERATION MATRIX IS SINGULAR
+660   MSG = 'DASRT--  AT T (=R1) AND STEPSIZE H (=R2) THE'
+      CALL XERRWV(MSG,44,660,0,0,0,0,2,TN,H)
+      MSG = 'DASRT--  ITERATION MATRIX IS SINGULAR'
+      CALL XERRWV(MSG,37,661,0,0,0,0,0,0.0D0,0.0D0)
+      GO TO 690
+C
+C     CORRECTOR FAILURE PRECEEDED BY ERROR TEST FAILURES.
+670   MSG = 'DASRT--  AT T (=R1) AND STEPSIZE H (=R2) THE'
+      CALL XERRWV(MSG,44,670,0,0,0,0,2,TN,H)
+      MSG = 'DASRT--  CORRECTOR COULD NOT CONVERGE.  ALSO, THE'
+      CALL XERRWV(MSG,49,671,0,0,0,0,0,0.0D0,0.0D0)
+      MSG = 'DASRT--  ERROR TEST FAILED REPEATEDLY.'
+      CALL XERRWV(MSG,38,672,0,0,0,0,0,0.0D0,0.0D0)
+      GO TO 690
+C
+C     CORRECTOR FAILURE BECAUSE IRES = -1
+675   MSG = 'DASRT--  AT T (=R1) AND STEPSIZE H (=R2) THE'
+      CALL XERRWV(MSG,44,675,0,0,0,0,2,TN,H)
+      MSG = 'DASRT--  CORRECTOR COULD NOT CONVERGE BECAUSE'
+      CALL XERRWV(MSG,45,676,0,0,0,0,0,0.0D0,0.0D0)
+      MSG = 'DASRT--  IRES WAS EQUAL TO MINUS ONE'
+      CALL XERRWV(MSG,36,677,0,0,0,0,0,0.0D0,0.0D0)
+      GO TO 690
+C
+C     FAILURE BECAUSE IRES = -2
+680   MSG = 'DASRT--  AT T (=R1) AND STEPSIZE H (=R2)'
+      CALL XERRWV(MSG,40,680,0,0,0,0,2,TN,H)
+      MSG = 'DASRT--  IRES WAS EQUAL TO MINUS TWO'
+      CALL XERRWV(MSG,36,681,0,0,0,0,0,0.0D0,0.0D0)
+      GO TO 690
+C
+C     FAILED TO COMPUTE INITIAL YPRIME
+685   MSG = 'DASRT--  AT T (=R1) AND STEPSIZE H (=R2) THE'
+      CALL XERRWV(MSG,44,685,0,0,0,0,2,TN,HO)
+      MSG = 'DASRT--  INITIAL YPRIME COULD NOT BE COMPUTED'
+      CALL XERRWV(MSG,45,686,0,0,0,0,0,0.0D0,0.0D0)
+      GO TO 690
+690   CONTINUE
+      INFO(1)=-1
+      T=TN
+      RWORK(LTN)=TN
+      RWORK(LH)=H
+      RETURN
+C-----------------------------------------------------------------------
+C     THIS BLOCK HANDLES ALL ERROR RETURNS DUE
+C     TO ILLEGAL INPUT, AS DETECTED BEFORE CALLING
+C     DDASTP. FIRST THE ERROR MESSAGE ROUTINE IS
+C     CALLED. IF THIS HAPPENS TWICE IN
+C     SUCCESSION, EXECUTION IS TERMINATED
+C
+C-----------------------------------------------------------------------
+701   MSG = 'DASRT--  SOME ELEMENT OF INFO VECTOR IS NOT ZERO OR ONE'
+      CALL XERRWV(MSG,55,1,0,0,0,0,0,0.0D0,0.0D0)
+      GO TO 750
+702   MSG = 'DASRT--  NEQ (=I1) .LE. 0'
+      CALL XERRWV(MSG,25,2,0,1,NEQ,0,0,0.0D0,0.0D0)
+      GO TO 750
+703   MSG = 'DASRT--  MAXORD (=I1) NOT IN RANGE'
+      CALL XERRWV(MSG,34,3,0,1,MXORD,0,0,0.0D0,0.0D0)
+      GO TO 750
+704   MSG='DASRT--  RWORK LENGTH NEEDED, LENRW (=I1), EXCEEDS LRW (=I2)'
+      CALL XERRWV(MSG,60,4,0,2,LENRW,LRW,0,0.0D0,0.0D0)
+      GO TO 750
+705   MSG='DASRT--  IWORK LENGTH NEEDED, LENIW (=I1), EXCEEDS LIW (=I2)'
+      CALL XERRWV(MSG,60,5,0,2,LENIW,LIW,0,0.0D0,0.0D0)
+      GO TO 750
+706   MSG = 'DASRT--  SOME ELEMENT OF RTOL IS .LT. 0'
+      CALL XERRWV(MSG,39,6,0,0,0,0,0,0.0D0,0.0D0)
+      GO TO 750
+707   MSG = 'DASRT--  SOME ELEMENT OF ATOL IS .LT. 0'
+      CALL XERRWV(MSG,39,7,0,0,0,0,0,0.0D0,0.0D0)
+      GO TO 750
+708   MSG = 'DASRT--  ALL ELEMENTS OF RTOL AND ATOL ARE ZERO'
+      CALL XERRWV(MSG,47,8,0,0,0,0,0,0.0D0,0.0D0)
+      GO TO 750
+709   MSG='DASRT--  INFO(4) = 1 AND TSTOP (=R1) BEHIND TOUT (=R2)'
+      CALL XERRWV(MSG,54,9,0,0,0,0,2,TSTOP,TOUT)
+      GO TO 750
+710   MSG = 'DASRT--  HMAX (=R1) .LT. 0.0'
+      CALL XERRWV(MSG,28,10,0,0,0,0,1,HMAX,0.0D0)
+      GO TO 750
+711   MSG = 'DASRT--  TOUT (=R1) BEHIND T (=R2)'
+      CALL XERRWV(MSG,34,11,0,0,0,0,2,TOUT,T)
+      GO TO 750
+712   MSG = 'DASRT--  INFO(8)=1 AND H0=0.0'
+      CALL XERRWV(MSG,29,12,0,0,0,0,0,0.0D0,0.0D0)
+      GO TO 750
+713   MSG = 'DASRT--  SOME ELEMENT OF WT IS .LE. 0.0'
+      CALL XERRWV(MSG,39,13,0,0,0,0,0,0.0D0,0.0D0)
+      GO TO 750
+714   MSG='DASRT-- TOUT (=R1) TOO CLOSE TO T (=R2) TO START INTEGRATION'
+      CALL XERRWV(MSG,60,14,0,0,0,0,2,TOUT,T)
+      GO TO 750
+715   MSG = 'DASRT--  INFO(4)=1 AND TSTOP (=R1) BEHIND T (=R2)'
+      CALL XERRWV(MSG,49,15,0,0,0,0,2,TSTOP,T)
+      GO TO 750
+716   MSG = 'DASRT--  INFO(12)=1 AND MXSTP (=I1) .LT. 0'
+      CALL XERRWV(MSG,42,16,0,1,IWORK(LMXSTP),0,0,0.0D0,0.0D0)
+      GO TO 750
+717   MSG = 'DASRT--  ML (=I1) ILLEGAL. EITHER .LT. 0 OR .GT. NEQ'
+      CALL XERRWV(MSG,52,17,0,1,IWORK(LML),0,0,0.0D0,0.0D0)
+      GO TO 750
+718   MSG = 'DASRT--  MU (=I1) ILLEGAL. EITHER .LT. 0 OR .GT. NEQ'
+      CALL XERRWV(MSG,52,18,0,1,IWORK(LMU),0,0,0.0D0,0.0D0)
+      GO TO 750
+719   MSG = 'DASRT--  TOUT (=R1) IS EQUAL TO T (=R2)'
+      CALL XERRWV(MSG,39,19,0,0,0,0,2,TOUT,T)
+      GO TO 750
+730   MSG = 'DASRT--  NG (=I1) .LT. 0'
+      CALL XERRWV(MSG,24,30,1,1,NG,0,0,0.0D0,0.0D0)
+      GO TO 750
+732   MSG = 'DASRT--  ONE OR MORE COMPONENTS OF G HAS A ROOT'
+      CALL XERRWV(MSG,47,32,1,0,0,0,0,0.0D0,0.0D0)
+      MSG = '         TOO NEAR TO THE INITIAL POINT'
+      CALL XERRWV(MSG,38,32,1,0,0,0,0,0.0D0,0.0D0)
+750   IF(INFO(1).EQ.-1) GO TO 760
+      INFO(1)=-1
+      IDID=-33
+      RETURN
+760   MSG = 'DASRT--  REPEATED OCCURRENCES OF ILLEGAL INPUT'
+      CALL XERRWV(MSG,46,801,0,0,0,0,0,0.0D0,0.0D0)
+770   MSG = 'DASRT--  RUN TERMINATED. APPARENT INFINITE LOOP'
+      CALL XERRWV(MSG,47,802,1,0,0,0,0,0.0D0,0.0D0)
+      RETURN
+C-----------END OF SUBROUTINE DDASRT------------------------------------
+      END
new file mode 100644
--- /dev/null
+++ b/libcruft/dasrt/drchek.f
@@ -0,0 +1,172 @@
+      SUBROUTINE DRCHEK (JOB, G, NG, NEQ, TN, TOUT, Y, YP, PHI, PSI,
+     *  KOLD, G0, G1, GX, JROOT, IRT, UROUND, INFO3, RWORK, IWORK,
+     *  RPAR, IPAR)
+C
+C***BEGIN PROLOGUE  DRCHEK
+C***REFER TO DDASRT
+C***ROUTINES CALLED  DDATRP, DROOTS, DCOPY
+C***DATE WRITTEN   821001   (YYMMDD)
+C***REVISION DATE  900926   (YYMMDD)
+C***END PROLOGUE  DRCHEK
+C
+      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
+      PARAMETER (LNGE=16, LIRFND=18, LLAST=19, LIMAX=20,
+     *           LT0=41, LTLAST=42, LALPHR=43, LX2=44)
+      EXTERNAL G
+      INTEGER JOB, NG, NEQ, KOLD, JROOT, IRT, INFO3, IWORK, IPAR
+      DOUBLE PRECISION TN, TOUT, Y, YP, PHI, PSI, G0, G1, GX, UROUND,
+     *  RWORK, RPAR
+      DIMENSION  Y(*), YP(*), PHI(NEQ,*), PSI(*),
+     1  G0(*), G1(*), GX(*), JROOT(*), RWORK(*), IWORK(*)
+      INTEGER I, JFLAG
+      DOUBLE PRECISION H
+      DOUBLE PRECISION HMING, T1, TEMP1, TEMP2, X
+      LOGICAL ZROOT
+C-----------------------------------------------------------------------
+C THIS ROUTINE CHECKS FOR THE PRESENCE OF A ROOT IN THE
+C VICINITY OF THE CURRENT T, IN A MANNER DEPENDING ON THE
+C INPUT FLAG JOB.  IT CALLS SUBROUTINE DROOTS TO LOCATE THE ROOT
+C AS PRECISELY AS POSSIBLE.
+C
+C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, DRCHEK
+C USES THE FOLLOWING FOR COMMUNICATION..
+C JOB    = INTEGER FLAG INDICATING TYPE OF CALL..
+C          JOB = 1 MEANS THE PROBLEM IS BEING INITIALIZED, AND DRCHEK
+C                  IS TO LOOK FOR A ROOT AT OR VERY NEAR THE INITIAL T.
+C          JOB = 2 MEANS A CONTINUATION CALL TO THE SOLVER WAS JUST
+C                  MADE, AND DRCHEK IS TO CHECK FOR A ROOT IN THE
+C                  RELEVANT PART OF THE STEP LAST TAKEN.
+C          JOB = 3 MEANS A SUCCESSFUL STEP WAS JUST TAKEN, AND DRCHEK
+C                  IS TO LOOK FOR A ROOT IN THE INTERVAL OF THE STEP.
+C G0     = ARRAY OF LENGTH NG, CONTAINING THE VALUE OF G AT T = T0.
+C          G0 IS INPUT FOR JOB .GE. 2 AND ON OUTPUT IN ALL CASES.
+C G1,GX  = ARRAYS OF LENGTH NG FOR WORK SPACE.
+C IRT    = COMPLETION FLAG..
+C          IRT = 0  MEANS NO ROOT WAS FOUND.
+C          IRT = -1 MEANS JOB = 1 AND A ROOT WAS FOUND TOO NEAR TO T.
+C          IRT = 1  MEANS A LEGITIMATE ROOT WAS FOUND (JOB = 2 OR 3).
+C                   ON RETURN, T0 IS THE ROOT LOCATION, AND Y IS THE
+C                   CORRESPONDING SOLUTION VECTOR.
+C T0     = VALUE OF T AT ONE ENDPOINT OF INTERVAL OF INTEREST.  ONLY
+C          ROOTS BEYOND T0 IN THE DIRECTION OF INTEGRATION ARE SOUGHT.
+C          T0 IS INPUT IF JOB .GE. 2, AND OUTPUT IN ALL CASES.
+C          T0 IS UPDATED BY DRCHEK, WHETHER A ROOT IS FOUND OR NOT.
+C          STORED IN THE GLOBAL ARRAY RWORK.
+C TLAST  = LAST VALUE OF T RETURNED BY THE SOLVER (INPUT ONLY).
+C          STORED IN THE GLOBAL ARRAY RWORK.
+C TOUT   = FINAL OUTPUT TIME FOR THE SOLVER.
+C IRFND  = INPUT FLAG SHOWING WHETHER THE LAST STEP TAKEN HAD A ROOT.
+C          IRFND = 1 IF IT DID, = 0 IF NOT.
+C          STORED IN THE GLOBAL ARRAY IWORK.
+C INFO3  = COPY OF INFO(3) (INPUT ONLY).
+C-----------------------------------------------------------------------
+C     
+      H = PSI(1)
+      IRT = 0
+      DO 10 I = 1,NG
+ 10     JROOT(I) = 0
+      HMING = (DABS(TN) + DABS(H))*UROUND*100.0D0
+C
+      GO TO (100, 200, 300), JOB
+C
+C EVALUATE G AT INITIAL T (STORED IN RWORK(LT0)), AND CHECK FOR
+C ZERO VALUES.----------------------------------------------------------
+ 100  CONTINUE
+      CALL DDATRP(TN,RWORK(LT0),Y,YP,NEQ,KOLD,PHI,PSI)
+      CALL G (NEQ, RWORK(LT0), Y, NG, G0, RPAR, IPAR)
+      IWORK(LNGE) = 1
+      ZROOT = .FALSE.
+      DO 110 I = 1,NG
+ 110    IF (DABS(G0(I)) .LE. 0.0D0) ZROOT = .TRUE.
+      IF (.NOT. ZROOT) GO TO 190
+C G HAS A ZERO AT T.  LOOK AT G AT T + (SMALL INCREMENT). --------------
+      TEMP1 = DSIGN(HMING,H)
+      RWORK(LT0) = RWORK(LT0) + TEMP1
+      TEMP2 = TEMP1/H
+      DO 120 I = 1,NEQ
+ 120    Y(I) = Y(I) + TEMP2*PHI(I,2)
+      CALL G (NEQ, RWORK(LT0), Y, NG, G0, RPAR, IPAR)
+      IWORK(LNGE) = IWORK(LNGE) + 1
+      ZROOT = .FALSE.
+      DO 130 I = 1,NG
+ 130    IF (DABS(G0(I)) .LE. 0.0D0) ZROOT = .TRUE.
+      IF (.NOT. ZROOT) GO TO 190
+C G HAS A ZERO AT T AND ALSO CLOSE TO T.  TAKE ERROR RETURN. -----------
+      IRT = -1
+      RETURN
+C
+ 190  CONTINUE
+      RETURN
+C
+C
+ 200  CONTINUE
+      IF (IWORK(LIRFND) .EQ. 0) GO TO 260
+C IF A ROOT WAS FOUND ON THE PREVIOUS STEP, EVALUATE G0 = G(T0). -------
+      CALL DDATRP (TN, RWORK(LT0), Y, YP, NEQ, KOLD, PHI, PSI)
+      CALL G (NEQ, RWORK(LT0), Y, NG, G0, RPAR, IPAR)
+      IWORK(LNGE) = IWORK(LNGE) + 1
+      ZROOT = .FALSE.
+      DO 210 I = 1,NG
+ 210    IF (DABS(G0(I)) .LE. 0.0D0) ZROOT = .TRUE.
+      IF (.NOT. ZROOT) GO TO 260
+C G HAS A ZERO AT T0.  LOOK AT G AT T + (SMALL INCREMENT). -------------
+      TEMP1 = DSIGN(HMING,H)
+      RWORK(LT0) = RWORK(LT0) + TEMP1
+      IF ((RWORK(LT0) - TN)*H .LT. 0.0D0) GO TO 230
+      TEMP2 = TEMP1/H
+      DO 220 I = 1,NEQ
+ 220    Y(I) = Y(I) + TEMP2*PHI(I,2)
+      GO TO 240
+ 230  CALL DDATRP (TN, RWORK(LT0), Y, YP, NEQ, KOLD, PHI, PSI)
+ 240  CALL G (NEQ, RWORK(LT0), Y, NG, G0, RPAR, IPAR)
+      IWORK(LNGE) = IWORK(LNGE) + 1
+      ZROOT = .FALSE.
+      DO 250 I = 1,NG
+        IF (DABS(G0(I)) .GT. 0.0D0) GO TO 250
+        JROOT(I) = 1
+        ZROOT = .TRUE.
+ 250    CONTINUE
+      IF (.NOT. ZROOT) GO TO 260
+C G HAS A ZERO AT T0 AND ALSO CLOSE TO T0.  RETURN ROOT. ---------------
+      IRT = 1
+      RETURN
+C     HERE, G0 DOES NOT HAVE A ROOT
+C G0 HAS NO ZERO COMPONENTS.  PROCEED TO CHECK RELEVANT INTERVAL. ------
+ 260  IF (TN .EQ. RWORK(LTLAST)) GO TO 390
+C
+ 300  CONTINUE
+C SET T1 TO TN OR TOUT, WHICHEVER COMES FIRST, AND GET G AT T1. --------
+      IF (INFO3 .EQ. 1) GO TO 310
+      IF ((TOUT - TN)*H .GE. 0.0D0) GO TO 310
+      T1 = TOUT
+      IF ((T1 - RWORK(LT0))*H .LE. 0.0D0) GO TO 390
+      CALL DDATRP (TN, T1, Y, YP, NEQ, KOLD, PHI, PSI)
+      GO TO 330
+ 310  T1 = TN
+      DO 320 I = 1,NEQ
+ 320    Y(I) = PHI(I,1)
+ 330  CALL G (NEQ, T1, Y, NG, G1, RPAR, IPAR)
+      IWORK(LNGE) = IWORK(LNGE) + 1
+C CALL DROOTS TO SEARCH FOR ROOT IN INTERVAL FROM T0 TO T1. ------------
+      JFLAG = 0
+ 350  CONTINUE
+      CALL DROOTS (NG, HMING, JFLAG, RWORK(LT0), T1, G0, G1, GX, X,
+     *             JROOT, IWORK(LIMAX), IWORK(LLAST), RWORK(LALPHR),
+     *             RWORK(LX2))
+      IF (JFLAG .GT. 1) GO TO 360
+      CALL DDATRP (TN, X, Y, YP, NEQ, KOLD, PHI, PSI)
+      CALL G (NEQ, X, Y, NG, GX, RPAR, IPAR)
+      IWORK(LNGE) = IWORK(LNGE) + 1
+      GO TO 350
+ 360  RWORK(LT0) = X
+      CALL DCOPY (NG, GX, 1, G0, 1)
+      IF (JFLAG .EQ. 4) GO TO 390
+C FOUND A ROOT.  INTERPOLATE TO X AND RETURN. --------------------------
+      CALL DDATRP (TN, X, Y, YP, NEQ, KOLD, PHI, PSI)
+      IRT = 1
+      RETURN
+C
+ 390  CONTINUE
+      RETURN
+C---------------------- END OF SUBROUTINE DRCHEK -----------------------
+      END
new file mode 100644
--- /dev/null
+++ b/libcruft/dasrt/droots.f
@@ -0,0 +1,217 @@
+      SUBROUTINE DROOTS (NG, HMIN, JFLAG, X0, X1, G0, G1, GX, X, JROOT,
+     *                   IMAX, LAST, ALPHA, X2)
+C
+C***BEGIN PROLOGUE  DROOTS
+C***REFER TO DDASRT
+C***ROUTINES CALLED  DCOPY
+C***DATE WRITTEN   821001   (YYMMDD)
+C***REVISION DATE  900926   (YYMMDD)
+C***END PROLOGUE  DROOTS
+C
+      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
+      INTEGER NG, JFLAG, JROOT, IMAX, LAST
+      DOUBLE PRECISION HMIN, X0, X1, G0, G1, GX, X, ALPHA, X2
+      DIMENSION G0(NG), G1(NG), GX(NG), JROOT(NG)
+C-----------------------------------------------------------------------
+C THIS SUBROUTINE FINDS THE LEFTMOST ROOT OF A SET OF ARBITRARY
+C FUNCTIONS GI(X) (I = 1,...,NG) IN AN INTERVAL (X0,X1).  ONLY ROOTS
+C OF ODD MULTIPLICITY (I.E. CHANGES OF SIGN OF THE GI) ARE FOUND.
+C HERE THE SIGN OF X1 - X0 IS ARBITRARY, BUT IS CONSTANT FOR A GIVEN
+C PROBLEM, AND -LEFTMOST- MEANS NEAREST TO X0.
+C THE VALUES OF THE VECTOR-VALUED FUNCTION G(X) = (GI, I=1...NG)
+C ARE COMMUNICATED THROUGH THE CALL SEQUENCE OF DROOTS.
+C THE METHOD USED IS THE ILLINOIS ALGORITHM.
+C
+C REFERENCE..
+C KATHIE L. HIEBERT AND LAWRENCE F. SHAMPINE, IMPLICITLY DEFINED
+C OUTPUT POINTS FOR SOLUTIONS OF ODE-S, SANDIA REPORT SAND80-0180,
+C FEBRUARY, 1980.
+C
+C DESCRIPTION OF PARAMETERS.
+C
+C NG     = NUMBER OF FUNCTIONS GI, OR THE NUMBER OF COMPONENTS OF
+C          THE VECTOR VALUED FUNCTION G(X).  INPUT ONLY.
+C
+C HMIN   = RESOLUTION PARAMETER IN X.  INPUT ONLY.  WHEN A ROOT IS
+C          FOUND, IT IS LOCATED ONLY TO WITHIN AN ERROR OF HMIN IN X.
+C          TYPICALLY, HMIN SHOULD BE SET TO SOMETHING ON THE ORDER OF
+C               100 * UROUND * MAX(ABS(X0),ABS(X1)),
+C          WHERE UROUND IS THE UNIT ROUNDOFF OF THE MACHINE.
+C
+C JFLAG  = INTEGER FLAG FOR INPUT AND OUTPUT COMMUNICATION.
+C
+C          ON INPUT, SET JFLAG = 0 ON THE FIRST CALL FOR THE PROBLEM,
+C          AND LEAVE IT UNCHANGED UNTIL THE PROBLEM IS COMPLETED.
+C          (THE PROBLEM IS COMPLETED WHEN JFLAG .GE. 2 ON RETURN.)
+C
+C          ON OUTPUT, JFLAG HAS THE FOLLOWING VALUES AND MEANINGS..
+C          JFLAG = 1 MEANS DROOTS NEEDS A VALUE OF G(X).  SET GX = G(X)
+C                    AND CALL DROOTS AGAIN.
+C          JFLAG = 2 MEANS A ROOT HAS BEEN FOUND.  THE ROOT IS
+C                    AT X, AND GX CONTAINS G(X).  (ACTUALLY, X IS THE
+C                    RIGHTMOST APPROXIMATION TO THE ROOT ON AN INTERVAL
+C                    (X0,X1) OF SIZE HMIN OR LESS.)
+C          JFLAG = 3 MEANS X = X1 IS A ROOT, WITH ONE OR MORE OF THE GI
+C                    BEING ZERO AT X1 AND NO SIGN CHANGES IN (X0,X1).
+C                    GX CONTAINS G(X) ON OUTPUT.
+C          JFLAG = 4 MEANS NO ROOTS (OF ODD MULTIPLICITY) WERE
+C                    FOUND IN (X0,X1) (NO SIGN CHANGES).
+C
+C X0,X1  = ENDPOINTS OF THE INTERVAL WHERE ROOTS ARE SOUGHT.
+C          X1 AND X0 ARE INPUT WHEN JFLAG = 0 (FIRST CALL), AND
+C          MUST BE LEFT UNCHANGED BETWEEN CALLS UNTIL THE PROBLEM IS
+C          COMPLETED.  X0 AND X1 MUST BE DISTINCT, BUT X1 - X0 MAY BE
+C          OF EITHER SIGN.  HOWEVER, THE NOTION OF -LEFT- AND -RIGHT-
+C          WILL BE USED TO MEAN NEARER TO X0 OR X1, RESPECTIVELY.
+C          WHEN JFLAG .GE. 2 ON RETURN, X0 AND X1 ARE OUTPUT, AND
+C          ARE THE ENDPOINTS OF THE RELEVANT INTERVAL.
+C
+C G0,G1  = ARRAYS OF LENGTH NG CONTAINING THE VECTORS G(X0) AND G(X1),
+C          RESPECTIVELY.  WHEN JFLAG = 0, G0 AND G1 ARE INPUT AND
+C          NONE OF THE G0(I) SHOULD BE BE ZERO.
+C          WHEN JFLAG .GE. 2 ON RETURN, G0 AND G1 ARE OUTPUT.
+C
+C GX     = ARRAY OF LENGTH NG CONTAINING G(X).  GX IS INPUT
+C          WHEN JFLAG = 1, AND OUTPUT WHEN JFLAG .GE. 2.
+C
+C X      = INDEPENDENT VARIABLE VALUE.  OUTPUT ONLY.
+C          WHEN JFLAG = 1 ON OUTPUT, X IS THE POINT AT WHICH G(X)
+C          IS TO BE EVALUATED AND LOADED INTO GX.
+C          WHEN JFLAG = 2 OR 3, X IS THE ROOT.
+C          WHEN JFLAG = 4, X IS THE RIGHT ENDPOINT OF THE INTERVAL, X1.
+C
+C JROOT  = INTEGER ARRAY OF LENGTH NG.  OUTPUT ONLY.
+C          WHEN JFLAG = 2 OR 3, JROOT INDICATES WHICH COMPONENTS
+C          OF G(X) HAVE A ROOT AT X.  JROOT(I) IS 1 IF THE I-TH
+C          COMPONENT HAS A ROOT, AND JROOT(I) = 0 OTHERWISE.
+C
+C IMAX, LAST, ALPHA, X2 =
+C          BOOKKEEPING VARIABLES WHICH MUST BE SAVED FROM CALL
+C          TO CALL.  THEY ARE SAVED INSIDE THE CALLING ROUTINE,
+C          BUT THEY ARE USED ONLY WITHIN THIS ROUTINE.
+C-----------------------------------------------------------------------
+      INTEGER I, IMXOLD, NXLAST
+      DOUBLE PRECISION T2, TMAX, ZERO
+      LOGICAL ZROOT, SGNCHG, XROOT
+      DATA ZERO/0.0D0/
+C
+      IF (JFLAG .EQ. 1) GO TO 200
+C JFLAG .NE. 1.  CHECK FOR CHANGE IN SIGN OF G OR ZERO AT X1. ----------
+      IMAX = 0
+      TMAX = ZERO
+      ZROOT = .FALSE.
+      DO 120 I = 1,NG
+        IF (DABS(G1(I)) .GT. ZERO) GO TO 110
+        ZROOT = .TRUE.
+        GO TO 120
+C AT THIS POINT, G0(I) HAS BEEN CHECKED AND CANNOT BE ZERO. ------------
+ 110    IF (DSIGN(1.0D0,G0(I)) .EQ. DSIGN(1.0D0,G1(I))) GO TO 120
+          T2 = DABS(G1(I)/(G1(I)-G0(I)))
+          IF (T2 .LE. TMAX) GO TO 120
+            TMAX = T2
+            IMAX = I
+ 120    CONTINUE
+      IF (IMAX .GT. 0) GO TO 130
+      SGNCHG = .FALSE.
+      GO TO 140
+ 130  SGNCHG = .TRUE.
+ 140  IF (.NOT. SGNCHG) GO TO 400
+C THERE IS A SIGN CHANGE.  FIND THE FIRST ROOT IN THE INTERVAL. --------
+      XROOT = .FALSE.
+      NXLAST = 0
+      LAST = 1
+C
+C REPEAT UNTIL THE FIRST ROOT IN THE INTERVAL IS FOUND.  LOOP POINT. ---
+ 150  CONTINUE
+      IF (XROOT) GO TO 300
+      IF (NXLAST .EQ. LAST) GO TO 160
+      ALPHA = 1.0D0
+      GO TO 180
+ 160  IF (LAST .EQ. 0) GO TO 170
+      ALPHA = 0.5D0*ALPHA
+      GO TO 180
+ 170  ALPHA = 2.0D0*ALPHA
+ 180  X2 = X1 - (X1-X0)*G1(IMAX)/(G1(IMAX) - ALPHA*G0(IMAX))
+      IF ((DABS(X2-X0) .LT. HMIN) .AND.
+     1   (DABS(X1-X0) .GT. 10.0D0*HMIN)) X2 = X0 + 0.1D0*(X1-X0)
+      JFLAG = 1
+      X = X2
+C RETURN TO THE CALLING ROUTINE TO GET A VALUE OF GX = G(X). -----------
+      RETURN
+C CHECK TO SEE IN WHICH INTERVAL G CHANGES SIGN. -----------------------
+ 200  IMXOLD = IMAX
+      IMAX = 0
+      TMAX = ZERO
+      ZROOT = .FALSE.
+      DO 220 I = 1,NG
+        IF (DABS(GX(I)) .GT. ZERO) GO TO 210
+        ZROOT = .TRUE.
+        GO TO 220
+C NEITHER G0(I) NOR GX(I) CAN BE ZERO AT THIS POINT. -------------------
+ 210    IF (DSIGN(1.0D0,G0(I)) .EQ. DSIGN(1.0D0,GX(I))) GO TO 220
+          T2 = DABS(GX(I)/(GX(I) - G0(I)))
+          IF (T2 .LE. TMAX) GO TO 220
+            TMAX = T2
+            IMAX = I
+ 220    CONTINUE
+      IF (IMAX .GT. 0) GO TO 230
+      SGNCHG = .FALSE.
+      IMAX = IMXOLD
+      GO TO 240
+ 230  SGNCHG = .TRUE.
+ 240  NXLAST = LAST
+      IF (.NOT. SGNCHG) GO TO 250
+C SIGN CHANGE BETWEEN X0 AND X2, SO REPLACE X1 WITH X2. ----------------
+      X1 = X2
+      CALL DCOPY (NG, GX, 1, G1, 1)
+      LAST = 1
+      XROOT = .FALSE.
+      GO TO 270
+ 250  IF (.NOT. ZROOT) GO TO 260
+C ZERO VALUE AT X2 AND NO SIGN CHANGE IN (X0,X2), SO X2 IS A ROOT. -----
+      X1 = X2
+      CALL DCOPY (NG, GX, 1, G1, 1)
+      XROOT = .TRUE.
+      GO TO 270
+C NO SIGN CHANGE BETWEEN X0 AND X2.  REPLACE X0 WITH X2. ---------------
+ 260  CONTINUE
+      CALL DCOPY (NG, GX, 1, G0, 1)
+      X0 = X2
+      LAST = 0
+      XROOT = .FALSE.
+ 270  IF (DABS(X1-X0) .LE. HMIN) XROOT = .TRUE.
+      GO TO 150
+C
+C RETURN WITH X1 AS THE ROOT.  SET JROOT.  SET X = X1 AND GX = G1. -----
+ 300  JFLAG = 2
+      X = X1
+      CALL DCOPY (NG, G1, 1, GX, 1)
+      DO 320 I = 1,NG
+        JROOT(I) = 0
+        IF (DABS(G1(I)) .GT. ZERO) GO TO 310
+          JROOT(I) = 1
+          GO TO 320
+ 310    IF (DSIGN(1.0D0,G0(I)) .NE. DSIGN(1.0D0,G1(I))) JROOT(I) = 1
+ 320    CONTINUE
+      RETURN
+C
+C NO SIGN CHANGE IN THE INTERVAL.  CHECK FOR ZERO AT RIGHT ENDPOINT. ---
+ 400  IF (.NOT. ZROOT) GO TO 420
+C
+C ZERO VALUE AT X1 AND NO SIGN CHANGE IN (X0,X1).  RETURN JFLAG = 3. ---
+      X = X1
+      CALL DCOPY (NG, G1, 1, GX, 1)
+      DO 410 I = 1,NG
+        JROOT(I) = 0
+        IF (DABS(G1(I)) .LE. ZERO) JROOT (I) = 1
+ 410  CONTINUE
+      JFLAG = 3
+      RETURN
+C
+C NO SIGN CHANGES IN THIS INTERVAL.  SET X = X1, RETURN JFLAG = 4. -----
+ 420  CALL DCOPY (NG, G1, 1, GX, 1)
+      X = X1
+      JFLAG = 4
+      RETURN
+C---------------------- END OF SUBROUTINE DROOTS -----------------------
+      END
new file mode 100644
--- /dev/null
+++ b/libcruft/dasrt/xerrwv.f
@@ -0,0 +1,68 @@
+      SUBROUTINE XERRWV (MSG, NMES, NERR, LEVEL, NI, I1, I2, NR, R1, R2)
+      INTEGER NMES, NERR, LEVEL, NI, I1, I2, NR
+      DOUBLE PRECISION R1, R2
+      CHARACTER*1 MSG(NMES)
+C-----------------------------------------------------------------------
+C Subroutine XERRWV, as given here, constitutes a simplified version of
+C the SLATEC error handling package.
+C Written by A. C. Hindmarsh and P. N. Brown at LLNL.
+C Modified 1/8/90 by Clement Ulrich at LLNL.
+C Version of 8 January, 1990.
+C This version is in double precision.
+C
+C All arguments are input arguments.
+C
+C MSG    = The message (character array).
+C NMES   = The length of MSG (number of characters).
+C NERR   = The error number (not used).
+C LEVEL  = The error level..
+C          0 or 1 means recoverable (control returns to caller).
+C          2 means fatal (run is aborted--see note below).
+C NI     = Number of integers (0, 1, or 2) to be printed with message.
+C I1,I2  = Integers to be printed, depending on NI.
+C NR     = Number of reals (0, 1, or 2) to be printed with message.
+C R1,R2  = Reals to be printed, depending on NR.
+C
+C Note..  this routine is compatible with ANSI-77; however the
+C following assumptions may not be valid for some machines:
+C
+C 1. The argument MSG is assumed to be of type CHARACTER, and
+C    the message is printed with a format of (1X,80A1).
+C 2. The message is assumed to take only one line.
+C    Multi-line messages are generated by repeated calls.
+C 3. If LEVEL = 2, control passes to the statement STOP
+C    to abort the run.  For a different run-abort command,
+C    change the statement following statement 100 at the end.
+C 4. R1 and R2 are assumed to be in double precision and are printed
+C    in E21.13 format.
+C 5. The logical unit number 6 is standard output.
+C    For a different default logical unit number, change the assignment
+C    statement for LUNIT below.
+C
+C-----------------------------------------------------------------------
+C Subroutines called by XERRWV.. None
+C Function routines called by XERRWV.. None
+C-----------------------------------------------------------------------
+C
+      INTEGER I, LUNIT, MESFLG
+C
+C Define message print flag and logical unit number. -------------------
+      MESFLG = 1
+      LUNIT = 6
+      IF (MESFLG .EQ. 0) GO TO 100
+C Write the message. ---------------------------------------------------
+      WRITE (LUNIT,10) (MSG(I),I=1,NMES)
+ 10   FORMAT(1X,80A1)
+      IF (NI .EQ. 1) WRITE (LUNIT, 20) I1
+ 20   FORMAT(6X,'In above message,  I1 =',I10)
+      IF (NI .EQ. 2) WRITE (LUNIT, 30) I1,I2
+ 30   FORMAT(6X,'In above message,  I1 =',I10,3X,'I2 =',I10)
+      IF (NR .EQ. 1) WRITE (LUNIT, 40) R1
+ 40   FORMAT(6X,'In above message,  R1 =',E21.13)
+      IF (NR .EQ. 2) WRITE (LUNIT, 50) R1,R2
+ 50   FORMAT(6X,'In above,  R1 =',E21.13,3X,'R2 =',E21.13)
+C Abort the run if LEVEL = 2. ------------------------------------------
+ 100  IF (LEVEL .NE. 2) RETURN
+      STOP
+C----------------------- End of Subroutine XERRWV ----------------------
+      END