      SUBROUTINE DDASSL (RES,NEQ,T,Y,YPRIME,TOUT,
     *  INFO,RTOL,ATOL,IDID,
     *  RWORK,LRW,IWORK,LIW,RPAR,IPAR,
     *  JAC)
C
C***BEGIN PROLOGUE  DDASSL
C***DATE WRITTEN   830315   (YYMMDD)
C***REVISION DATE  830315   (YYMMDD)
C***CATEGORY NO.  D2A2
C***KEYWORDS  DIFFERENTIAL/ALGEBRAIC,BACKWARD DIFFERENTIATION FORMULAS
C             IMPLICIT DIFFERENTIAL SYSTEMS
C***AUTHOR  PETZOLD,LINDA R.
C             APPLIED MATHEMATICS DIVISION 8331
C             SANDIA NATIONAL LABORATORIES
C             LIVERMORE, CA.    94550
C***PURPOSE  DIFFERENTIAL/ALGEBRAIC SYSTEM SOLVER
C***DESCRIPTION
C  ---------------------------------------------------------------------
C
C  this code solves a system of differential/
C  algebraic equations of the form
C  g(t,y,yprime) = 0.
C
C  subroutine ddassl uses the backward
C  differentiation formulas of orders one
C  through five to solve a system of the above
C  form for y and yprime. values for y
C  and yprime at the initial time must
C  be given as input. these values must
C  be consistent, (that is. if t,y,yprime
C  are the given initial values, they must
C  satisfy g(t,y,yprime) = 0.)
C  the subroutine solves the system from t to tout. it is
C  easy to continue the solution to get results
C  at additional tout. this is the interval
C  mode of operation. intermediate results can
C  also be obtained easily by using the intermediate-
C  output capability.
C
C
C  ------------description of arguments to ddassl-----------------------
C  ------------(an overview)--------------------------------------------
C
C  the parameters are
C
C  res -- this is a subroutine which you provide
C         to define the differential/algebraic
C         system
C
C  neq -- this is the number of equations
C         to be solved
C
C  t -- this is the current value of the
C       independent variable.
C
C  tout -- this is a point at which a solution
C      is desired.
C
C  info(*) -- the basic task of the code is
C             to solve the system from t to
C             tout and return an answer at tout.
C             info(*) is an integer array which is
C             used to communicate exactly how you
C             want this task to be carried out.
C
C  y(*) -- this array contains the solution
C          components at t
C
C  yprime(*) -- this array contains the derivatives
C               of the solution components at t
C
C  rtol,atol -- these quantities represent
C               absolute and relative error
C               tolerances which you provide to indicate
C               how accurately you wish the solution
C               to be computed. you may choose them
C               to be both scalars or else both
C               vectors.
C
C  idid -- this scalar quantity is an indicator reporting
C          what the code did. you must monitor this
C          integer variable to decide what action to
C          take next.
C
C  rwork(*),lrw -- rwork(*) is a real work array of
C                  length lrw which provides the code
C                  with needed storage space.
C
C  iwork(*),liw -- iwork(*) is an integer work array
C                  of length liw which provides the code
C                  with needed storage space.
C
C  rpar,ipar -- these are real and integer parameter
C               arrays which you can use for
C               communication between your calling
C               program and the res subroutine
C               (and the jac subroutine)
C
C  jac -- this is the name of a subroutine which you
C         may choose to provide for defining
C         a matrix of partial derivatives
C         described below.
C
C  quantities which are used as input items are
C     neq,t,y(*),yprime(*),tout,info(*),
C     rtol,atol,rwork(1),rwork(2),rwork(3),lrw,iwork(1),
C     iwork(2),iwork(3),and liw.
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  ----------input-what to do on the first call to ddassl---------------
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 differential/algebraic
C         system
C             delta = g(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 ddassl.
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 ddassl will try to solve the problem
C         without getting ires = -1.  if ires = -2, ddassl
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 ddassl. 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 ddassl uses
C            only the first nine 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) - ddassl will perform much better if the matrix of
C              partial derivatives, dg/dy + cj*dg/dyprime,
C              (here cj is a scalar determined by ddassl)
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
C              the code works with a full matrix of neq**2 elements
C              (stored in the conventional way). computations with
C              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 ho.
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 ho by setting
C                      rwork(3)=ho ****
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 will
C               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) --ddassl normally requires the initial t,
C               y, and yprime to be consistent.  that is,
C               you must have g(t,y,yprime) = 0 at the initial
C               time.  if you do not know the initial
C               derivative precisely, you can let ddassl 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                       g(t,y,yprime) = 0.)
C
C   rtol, atol -- you must assign relative (rtol) and absolute (atol
C               error tolerances to tell the code how accurately you want
C               the solution to be computed. they must be defined as
C               variables because the code may change them. you have two
C               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 test
C               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 true
C               solution of the initial value problem and the computed
C               approximation. practically all present day codes.
C               including this one, control the local error at each step
C               and do not even attempt to control the global error
C               directly.
C               usually, but not always, the true accuracy of
C               the computed y is comparable to the error tolerances. this
C               code will usually, but not always, deliver a more accurate
C               solution if you reduce the tolerances and integrate again.
C               by comparing two such solutions you can get a fairly
C               reliable idea of the true error in the solution at the
C               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 absolute
C               error test on that component. a mixed test with non-zero
C               rtol and atol corresponds roughly to a relative error
C               test when the solution component is much bigger than atol
C               and to an absolute error test when the solution component
C               is smaller than the threshold atol.
C
C               the code will not attempt to compute a solution at an
C               accuracy unreasonable for the machine being used. it will
C               advise you if you ask for too much accuracy and inform
C               you as to the maximum accuracy it believes 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. 40+(maxord+4)*neq+neq**2
C               for the full (dense) jacobian case (when info(6)=0),  or
C                    lrw .ge. 40+(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. 40+(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 ddassl and the
C               res subroutine (and the jac subroutine). they are not
C               altered by ddassl. if you do not need rpar or ipar, ignore
C               these parameters by treating them as dummy arguments. if
C               you do choose to use them, dimension them in your calling
C               program and in res (and in jac) as arrays of appropriate
C               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               compowent, 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               ddassl. 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                   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) = * dg(i)/dy(j)+cj*dg(i)/dyprime(j)*
C               *** info(6)=1 -- banded jacobian with ml lower and mu
C                   upper diagonal bands (refer to info(6) description of
C                   ml and mu) ***
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) = *dg(i)/dy(j)+cj*dg(i)/dyprime(j)*
C               rpar and ipar are real and integer parameter arrays which
C               you can use for communication between your calling
C               program and your jacobian subroutine jac. they are not
C               altered by ddassl. if you do not need rpar or ipar, ignore
C               these parameters by treating them as dummy arguments. if
C               you do choose to use them, dimension them in your calling
C               program and in jac as arrays of appropriate length.
C
C
C
C  optionally replaceable norm routine:
C  ddassl 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 ddassl
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 ddassl----
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 tout was successfully
C                   completed (t=tout) by stepping exactly to tout.
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                    *** task interrupted ***
C                reported by negative values of idid
C
C           idid = -1 -- a large amount of work has been expended.
C                   (about 500 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 -- ddassl 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 -- ddassl 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 be
C               appropriate for continuing the integration. however, the
C               reported solution at t was obtained using the input values
C               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 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 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
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 beyound 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                    *** 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 taken about 500 steps.
C                  if you want to continue, set info(1) = 1 and
C                  call the code again. an additional 500 steps
C                  will be allowed.
C
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 ddassl. 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 ddassl. an inaccurate or
C                  illconditioned jacobian may be the problem. if you
C                  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                  ddassl cannot solve the problem as stated.
C                  it is possible that the redundant equations
C                  could be removed, and then ddassl 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, ddassl 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, ddassl 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, ddassl 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***REFERENCES  A DESCRIPTION OF DASSL: A DIFFERENTIAL/ALGEBRAIC
C                  SYSTEM SOLVER, L. R. PETZOLD, SAND82-8637,
C                  SANDIA NATIONAL LABORATORIES, SEPTEMBER 1982.
C***ROUTINES CALLED  DDASTP,DDAINI,DDANRM,DDAWTS,DDATRP,XERRWV,D1MACH
C***COMMON BLOCKS    DDA001
C***END PROLOGUE DDASSL
C
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL DONE
      EXTERNAL RES,JAC
      DIMENSION Y(1),YPRIME(1)
      DIMENSION INFO(15)
      DIMENSION RWORK(1),IWORK(1)
      DIMENSION RTOL(1),ATOL(1)
      DIMENSION RPAR(1),IPAR(1)
      COMMON/DDA001/NPD,NTEMP,
     *   LML,LMU,LMXORD,LMTYPE,
     *   LNST,LNRE,LNJE,LETF,LCTF,LIPVT
      DATA LTSTOP,LHMAX,LH,LTN,
     *   LCJ,LCJOLD,LHOLD,LS,LROUND,
     *   LALPHA,LBETA,LGAMMA,
     *   LPSI,LSIGMA,LDELTA
     *   /1,2,3,4,
     *   5,6,7,8,9,
     *   11,17,23,
     *   29,35,41/
      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,11
         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     set pointers into iwork
      LML=1
      LMU=2
      LMXORD=3
      LMTYPE=4
      LJCALC=5
      LPHASE=6
      LK=7
      LKOLD=8
      LNS=9
      LNSTL=10
      LNST=11
      LNRE=12
      LNJE=13
      LETF=14
      LCTF=15
      LIPVT=21
      LIWM=1
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=40+(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=40+(IWORK(LMXORD)+4)*NEQ+LENPD+2*MSAVE
         GO TO 60
50       IWORK(LMTYPE)=4
         LENRW=40+(IWORK(LMXORD)+4)*NEQ+LENPD
C
C     check lengths of rwork and iwork
60    LENIW=20+NEQ
      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
      IF(TOUT .EQ. T)GO TO 719
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     initialize counters
      IWORK(LNST)=0
      IWORK(LNRE)=0
      IWORK(LNJE)=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.
      CALL XERRWV(
     *49HDASSL--  THE LAST STEP TERMINATED WITH A NEGATIVE,
     *49,201,0,0,0,0,0,0.0D0,0.0D0)
      CALL XERRWV(
     *47HDASSL--  VALUE (=I1) OF IDID AND NO APPROPRIATE,
     *47,202,0,1,IDID,0,0,0.0D0,0.0D0)
      CALL XERRWV(
     *41HDASSL--  ACTION WAS TAKEN. RUN TERMINATED,
     *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.
      LE=LDELTA+NEQ
      LWT=LE+NEQ
      LPHI=LWT+NEQ
      LPD=LPHI+(IWORK(LMXORD)+1)*NEQ
      LWM=LPD
      NPD=1
      NTEMP=NPD+LENPD
      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 ho, 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)/HMAX
         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, if applicable
340   IF (INFO(11) .EQ. 0) GO TO 350
      CALL DDAINI(T,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))
      IF (IDID .LT. 0) GO TO 390
C
C     load h with ho.  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
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
      DONE = .FALSE.
      TN=RWORK(LTN)
      H=RWORK(LH)
      IF(INFO(7) .EQ. 0) GO TO 410
         RH = DABS(H)/HMAX
         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
      IDID=2
      T=TSTOP
      DONE = .TRUE.
      GO TO 490
460   TNEXT=TN+H*(1.0D0+4.0D0*UROUND)
      IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 490
      H=(TSTOP-TN)*(1.0D0-4.0D0*UROUND)
      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
C     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.500)
     *   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
      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))
527   IF(IDID.LT.0)GO TO 600
C
C------------------------------------------------------
C     this block handles the case of a successful
C     return from ddastp (idid=1) test for
C     stop conditions.
C--------------------------------------------------------
C
      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*(1.0D0+4.0D0*UROUND)
      IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 500
      H=(TSTOP-TN)*(1.0D0-4.0D0*UROUND)
      GO TO 500
545   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   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 ddassl are made from
C     this block.
C--------------------------------------------------------
C
590   CONTINUE
      RWORK(LTN)=TN
      RWORK(LH)=H
      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   CALL XERRWV(
     *38HDASSL--  AT CURRENT T (=R1)  500 STEPS,
     *38,610,0,0,0,0,1,TN,0.0D0)
      CALL XERRWV(48HDASSL--  TAKEN ON THIS CALL BEFORE REACHING TOUT,
     *48,611,0,0,0,0,0,0.0D0,0.0D0)
      GO TO 690
C
C     too much accuracy for machine precision
620   CALL XERRWV(
     *47HDASSL--  AT T (=R1) TOO MUCH ACCURACY REQUESTED,
     *47,620,0,0,0,0,1,TN,0.0D0)
      CALL XERRWV(
     *48HDASSL--  FOR PRECISION OF MACHINE. RTOL AND ATOL,
     *48,621,0,0,0,0,0,0.0D0,0.0D0)
      CALL XERRWV(
     *45HDASSL--  WERE INCREASED TO APPROPRIATE VALUES,
     *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   CALL XERRWV(
     *38HDASSL--  AT T (=R1) SOME ELEMENT OF WT,
     *38,630,0,0,0,0,1,TN,0.0D0)
      CALL XERRWV(28HDASSL--  HAS BECOME .LE. 0.0,
     *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   CALL XERRWV(
     *44HDASSL--  AT T (=R1) AND STEPSIZE H (=R2) THE,
     *44,640,0,0,0,0,2,TN,H)
      CALL XERRWV(
     *57HDASSL--  ERROR TEST FAILED REPEATEDLY OR WITH ABS(H)=HMIN,
     *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   CALL XERRWV(
     *44HDASSL--  AT T (=R1) AND STEPSIZE H (=R2) THE,
     *44,650,0,0,0,0,2,TN,H)
      CALL XERRWV(
     *48HDASSL--  CORRECTOR FAILED TO CONVERGE REPEATEDLY,
     *48,651,0,0,0,0,0,0.0D0,0.0D0)
      CALL XERRWV(
     *28HDASSL--  OR WITH ABS(H)=HMIN,
     *28,652,0,0,0,0,0,0.0D0,0.0D0)
      GO TO 690
C
C     the iteration matrix is singular
660   CALL XERRWV(
     *44HDASSL--  AT T (=R1) AND STEPSIZE H (=R2) THE,
     *44,660,0,0,0,0,2,TN,H)
      CALL XERRWV(
     *37HDASSL--  ITERATION MATRIX IS SINGULAR,
     *37,661,0,0,0,0,0,0.0D0,0.0D0)
      GO TO 690
C
C     corrector failure preceeded by error test failures.
670   CALL XERRWV(
     *44HDASSL--  AT T (=R1) AND STEPSIZE H (=R2) THE,
     *44,670,0,0,0,0,2,TN,H)
      CALL XERRWV(
     *49HDASSL--  CORRECTOR COULD NOT CONVERGE.  ALSO, THE,
     *49,671,0,0,0,0,0,0.0D0,0.0D0)
      CALL XERRWV(
     *38HDASSL--  ERROR TEST FAILED REPEATEDLY.,
     *38,672,0,0,0,0,0,0.0D0,0.0D0)
      GO TO 690
C
C     corrector failure because ires = -1
675   CALL XERRWV(
     *44HDASSL--  AT T (=R1) AND STEPSIZE H (=R2) THE,
     *44,675,0,0,0,0,2,TN,H)
      CALL XERRWV(
     *45HDASSL--  CORRECTOR COULD NOT CONVERGE BECAUSE,
     *455,676,0,0,0,0,0,0.0D0,0.0D0)
      CALL XERRWV(
     *36HDASSL--  IRES WAS EQUAL TO MINUS ONE,
     *36,677,0,0,0,0,0,0.0D0,0.0D0)
      GO TO 690
C
C     failure because ires = -2
680   CALL XERRWV(
     *40HDASSL--  AT T (=R1) AND STEPSIZE H (=R2),
     *40,680,0,0,0,0,2,TN,H)
      CALL XERRWV(
     *36HDASSL--  IRES WAS EQUAL TO MINUS TWO,
     *36,681,0,0,0,0,0,0.0D0,0.0D0)
      GO TO 690
C
C     failed to compute initial yprime
685   CALL XERRWV(
     *44HDASSL--  AT T (=R1) AND STEPSIZE H (=R2) THE,
     *44,685,0,0,0,0,2,TN,HO)
      CALL XERRWV(
     *45HDASSL--  INITIAL YPRIME COULD NOT BE COMPUTED,
     *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   CALL XERRWV(
     *55HDASSL--  SOME ELEMENT OF INFO VECTOR IS NOT ZERO OR ONE,
     *55,1,0,0,0,0,0,0.0D0,0.0D0)
      GO TO 750
702   CALL XERRWV(25HDASSL--  NEQ (=I1) .LE. 0,
     *25,2,0,1,NEQ,0,0,0.0D0,0.0D0)
      GO TO 750
703   CALL XERRWV(34HDASSL--  MAXORD (=I1) NOT IN RANGE,
     *34,3,0,1,MXORD,0,0,0.0D0,0.0D0)
      GO TO 750
704   CALL XERRWV(
     *60HDASSL--  RWORK LENGTH NEEDED, LENRW (=I1), EXCEEDS LRW (=I2),
     *60,4,0,2,LENRW,LRW,0,0.0D0,0.0D0)
      GO TO 750
705   CALL XERRWV(
     *60HDASSL--  IWORK LENGTH NEEDED, LENIW (=I1), EXCEEDS LIW (=I2),
     *60,5,0,2,LENIW,LIW,0,0.0D0,0.0D0)
      GO TO 750
706   CALL XERRWV(
     *39HDASSL--  SOME ELEMENT OF RTOL IS .LT. 0,
     *39,6,0,0,0,0,0,0.0D0,0.0D0)
      GO TO 750
707   CALL XERRWV(
     *39HDASSL--  SOME ELEMENT OF ATOL IS .LT. 0,
     *39,7,0,0,0,0,0,0.0D0,0.0D0)
      GO TO 750
708   CALL XERRWV(
     *47HDASSL--  ALL ELEMENTS OF RTOL AND ATOL ARE ZERO,
     *47,8,0,0,0,0,0,0.0D0,0.0D0)
      GO TO 750
709   CALL XERRWV(
     *54HDASSL--  INFO(4) = 1 AND TSTOP (=R1) BEHIND TOUT (=R2),
     *54,9,0,0,0,0,2,TSTOP,TOUT)
      GO TO 750
710   CALL XERRWV(28HDASSL--  HMAX (=R1) .LT. 0.0,
     *28,10,0,0,0,0,1,HMAX,0.0D0)
      GO TO 750
711   CALL XERRWV(34HDASSL--  TOUT (=R1) BEHIND T (=R2),
     *34,11,0,0,0,0,2,TOUT,T)
      GO TO 750
712   CALL XERRWV(29HDASSL--  INFO(8)=1 AND H0=0.0,
     *29,12,0,0,0,0,0,0.0D0,0.0D0)
      GO TO 750
713   CALL XERRWV(39HDASSL--  SOME ELEMENT OF WT IS .LE. 0.0,
     *39,13,0,0,0,0,0,0.0D0,0.0D0)
      GO TO 750
714   CALL XERRWV(
     *61HDASSL--  TOUT (=R1) TOO CLOSE TO T (=R2) TO START INTEGRATION,
     *61,14,0,0,0,0,2,TOUT,T)
      GO TO 750
715   CALL XERRWV(
     *49HDASSL--  INFO(4)=1 AND TSTOP (=R1) BEHIND T (=R2),
     *49,15,0,0,0,0,2,TSTOP,T)
      GO TO 750
717   CALL XERRWV(
     *52HDASSL--  ML (=I1) ILLEGAL. EITHER .LT. 0 OR .GT. NEQ,
     *52,17,0,1,IWORK(LML),0,0,0.0D0,0.0D0)
      GO TO 750
718   CALL XERRWV(
     *52HDASSL--  MU (=I1) ILLEGAL. EITHER .LT. 0 OR .GT. NEQ,
     *52,18,0,1,IWORK(LMU),0,0,0.0D0,0.0D0)
      GO TO 750
719   CALL XERRWV(
     *39HDASSL--  TOUT (=R1) IS EQUAL TO T (=R2),
     *39,19,0,0,0,0,2,TOUT,T)
      GO TO 750
750   IF(INFO(1).EQ.-1) GO TO 760
      INFO(1)=-1
      IDID=-33
      RETURN
760   CALL XERRWV(
     *46HDASSL--  REPEATED OCCURRENCES OF ILLEGAL INPUT,
     *46,801,0,0,0,0,0,0.0D0,0.0D0)
770   CALL XERRWV(
     *47HDASSL--  RUN TERMINATED. APPARENT INFINITE LOOP,
     *47,802,1,0,0,0,0,0.0D0,0.0D0)
      RETURN
C-----------end of subroutine ddassl-------------------------------------
      END
      SUBROUTINE DDASTP(X,Y,YPRIME,NEQ,
     *  RES,JAC,H,WT,JSTART,IDID,RPAR,IPAR,
     *  PHI,DELTA,E,WM,IWM,
     *  ALPHA,BETA,GAMMA,PSI,SIGMA,
     *  CJ,CJOLD,HOLD,S,HMIN,UROUND,
     *  IPHASE,JCALC,K,KOLD,NS,NONNEG)
C
C***BEGIN PROLOGUE  DDASTP
C***REFER TO  DDASSL
C***ROUTINES CALLED  DDANRM,DDAJAC,DDASLV,DDATRP
C***COMMON BLOCKS    DDA001
C***DATE WRITTEN   830315   (YYMMDD)
C***REVISION DATE  830315   (YYMMDD)
C***END PROLOGUE  DDASTP
C
C
C-----------------------------------------------------------------------
C     dastep solves a system of differential/
C     algebraic equations of the form
C     g(x,y,yprime) = 0,  for one step (normally
C     from x to x+h).
C
C     the methods used are modified divided
C     difference,fixed leading coefficient
C     forms of backward differentiation
C     formulas. the code adjusts the stepsize
C     and order to control the local error per
C     step.
C
C
C     the parameters represent
C     x  --        independent variable
C     y  --        solution vector at x
C     yprime --    derivative of solution vector
C                  after successful step
C     neq --       number of equations to be integrated
C     res --       external user-supplied subroutine
C                  to evaluate the residual.  the call is
C                  call res(x,y,yprime,delta,ires,rpar,ipar)
C                  x,y,yprime are input.  delta is output.
C                  on input, ires=0.  res should alter ires only
C                  if it encounters an illegal value of y or a
C                  stop condition.  set ires=-1 if an input value
C                  of y is illegal, and dastep will try to solve
C                  the problem without getting ires = -1.  if
C                  ires=-2, dastep returns control to the calling
C                  program with idid = -11.
C     jac --       external user-supplied routine to evaluate
C                  the iteration matrix (this is optional)
C                  the call is of the form
C                  call jac(x,y,yprime,pd,cj,rpar,ipar)
C                  pd is the matrix of partial derivatives,
C                  pd=dg/dy+cj*dg/dyprime
C     h --         appropriate step size for next step.
C                  normally determined by the code
C     wt --        vector of weights for error criterion.
C     jstart --    integer variable set 0 for
C                  first step, 1 otherwise.
C     idid --      completion code with the following meanings%
C                  idid= 1 -- the step was completed successfully
C                  idid=-6 -- the error test failed repeatedly
C                  idid=-7 -- the corrector could not converge
C                  idid=-8 -- the iteration matrix is singular
C                  idid=-9 -- the corrector could not converge.
C                             there were repeated error test
C                             failures on this step.
C                  idid=-10-- the corrector could not converge
C                             because ires was equal to minus one
C                  idid=-11-- ires equal to -2 was encountered,
C                             and control is being returned to
C                             the calling program
C     rpar,ipar -- real and integer parameter arrays that
C                  are used for communication between the
C                  calling program and external user routines
C                  they are not altered by dastep
C     phi --       array of divided differences used by
C                  dastep. the length is neq*(k+1),where
C                  k is the maximum order
C     delta,e --   work vectors for dastep of length neq
C     wm,iwm --    real and integer arrays storing
C                  matrix information such as the matrix
C                  of partial derivatives,permutation
C                  vector,and various other information.
C
C     the other parameters are information
C     which is needed internally by dastep to
C     continue from step to step.
C
C-----------------------------------------------------------------------
C
C
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL CONVGD
      DIMENSION Y(1),YPRIME(1),WT(1)
      DIMENSION PHI(NEQ,1),DELTA(1),E(1)
      DIMENSION WM(1),IWM(1)
      DIMENSION PSI(1),ALPHA(1),BETA(1),GAMMA(1),SIGMA(1)
      DIMENSION RPAR(1),IPAR(1)
      EXTERNAL RES,JAC
      COMMON/DDA001/NPD,NTEMP,
     *   LML,LMU,LMXORD,LMTYPE,
     *   LNST,LNRE,LNJE,LETF,LCTF,LIPVT
      DATA MAXIT/4/
      DATA XRATE/0.25D0/
C
C
C
C
C
C-----------------------------------------------------------------------
C     block 1.
C     initialize. on the first call,set
C     the order to 1 and initialize
C     other variables.
C-----------------------------------------------------------------------
C
C     initializations for all calls
      IDID=1
      XOLD=X
      NCF=0
      NSF=0
      NEF=0
      IF(JSTART .NE. 0) GO TO 120
C
C     if this is the first step,perform
C     other initializations
      IWM(LETF) = 0
      IWM(LCTF) = 0
      K=1
      KOLD=0
      HOLD=0.0D0
      JSTART=1
      PSI(1)=H
      CJOLD = 1.0D0/H
      CJ = CJOLD
      S = 100.D0
      JCALC = -1
      DELNRM=1.0D0
      IPHASE = 0
      NS=0
120   CONTINUE
C
C
C
C
C
C-----------------------------------------------------------------------
C     block 2
C     compute coefficients of formulas for
C     this step.
C-----------------------------------------------------------------------
200   CONTINUE
      KP1=K+1
      KP2=K+2
      KM1=K-1
      XOLD=X
      IF(H.NE.HOLD.OR.K .NE. KOLD) NS = 0
      NS=MIN0(NS+1,KOLD+2)
      NSP1=NS+1
      IF(KP1 .LT. NS)GO TO 230
C
      BETA(1)=1.0D0
      ALPHA(1)=1.0D0
      TEMP1=H
      GAMMA(1)=0.0D0
      SIGMA(1)=1.0D0
      DO 210 I=2,KP1
         TEMP2=PSI(I-1)
         PSI(I-1)=TEMP1
         BETA(I)=BETA(I-1)*PSI(I-1)/TEMP2
         TEMP1=TEMP2+H
         ALPHA(I)=H/TEMP1
         SIGMA(I)=DFLOAT(I-1)*SIGMA(I-1)*ALPHA(I)
         GAMMA(I)=GAMMA(I-1)+ALPHA(I-1)/H
210      CONTINUE
      PSI(KP1)=TEMP1
230   CONTINUE
C
C     compute alphas, alpha0
      ALPHAS = 0.0D0
      ALPHA0 = 0.0D0
      DO 240 I = 1,K
        ALPHAS = ALPHAS - 1.0D0/DFLOAT(I)
        ALPHA0 = ALPHA0 - ALPHA(I)
240     CONTINUE
C
C     compute leading coefficient cj
      CJLAST = CJ
      CJ = -ALPHAS/H
C
C     compute variable stepsize error coefficient ck
      CK = DABS(ALPHA(KP1) + ALPHAS - ALPHA0)
      CK = DMAX1(CK,ALPHA(KP1))
C
C     decide whether new jacobian is needed
      TEMP1 = (1.0D0 - XRATE)/(1.0D0 + XRATE)
      TEMP2 = 1.0D0/TEMP1
      IF (CJ/CJOLD .LT. TEMP1 .OR. CJ/CJOLD .GT. TEMP2) JCALC = -1
      IF (CJ .NE. CJLAST) S = 100.D0
C
C     change phi to phi star
      IF(KP1 .LT. NSP1) GO TO 280
      DO 270 J=NSP1,KP1
         DO 260 I=1,NEQ
260         PHI(I,J)=BETA(J)*PHI(I,J)
270      CONTINUE
280   CONTINUE
C
C     update time
      X=X+H
C
C
C
C
C
C-----------------------------------------------------------------------
C     block 3
C     predict the solution and derivative,
C     and solve the corrector equation
C-----------------------------------------------------------------------
C
C     first,predict the solution and derivative
300   CONTINUE
      DO 310 I=1,NEQ
         Y(I)=PHI(I,1)
310      YPRIME(I)=0.0D0
      DO 330 J=2,KP1
         DO 320 I=1,NEQ
            Y(I)=Y(I)+PHI(I,J)
320         YPRIME(I)=YPRIME(I)+GAMMA(J)*PHI(I,J)
330   CONTINUE
      PNORM = DDANRM (NEQ,Y,WT,RPAR,IPAR)
C
C
C
C     solve the corrector equation using a
C     modified newton scheme.
      CONVGD= .TRUE.
      M=0
      IWM(LNRE)=IWM(LNRE)+1
      IRES = 0
      CALL RES(X,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
      IF (IRES .LT. 0) GO TO 380
C
C
C     if indicated,reevaluate the
C     iteration matrix pd = dg/dy + cj*dg/dyprime
C     (where g(x,y,yprime)=0). set
C     jcalc to 0 as an indicator that
C     this has been done.
      IF(JCALC .NE. -1)GO TO 340
      IWM(LNJE)=IWM(LNJE)+1
      JCALC=0
      CALL DDAJAC(NEQ,X,Y,YPRIME,DELTA,CJ,H,
     * IER,WT,E,WM,IWM,RES,IRES,UROUND,JAC,RPAR,IPAR)
      CJOLD=CJ
      S = 100.D0
      IF (IRES .LT. 0) GO TO 380
      IF(IER .NE. 0)GO TO 380
      NSF=0
C
C
C     initialize the error accumulation vector e.
340   CONTINUE
      DO 345 I=1,NEQ
345      E(I)=0.0D0
C
      S = 100.E0
C
C
C     corrector loop.
350   CONTINUE
C
C     multiply residual by temp1 to accelerate convergence
      TEMP1 = 2.0D0/(1.0D0 + CJ/CJOLD)
      DO 355 I = 1,NEQ
355     DELTA(I) = DELTA(I) * TEMP1
C
C     compute a new iterate (back-substitution).
C     store the correction in delta.
      CALL DDASLV(NEQ,DELTA,WM,IWM)
C
C     update y,e,and yprime
      DO 360 I=1,NEQ
         Y(I)=Y(I)-DELTA(I)
         E(I)=E(I)-DELTA(I)
360      YPRIME(I)=YPRIME(I)-CJ*DELTA(I)
C
C     test for convergence of the iteration
      DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
      IF (DELNRM .LE. 100.D0*UROUND*PNORM) GO TO 375
      IF (M .GT. 0) GO TO 365
         OLDNRM = DELNRM
         GO TO 367
365   RATE = (DELNRM/OLDNRM)**(1.0D0/DFLOAT(M))
      IF (RATE .GT. 0.90D0) GO TO 370
      S = RATE/(1.0D0 - RATE)
367   IF (S*DELNRM .LE. 0.33D0) GO TO 375
C
C     the corrector has not yet converged.
C     update m and test whether the
C     maximum number of iterations have
C     been tried.
      M=M+1
      IF(M.GE.MAXIT)GO TO 370
C
C     evaluate the residual
C     and go back to do another iteration
      IWM(LNRE)=IWM(LNRE)+1
      IRES = 0
      CALL RES(X,Y,YPRIME,DELTA,IRES,
     *  RPAR,IPAR)
      IF (IRES .LT. 0) GO TO 380
      GO TO 350
C
C
C     the corrector failed to converge in maxit
C     iterations. if the iteration matrix
C     is not current,re-do the step with
C     a new iteration matrix.
370   CONTINUE
      IF(JCALC.EQ.0)GO TO 380
      JCALC=-1
      GO TO 300
C
C
C     the iteration has converged.  if nonnegativity of solution is
C     required, set the solution nonnegative, if the perturbation
C     to do it is small enough.  if the change is too large, then
C     consider the corrector iteration to have failed.
375   IF(NONNEG .EQ. 0) GO TO 390
      DO 377 I = 1,NEQ
377      DELTA(I) = DMIN1(Y(I),0.0D0)
      DELNRM = DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
      IF(DELNRM .GT. 0.33D0) GO TO 380
      DO 378 I = 1,NEQ
378      E(I) = E(I) - DELTA(I)
      GO TO 390
C
C
C     exits from block 3
C     no convergence with current iteration
C     matrix,or singular iteration matrix
380   CONVGD= .FALSE.
390   JCALC = 1
      IF(.NOT.CONVGD)GO TO 600
C
C
C
C
C
C-----------------------------------------------------------------------
C     block 4
C     estimate the errors at orders k,k-1,k-2
C     as if constant stepsize was used. estimate
C     the local error at order k and test
C     whether the current step is successful.
C-----------------------------------------------------------------------
C
C     estimate errors at orders k,k-1,k-2
      ENORM = DDANRM(NEQ,E,WT,RPAR,IPAR)
      ERK = SIGMA(K+1)*ENORM
      TERK = FLOAT(K+1)*ERK
      EST = ERK
      KNEW=K
      IF(K .EQ. 1)GO TO 430
      DO 405 I = 1,NEQ
405     DELTA(I) = PHI(I,KP1) + E(I)
      ERKM1=SIGMA(K)*DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
      TERKM1 = FLOAT(K)*ERKM1
      IF(K .GT. 2)GO TO 410
      IF(TERKM1 .LE. 0.5*TERK)GO TO 420
      GO TO 430
410   CONTINUE
      DO 415 I = 1,NEQ
415     DELTA(I) = PHI(I,K) + DELTA(I)
      ERKM2=SIGMA(K-1)*DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
      TERKM2 = FLOAT(K-1)*ERKM2
      IF(DMAX1(TERKM1,TERKM2).GT.TERK)GO TO 430
C     lower the order
420   CONTINUE
      KNEW=K-1
      EST = ERKM1
C
C
C     calculate the local error for the current step
C     to see if the step was successful
430   CONTINUE
      ERR = CK * ENORM
      IF(ERR .GT. 1.0D0)GO TO 600
C
C
C
C
C
C-----------------------------------------------------------------------
C     block 5
C     the step is successful. determine
C     the best order and stepsize for
C     the next step. update the differences
C     for the next step.
C-----------------------------------------------------------------------
      IDID=1
      IWM(LNST)=IWM(LNST)+1
      KDIFF=K-KOLD
      KOLD=K
      HOLD=H
C
C
C     estimate the error at order k+1 unless%
C        already decided to lower order, or
C        already using maximum order, or
C        stepsize not constant, or
C        order raised in previous step
      IF(KNEW.EQ.KM1.OR.K.EQ.IWM(LMXORD))IPHASE=1
      IF(IPHASE .EQ. 0)GO TO 545
      IF(KNEW.EQ.KM1)GO TO 540
      IF(K.EQ.IWM(LMXORD)) GO TO 550
      IF(KP1.GE.NS.OR.KDIFF.EQ.1)GO TO 550
      DO 510 I=1,NEQ
510      DELTA(I)=E(I)-PHI(I,KP2)
      ERKP1 = (1.0D0/DFLOAT(K+2))*DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
      TERKP1 = FLOAT(K+2)*ERKP1
      IF(K.GT.1)GO TO 520
      IF(TERKP1.GE.0.5D0*TERK)GO TO 550
      GO TO 530
520   IF(TERKM1.LE.DMIN1(TERK,TERKP1))GO TO 540
      IF(TERKP1.GE.TERK.OR.K.EQ.IWM(LMXORD))GO TO 550
C
C     raise order
530   K=KP1
      EST = ERKP1
      GO TO 550
C
C     lower order
540   K=KM1
      EST = ERKM1
      GO TO 550
C
C     if iphase = 0, increase order by one and multiply stepsize by
C     factor two
545   K = KP1
      HNEW = H*2.0D0
      H = HNEW
      GO TO 575
C
C
C     determine the appropriate stepsize for
C     the next step.
550   HNEW=H
      TEMP2=K+1
      R=(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2)
      IF(R .LT. 2.0D0) GO TO 555
      HNEW = 2.0D0*H
      GO TO 560
555   IF(R .GT. 1.0D0) GO TO 560
      R = DMAX1(0.5D0,DMIN1(0.9D0,R))
      HNEW = H*R
560   H=HNEW
C
C
C     update differences for next step
575   CONTINUE
      IF(KOLD.EQ.IWM(LMXORD))GO TO 585
      DO 580 I=1,NEQ
580      PHI(I,KP2)=E(I)
585   CONTINUE
      DO 590 I=1,NEQ
590      PHI(I,KP1)=PHI(I,KP1)+E(I)
      DO 595 J1=2,KP1
         J=KP1-J1+1
         DO 595 I=1,NEQ
595      PHI(I,J)=PHI(I,J)+PHI(I,J+1)
      RETURN
C
C
C
C
C
C-----------------------------------------------------------------------
C     block 6
C     the step is unsuccessful. restore x,psi,phi
C     determine appropriate stepsize for
C     continuing the integration, or exit with
C     an error flag if there have been many
C     failures.
C-----------------------------------------------------------------------
600   IPHASE = 1
C
C     restore x,phi,psi
      X=XOLD
      IF(KP1.LT.NSP1)GO TO 630
      DO 620 J=NSP1,KP1
         TEMP1=1.0D0/BETA(J)
         DO 610 I=1,NEQ
610         PHI(I,J)=TEMP1*PHI(I,J)
620      CONTINUE
630   CONTINUE
      DO 640 I=2,KP1
640      PSI(I-1)=PSI(I)-H
C
C
C     test whether failure is due to corrector iteration
C     or error test
      IF(CONVGD)GO TO 660
      IWM(LCTF)=IWM(LCTF)+1
C
C
C     the newton iteration failed to converge with
C     a current iteration matrix.  determine the cause
C     of the failure and take appropriate action.
      IF(IER.EQ.0)GO TO 650
C
C     the iteration matrix is singular. reduce
C     the stepsize by a factor of 4. if
C     this happens three times in a row on
C     the same step, return with an error flag
      NSF=NSF+1
      R = 0.25D0
      H=H*R
      IF (NSF .LT. 3 .AND. DABS(H) .GE. HMIN) GO TO 690
      IDID=-8
      GO TO 675
C
C
C     the newton iteration failed to converge for a reason
C     other than a singular iteration matrix.  if ires = -2, then
C     return.  otherwise, reduce the stepsize and try again, unless
C     too many failures have occured.
650   CONTINUE
      IF (IRES .GT. -2) GO TO 655
      IDID = -11
      GO TO 675
655   NCF = NCF + 1
      R = 0.25D0
      H = H*R
      IF (NCF .LT. 10 .AND. DABS(H) .GE. HMIN) GO TO 690
      IDID = -7
      IF (IRES .LT. 0) IDID = -10
      IF (NEF .GE. 3) IDID = -9
      GO TO 675
C
C
C     the newton scheme converged,and the cause
C     of the failure was the error estimate
C     exceeding the tolerance.
660   NEF=NEF+1
      IWM(LETF)=IWM(LETF)+1
      IF (NEF .GT. 1) GO TO 665
C
C     on first error test failure, keep current order or lower
C     order by one.  compute new stepsize based on differences
C     of the solution.
      K = KNEW
      TEMP2 = K + 1
      R = 0.90D0*(2.0D0*EST+0.0001D0)**(-1.0D0/TEMP2)
      R = DMAX1(0.25D0,DMIN1(0.9D0,R))
      H = H*R
      IF (DABS(H) .GE. HMIN) GO TO 690
      IDID = -6
      GO TO 675
C
C     on second error test failure, use the current order or
C     decrease order by one.  reduce the stepsize by a factor of
C     one quarter.
665   IF (NEF .GT. 2) GO TO 670
      K = KNEW
      H = 0.25D0*H
      IF (DABS(H) .GE. HMIN) GO TO 690
      IDID = -6
      GO TO 675
C
C     on third and subsequent error test failures, set the order to
C     one and reduce the stepsize by a factor of one quarter
670   K = 1
      H = 0.25D0*H
      IF (DABS(H) .GE. HMIN) GO TO 690
      IDID = -6
      GO TO 675
C
C
C
C
C     for all crashes, restore y to its last value,
C     interpolate to find yprime at last x, and return
675   CONTINUE
      CALL DDATRP(X,X,Y,YPRIME,NEQ,K,PHI,PSI)
      RETURN
C
C
C     go back and try this step again
690   GO TO 200
C
C------end of subroutine dastep------
      END
      SUBROUTINE DDAINI(X,Y,YPRIME,NEQ,
     *   RES,JAC,H,WT,IDID,RPAR,IPAR,
     *   PHI,DELTA,E,WM,IWM,
     *   HMIN,UROUND,NONNEG)
C
C***BEGIN PROLOGUE  DDAINI
C***REFER TO  DDASSL
C***ROUTINES CALLED  DDANRM,DDAJAC,DDASLV
C***COMMON BLOCKS    DDA001
C***DATE WRITTEN   830315   (YYMMDD)
C***REVISION DATE  830315   (YYMMDD)
C***END PROLOGUE DDAINI
C
C-------------------------------------------------------
C     ddaini takes one step of size h or smaller
C     with the backward euler method, to
C     find yprime at the initial time x. a modified
C     damped newton iteration is used to
C     solve the corrector iteration.
C
C     the initial guess yprime is used in the
C     prediction, and in forming the iteration
C     matrix, but is not involved in the
C     error test. this may have trouble
C     converging if the initial guess is no
C     good, or if g(xy,yprime) depends
C     nonlinearly on yprime.
C
C     the parameters represent:
C     x --         independent variable
C     y --         solution vector at x
C     yprime --    derivative of solution vector
C     neq --       number of equations
C     h --         stepsize. imder may use a stepsize
C                  smaller than h.
C     wt --        vector of weights for error
C                  criterion
C     idid --      completion code with the following meanings
C                  idid= 1 -- yprime was found successfully
C                  idid=-12 -- ddaini failed to find yprime
C     rpar,ipar -- real and integer parameter arrays
C                  that are not altered by ddaini
C     phi --       work space for ddaini
C     delta,e --   work space for ddaini
C     wm,iwm --    real and integer arrays storing
C                  matrix information
C
C-----------------------------------------------------------------
C
C

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL CONVGD
      DIMENSION Y(1),YPRIME(1),WT(1)
      DIMENSION PHI(NEQ,1),DELTA(1),E(1)
      DIMENSION WM(1),IWM(1)
      DIMENSION RPAR(1),IPAR(1)
      EXTERNAL RES,JAC
      COMMON/DDA001/NPD,NTEMP,
     *  LML,LMU,LMXORD,LMTYPE,
     *  LNST,LNRE,LNJE,LETF,LCTF,LIPVT

      DATA MAXIT/10/,MJAC/5/
      DATA DAMP/0.75D0/

C
C
C---------------------------------------------------
C     block 1.
C     initializations.
C---------------------------------------------------
C
      IDID=1
      NEF=0
      NCF=0
      NSF=0
      YNORM=DDANRM(NEQ,Y,WT,RPAR,IPAR)
C
C     save y and yprime in phi
      DO 100 I=1,NEQ
         PHI(I,1)=Y(I)
100      PHI(I,2)=YPRIME(I)

C
C
C----------------------------------------------------
C     block 2.
C     do one backward euler step.
C----------------------------------------------------
C
C     set up for start of corrector iteration
200   CJ=1.0D0/H
      XNEW=X+H
C
C     predict solution and derivative

      DO 250 I=1,NEQ
250     Y(I)=Y(I)+H*YPRIME(I)
C
      JCALC=-1
      M=0
      CONVGD=.TRUE.
C
C
C     corrector loop.
300   IWM(LNRE)=IWM(LNRE)+1
      IRES=0

      CALL RES(XNEW,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
      IF (IRES.LT.0) GO TO 430
C
C
C     evaluate the iteration matrix
      IF (JCALC.NE.-1) GO TO 310
      IWM(LNJE)=IWM(LNJE)+1
      JCALC=0
      CALL DDAJAC(NEQ,XNEW,Y,YPRIME,DELTA,CJ,H,
     *   IER,WT,E,WM,IWM,RES,IRES,
     *   UROUND,JAC,RPAR,IPAR)

      S=1000000.D0
      IF (IRES.LT.0) GO TO 430
      IF (IER.NE.0) GO TO 430
      NSF=0

C
C
C
C     multiply residual by damping factor
310   CONTINUE
      DO 320 I=1,NEQ
320      DELTA(I)=DELTA(I)*DAMP

C
C     compute a new iterate (back substitution)
C     store the correction in delta

      CALL DDASLV(NEQ,DELTA,WM,IWM)

C
C     update y and yprime

      DO 330 I=1,NEQ
         Y(I)=Y(I)-DELTA(I)
330      YPRIME(I)=YPRIME(I)-CJ*DELTA(I)

C
C     test for convergence of the iteration.

      DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
      IF (DELNRM.LE.100.D0*UROUND*YNORM)
     *   GO TO 400

      IF (M.GT.0) GO TO 340
         OLDNRM=DELNRM
         GO TO 350

340   RATE=(DELNRM/OLDNRM)**(1.0D0/DFLOAT(M))
      IF (RATE.GT.0.90D0) GO TO 430
      S=RATE/(1.0D0-RATE)

350   IF (S*DELNRM .LE. 0.33D0) GO TO 400
C
C
C     the corrector has not yet converged. update
C     m and and test whether the maximum
C     number of iterations have been tried.
C     every mjac iterations, get a new
C     iteration matrix.

      M=M+1
      IF (M.GE.MAXIT) GO TO 430

      IF ((M/MJAC)*MJAC.EQ.M) JCALC=-1

      GO TO 300

C
C
C     the iteration has converged.
C     check nonnegativity constraints
400   IF (NONNEG.EQ.0) GO TO 450
      DO 410 I=1,NEQ
410      DELTA(I)=DMIN1(Y(I),0.0D0)

      DELNRM=DDANRM(NEQ,DELTA,WT,RPAR,IPAR)
      IF (DELNRM.GT.0.33D0) GO TO 430

      DO 420 I=1,NEQ
         Y(I)=Y(I)-DELTA(I)
420      YPRIME(I)=YPRIME(I)-CJ*DELTA(I)
      GO TO 450
C
C
C     exits from corrector loop.
430   CONVGD=.FALSE.
450   IF (.NOT.CONVGD) GO TO 600
C
C
C
C-----------------------------------------------------
C     block 3.
C     the corrector iteration converged.
C     do error test.
C-----------------------------------------------------
C

      DO 510 I=1,NEQ
510      E(I)=Y(I)-PHI(I,1)

      ERR=DDANRM(NEQ,E,WT,RPAR,IPAR)

      IF (ERR.LE.1.0D0) RETURN

C
C
C
C--------------------------------------------------------
C     block 4.
C     the backward euler step failed. restore y
C     and yprime to their original values.
C     reduce stepsize and try again, if
C     possible.
C---------------------------------------------------------
C

600   CONTINUE
      DO 610 I=1,NEQ
         Y(I)=PHI(I,1)
610      YPRIME(I)=PHI(I,2)

      IF (CONVGD) GO TO 640
      IF (IER.EQ.0) GO TO 620
         NSF=NSF+1
         H=H*0.25D0
         IF (NSF.LT.3.AND.DABS(H).GE.HMIN) GO TO 690
         IDID=-12
         RETURN
620   IF (IRES.GT.-2) GO TO 630
         IDID=-12
         RETURN
630   NCF=NCF+1
      H=H*0.25D0
      IF (NCF.LT.10.AND.DABS(H).GE.HMIN) GO TO 690
         IDID=-12
         RETURN

640   NEF=NEF+1
      R=0.90D0/(2.0D0*ERR+0.0001D0)
      R=DMAX1(0.1D0,DMIN1(0.5D0,R))
      H=H*R
      IF (DABS(H).GE.HMIN.AND.NEF.LT.10) GO TO 690
         IDID=-12
         RETURN
690      GO TO 200

C-------------end of subroutine ddaini----------------------
      END
      DOUBLE PRECISION FUNCTION DDANRM(NEQ,V,WT,RPAR,IPAR)
C
C***BEGIN PROLOGUE  DDANRM
C***REFER TO  DDASSL
C***ROUTINES CALLED  (NONE)
C***DATE WRITTEN   830315   (YYMMDD)
C***REVISION DATE  830315   (YYMMDD)
C***END PROLOGUE  DDANRM
C-----------------------------------------------------------------------
C     this function routine computes the weighted
C     root-mean-square norm of the vector of length
C     neq contained in the array v,with weights
C     contained in the array wt of length neq.
C        ddanrm=sqrt((1/neq)*sum(v(i)/wt(i))**2)
C-----------------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION V(NEQ),WT(NEQ)
      DIMENSION RPAR(1),IPAR(1)
      DDANRM = 0.0D0
      VMAX = 0.0D0
      DO 10 I = 1,NEQ
10      IF(DABS(V(I)/WT(I)) .GT. VMAX) VMAX = DABS(V(I)/WT(I))
      IF(VMAX .LE. 0.0D0) GO TO 30
      SUM = 0.0D0
      DO 20 I = 1,NEQ
20      SUM = SUM + ((V(I)/WT(I))/VMAX)**2
      DDANRM = VMAX*DSQRT(SUM/DFLOAT(NEQ))
30    CONTINUE
      RETURN
C------end of function ddanrm------
      END
      SUBROUTINE DDAWTS(NEQ,IWT,RTOL,ATOL,Y,WT,RPAR,IPAR)
C
C***BEGIN PROLOGUE  DDAWTS
C***REFER TO  DDASSL
C***ROUTINES CALLED  (NONE)
C***DATE WRITTEN   830315   (YYMMDD)
C***REVISION DATE  830315   (YYMMDD)
C***END PROLOGUE  DDAWTS
C-----------------------------------------------------------------------
C     this subroutine sets the error weight vector
C     wt according to wt(i)=rtol(i)*abs(y(i))+atol(i),
C     i=1,-,n.
C     rtol and atol are scalars if iwt = 0,
C     and vectors if iwt = 1.
C-----------------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION RTOL(1),ATOL(1),Y(1),WT(1)
      DIMENSION RPAR(1),IPAR(1)
      RTOLI=RTOL(1)
      ATOLI=ATOL(1)
      DO 20 I=1,NEQ
         IF (IWT .EQ.0) GO TO 10
           RTOLI=RTOL(I)
           ATOLI=ATOL(I)
10         WT(I)=RTOLI*DABS(Y(I))+ATOLI
20         CONTINUE
      RETURN
C-----------end of subroutine ddawts-------------------------------------
      END
      SUBROUTINE DDATRP(X,XOUT,YOUT,YPOUT,NEQ,KOLD,PHI,PSI)
C
C***BEGIN PROLOGUE  DDATRP
C***REFER TO  DDASSL
C***ROUTINES CALLED  (NONE)
C***DATE WRITTEN   830315   (YYMMDD)
C***REVISION DATE  830315   (YYMMDD)
C***END PROLOGUE  DDATRP
C
C-----------------------------------------------------------------------
C     the methods in subroutine dastep use polynomials
C     to approximate the solution. ddatrp approximates the
C     solution and its derivative at time xout by evaluating
C     one of these polynomials,and its derivative,there.
C     information defining this polynomial is passed from
C     dastep, so ddatrp cannot be used alone.
C
C     the parameters are%
C     x     the current time in the integration.
C     xout  the time at which the solution is desired
C     yout  the interpolated approximation to y at xout
C           (this is output)
C     ypout the interpolated approximation to yprime at xout
C           (this is output)
C     neq   number of equations
C     kold  order used on last successful step
C     phi   array of scaled divided differences of y
C     psi   array of past stepsize history
C-----------------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION YOUT(1),YPOUT(1)
      DIMENSION PHI(NEQ,1),PSI(1)
      KOLDP1=KOLD+1
      TEMP1=XOUT-X
      DO 10 I=1,NEQ
         YOUT(I)=PHI(I,1)
10       YPOUT(I)=0.0D0
      C=1.0D0
      D=0.0D0
      GAMMA=TEMP1/PSI(1)
      DO 30 J=2,KOLDP1
         D=D*GAMMA+C/PSI(J-1)
         C=C*GAMMA
         GAMMA=(TEMP1+PSI(J-1))/PSI(J)
         DO 20 I=1,NEQ
            YOUT(I)=YOUT(I)+C*PHI(I,J)
20          YPOUT(I)=YPOUT(I)+D*PHI(I,J)
30       CONTINUE
      RETURN
C
C------end of subroutine ddatrp------
      END
      SUBROUTINE DDAJAC(NEQ,X,Y,YPRIME,DELTA,CJ,H,
     *  IER,WT,E,WM,IWM,RES,IRES,UROUND,JAC,RPAR,IPAR)
C
C***BEGIN PROLOGUE  DDAJAC
C***REFER TO  DDASSL
C***ROUTINES CALLED  DGEFA,DGBFA
C***COMMON BLOCKS    DDA001
C***DATE WRITTEN   830315   (YYMMDD)
C***REVISION DATE  830315   (YYMMDD)
C***END PROLOGUE  DDAJAC
C-----------------------------------------------------------------------
C     this routine computes the iteration matrix
C     pd=dg/dy+cj*dg/dyprime (where g(x,y,yprime)=0).
C     here pd is computed by the user-supplied
C     routine jac if iwm(mtype) is 1 or 4, and
C     it is computed by numerical finite differencing
C     if iwm(mtype)is 2 or 5
C     the parameters have the following meanings.
C     y        = array containing predicted values
C     yprime   = array containing predicted derivatives
C     delta    = residual evaluated at (x,y,yprime)
C                (used only if iwm(mtype)=2 or 5)
C     cj       = scalar parameter defining iteration matrix
C     h        = current stepsize in integration
C     ier      = variable which is .ne. 0
C                if iteration matrix is singular,
C                and 0 otherwise.
C     wt       = vector of weights for computing norms
C     e        = work space (temporary) of length neq
C     wm       = real work space for matrices. on
C                output it contains the lu decomposition
C                of the iteration matrix.
C     iwm      = integer work space containing
C                matrix information
C     res      = name of the external user-supplied routine
C                to evaluate the residual function g(x,y,yprime)
C     ires     = flag which is equal to zero if no illegal values
C                in res, and less than zero otherwise.  (if ires
C                is less than zero, the matrix was not completed)
C                in this case (if ires .lt. 0), then ier = 0.
C     uround   = the unit roundoff error of the machine being used.
C     jac      = name of the external user-supplied routine
C                to evaluate the iteration matrix (this routine
C                is only used if iwm(mtype) is 1 or 4)
C-----------------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      EXTERNAL RES,JAC
      DIMENSION Y(1),YPRIME(1),DELTA(1),WT(1),E(1)
      DIMENSION WM(1),IWM(1),RPAR(1),IPAR(1)
      COMMON/DDA001/NPD,NTEMP,
     *  LML,LMU,LMXORD,LMTYPE,
     *  LNST,LNRE,LNJE,LETF,LCTF,LIPVT
C
      IER = 0
      NPDM1=NPD-1
      MTYPE=IWM(LMTYPE)
      GO TO (100,200,300,400,500),MTYPE
C
C
C     dense user-supplied matrix
100   LENPD=NEQ*NEQ
      DO 110 I=1,LENPD
110      WM(NPDM1+I)=0.0D0
      CALL JAC(X,Y,YPRIME,WM(NPD),CJ,RPAR,IPAR)
      GO TO 230
C
C
C     dense finite-difference-generated matrix
200   IRES=0
      NROW=NPDM1
      SQUR = DSQRT(UROUND)
      DO 210 I=1,NEQ
         DEL=SQUR*DMAX1(DABS(Y(I)),DABS(H*YPRIME(I)),
     *     DABS(WT(I)))
         DEL=DSIGN(DEL,H*YPRIME(I))
         DEL=(Y(I)+DEL)-Y(I)
         YSAVE=Y(I)
         YPSAVE=YPRIME(I)
         Y(I)=Y(I)+DEL
         YPRIME(I)=YPRIME(I)+CJ*DEL
         CALL RES(X,Y,YPRIME,E,IRES,RPAR,IPAR)
         IF (IRES .LT. 0) RETURN
         DELINV=1.0D0/DEL
         DO 220 L=1,NEQ
220      WM(NROW+L)=(E(L)-DELTA(L))*DELINV
      NROW=NROW+NEQ
      Y(I)=YSAVE
      YPRIME(I)=YPSAVE
210   CONTINUE
C
C
C     do dense-matrix lu decomposition on pd
230      CALL DGEFA(WM(NPD),NEQ,NEQ,IWM(LIPVT),IER)
      RETURN
C
C
C     dummy section for iwm(mtype)=3
300   RETURN
C
C
C     banded user-supplied matrix
400   LENPD=(2*IWM(LML)+IWM(LMU)+1)*NEQ
      DO 410 I=1,LENPD
410      WM(NPDM1+I)=0.0D0
      CALL JAC(X,Y,YPRIME,WM(NPD),CJ,RPAR,IPAR)
      MEBAND=2*IWM(LML)+IWM(LMU)+1
      GO TO 550
C
C
C     banded finite-difference-generated matrix
500   MBAND=IWM(LML)+IWM(LMU)+1
      MBA=MIN0(MBAND,NEQ)
      MEBAND=MBAND+IWM(LML)
      MEB1=MEBAND-1
      MSAVE=(NEQ/MBAND)+1
      ISAVE=NTEMP-1
      IPSAVE=ISAVE+MSAVE
      IRES=0
      SQUR=DSQRT(UROUND)
      DO 540 J=1,MBA
         DO 510 N=J,NEQ,MBAND
          K= (N-J)/MBAND + 1
          WM(ISAVE+K)=Y(N)
          WM(IPSAVE+K)=YPRIME(N)
          DEL=SQUR*DMAX1(DABS(Y(N)),DABS(H*YPRIME(N)),
     *      DABS(WT(N)))
          DEL=DSIGN(DEL,H*YPRIME(N))
          DEL=(Y(N)+DEL)-Y(N)
          Y(N)=Y(N)+DEL
510       YPRIME(N)=YPRIME(N)+CJ*DEL
      CALL RES(X,Y,YPRIME,E,IRES,RPAR,IPAR)
      IF (IRES .LT. 0) RETURN
      DO 530 N=J,NEQ,MBAND
          K= (N-J)/MBAND + 1
          Y(N)=WM(ISAVE+K)
          YPRIME(N)=WM(IPSAVE+K)
          DEL=SQUR*DMAX1(DABS(Y(N)),DABS(H*YPRIME(N)),
     *      DABS(WT(N)))
          DEL=DSIGN(DEL,H*YPRIME(N))
          DEL=(Y(N)+DEL)-Y(N)
          DELINV=1.0D0/DEL
          I1=MAX0(1,(N-IWM(LMU)))
          I2=MIN0(NEQ,(N+IWM(LML)))
          II=N*MEB1-IWM(LML)+NPDM1
          DO 520 I=I1,I2
520         WM(II+I)=(E(I)-DELTA(I))*DELINV
530      CONTINUE
540   CONTINUE
C
C
C     do lu decomposition of banded pd
550   CALL DGBFA(WM(NPD),MEBAND,NEQ,
     *    IWM(LML),IWM(LMU),IWM(LIPVT),IER)
      RETURN
C------end of subroutine ddajac------
      END
      SUBROUTINE DDASLV(NEQ,DELTA,WM,IWM)
C
C***BEGIN PROLOGUE  DDASLV
C***REFER TO  DDASSL
C***ROUTINES CALLED DGESL,DGBSL
C***COMMON BLOCKS    DDA001
C***DATE WRITTEN   830315   (YYMMDD)
C***REVISION DATE  830315   (YYMMDD)
C***END PROLOGUE  DDASLV
C-----------------------------------------------------------------------
C     this routine manages the solution of the linear
C     system arising in the newton iteration.
C     matrices and real temporary storage and
C     real information are stored in the array wm.
C     integer matrix information is stored in
C     the array iwm.
C     for a dense matrix, the linpack routine
C     dgesl is called.
C     for a banded matrix,the linpack routine
C     dgbsl is called.
C-----------------------------------------------------------------------
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION DELTA(1),WM(1),IWM(1)
      COMMON/DDA001/NPD,NTEMP,LML,LMU,
     *  LMXORD,LMTYPE,
     *  LNST,LNRE,LNJE,LETF,LCTF,LIPVT
C
      MTYPE=IWM(LMTYPE)
      GO TO(100,100,300,400,400),MTYPE
C
C     dense matrix
100   CALL DGESL(WM(NPD),NEQ,NEQ,IWM(LIPVT),DELTA,0)
      RETURN
C
C     dummy section for mtype=3
300   CONTINUE
      RETURN
C
C     banded matrix
400   MEBAND=2*IWM(LML)+IWM(LMU)+1
      CALL DGBSL(WM(NPD),MEBAND,NEQ,IWM(LML),
     *  IWM(LMU),IWM(LIPVT),DELTA,0)
      RETURN
C------end of subroutine ddaslv------
      END
      SUBROUTINE XERRWV (MSG, NMES, NERR, IERT, NI, I1, I2, NR, R1, R2)
      INTEGER MSG, NMES, NERR, IERT, NI, I1, I2, NR,
     1   I, LUN, LUNIT, MESFLG, NWDS
      REAL R1, R2
      DIMENSION MSG(NMES)
C-----------------------------------------------------------------------
C subroutine xerrwv, as given here, constitutes
C a simplified version of the slatec error handling package.
C written by a. c. hindmarsh at lll.  version of january 23, 1980.
C modified by l. r. petzold, april 1982.
C this version is in single precision.
C
C all arguments are input arguments.
C
C msg    = the message (hollerith litteral or integer array).
C nmes   = the length of msg (number of characters).
C nerr   = the error number (not used).
C iert   = the error type..
C          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 ni.
C
C note..  this routine is machine-dependent and specialized for use
C in limited context, in the following ways..
C 1. the number of hollerith characters stored per word, denoted
C    by ncpw below, is set in a data statement below.
C 2. the value of nmes is assumed to be at most 60.
C    (multi-line messages are generated by repeated calls.)
C 3. if iert = 2, control passes to the statement   stop
C    to abort the run.  this statement may be machine-dependent.
C 4. r1 and r2 are assumed to be in real and are printed
C    in d21.13 format.
C 5. the data statement below contains default values of
C       mesflg = print control flag..
C                1 means print all messages (the default).
C                0 means no printing.
C       lunit  = logical unit number for messages.
C                the default is 6 (machine-dependent).
C                to change lunit, change the data statement
C                below.
C-----------------------------------------------------------------------
C the following are instructions for installing this routine
C in different machine environments.
C
C to change the default output unit, change the data statement
C below.
C
C for a different number of characters per word, change the
C data statement setting ncpw below.
C alternatives for various computers are shown in comment
C cards.
C
C for a different run-abort command, change the statement following
C statement 100 at the end.
C-----------------------------------------------------------------------
C the following value of ncpw is valid for the cdc-6600 and
C cdc-7600 computers.
C     data ncpw/10/
C the following is valid for the cray-1 computer.
C     data ncpw/8/
C the following is valid for the burroughs 6700 and 7800 computers.
C     data ncpw/6/
C the following is valid for the pdp-10 computer.
C     data ncpw/5/
C the following is valid for the vax computer with 4 bytes per integer,
C and for the ibm-360, ibm-303x, and ibm-43xx computers.
C     data ncpw/4/
C the following is valid for the pdp-11, or vax with 2-byte integers.
C     data ncpw/2/
C----------------------------------------------------------------------
      DIMENSION NFORM(13)
      DATA NFORM(1)/1H(/,NFORM(2)/1H1/,NFORM(3)/1HX/,NFORM(4)/1H,/,
     1  NFORM(7)/1HA/,NFORM(10)/1H,/,NFORM(11)/1HA/,NFORM(13)/1H)/
      DATA NCPW/4/
      DATA MESFLG/1/,LUNIT/6/
C
      IF (MESFLG .EQ. 0) GO TO 100
C get logical unit number. ---------------------------------------------
      LUN = LUNIT
C get number of words in message. --------------------------------------
      NCH = MIN0(NMES,60)
      NWDS = NCH/NCPW
      CALL S88FMT(2,NWDS,NFORM(5))
      CALL S88FMT(2,NCPW,NFORM(8))
      NREM = NCH - NWDS*NCPW
      IF (NREM .GT. 0) NWDS = NWDS + 1
      IF (NREM .LT. 1) NREM = 1
      CALL S88FMT(1,NREM,NFORM(12))
      WRITE(LUN,NFORM) (MSG(I),I=1,NWDS)
      IF (NI .EQ. 1) WRITE (LUN, 20) I1
 20   FORMAT(6X,23HIN ABOVE MESSAGE,  I1 =,I10)
      IF (NI .EQ. 2) WRITE (LUN, 30) I1,I2
 30   FORMAT(6X,23HIN ABOVE MESSAGE,  I1 =,I10,3X,4HI2 =,I10)
      IF (NR .EQ. 1) WRITE (LUN, 40) R1
 40   FORMAT(6X,23HIN ABOVE MESSAGE,  R1 =,D21.13)
      IF (NR .EQ. 2) WRITE (LUN, 50) R1,R2
 50   FORMAT(6X,15HIN ABOVE,  R1 =,D21.13,3X,4HR2 =,D21.13)
C abort the run if iert = 2. -------------------------------------------
 100  IF (IERT .NE. 2) RETURN
      STOP
C----------------------- end of subroutine xerrwv ----------------------
      END
      SUBROUTINE S88FMT(N,IVALUE,IFMT)
C***begin prologue  s88fmt
C***refer to  xerror
C     abstract
C        s88fmt replaces ifmt(1), ... ,ifmt(n) with the
C        characters corresponding to the n least significant
C        digits of ivalue.
C
C     taken from the bell laboratories port library error handler
C     latest revision ---  7 june 1978
C
C***references
C   jones r.e., *slatec common mathematical library error handling
C    package*, sand78-1189, sandia laboratories, 1978.
C***routines called  (none)
C***end prologue  s88fmt
C
      DIMENSION IFMT(N),IDIGIT(10)
      DATA IDIGIT(1),IDIGIT(2),IDIGIT(3),IDIGIT(4),IDIGIT(5),
     1     IDIGIT(6),IDIGIT(7),IDIGIT(8),IDIGIT(9),IDIGIT(10)
     2     /1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9/
C***first executable statement  s88fmt
      NT = N
      IT = IVALUE
   10    IF (NT .EQ. 0) RETURN
         INDEX = MOD(IT,10)
         IFMT(NT) = IDIGIT(INDEX+1)
         IT = IT/10
         NT = NT - 1
         GO TO 10
      END