text: !
  !    SLATEC: public domain
  !
  *DECK DDEABM
        SUBROUTINE DDEABM (DF, NEQ, T, Y, TOUT, INFO, RTOL, ATOL, IDID,
       +   RWORK, LRW, IWORK, LIW, RPAR, IPAR)
  C***BEGIN PROLOGUE  DDEABM
  C***PURPOSE  Solve an initial value problem in ordinary differential
  C            equations using an Adams-Bashforth method.
  C***LIBRARY   SLATEC (DEPAC)
  C***CATEGORY  I1A1B
  C***TYPE      DOUBLE PRECISION (DEABM-S, DDEABM-D)
  C***KEYWORDS  ADAMS-BASHFORTH METHOD, DEPAC, INITIAL VALUE PROBLEMS,
  C             ODE, ORDINARY DIFFERENTIAL EQUATIONS, PREDICTOR-CORRECTOR
  C***AUTHOR  Shampine, L. F., (SNLA)
  C           Watts, H. A., (SNLA)
  C***DESCRIPTION
  C
  C   This is the Adams code in the package of differential equation
  C   solvers DEPAC, consisting of the codes DDERKF, DDEABM, and DDEBDF.
  C   Design of the package was by L. F. Shampine and H. A. Watts.
  C   It is documented in
  C        SAND79-2374 , DEPAC - Design of a User Oriented Package of ODE
  C                              Solvers.
  C   DDEABM is a driver for a modification of the code ODE written by
  C             L. F. Shampine and M. K. Gordon
  C             Sandia Laboratories
  C             Albuquerque, New Mexico 87185
  C
  C **********************************************************************
  C * ABSTRACT *
  C ************
  C
  C   Subroutine DDEABM uses the Adams-Bashforth-Moulton
  C   Predictor-Corrector formulas of orders one through twelve to
  C   integrate a system of NEQ first order ordinary differential
  C   equations of the form
  C                         DU/DX = DF(X,U)
  C   when the vector Y(*) of initial values for U(*) at X=T is given.
  C   The subroutine integrates from T to TOUT. It is easy to continue the
  C   integration to get results at additional TOUT.  This is the interval
  C   mode of operation.  It is also easy for the routine to return with
  C   the solution at each intermediate step on the way to TOUT.  This is
  C   the intermediate-output mode of operation.
  C
  C   DDEABM uses subprograms DDES, DSTEPS, DINTP, DHSTRT, DHVNRM,
  C   D1MACH, and the error handling routine XERMSG.  The only machine
  C   dependent parameters to be assigned appear in D1MACH.
  C
  C **********************************************************************
  C * Description of The Arguments To DDEABM (An Overview) *
  C **********************************************************************
  C
  C   The Parameters are
  C
  C      DF -- This is the name of a subroutine which you provide to
  C             define the differential equations.
  C
  C      NEQ -- This is the number of (first order) differential
  C             equations to be integrated.
  C
  C      T -- This is a DOUBLE PRECISION value of the independent
  C           variable.
  C
  C      Y(*) -- This DOUBLE PRECISION array contains the solution
  C              components at T.
  C
  C      TOUT -- This is a DOUBLE PRECISION point at which a solution is
  C              desired.
  C
  C      INFO(*) -- The basic task of the code is to integrate the
  C             differential equations from T to TOUT and return an
  C             answer at TOUT.  INFO(*) is an INTEGER array which is used
  C             to communicate exactly how you want this task to be
  C             carried out.
  C
  C      RTOL, ATOL -- These DOUBLE PRECISION quantities represent
  C                    relative and absolute error tolerances which you
  C                    provide to indicate how accurately you wish the
  C                    solution to be computed.  You may choose them to be
  C                    both scalars or else both vectors.
  C
  C      IDID -- This scalar quantity is an indicator reporting what
  C             the code did.  You must monitor this INTEGER variable to
  C             decide what action to take next.
  C
  C      RWORK(*), LRW -- RWORK(*) is a DOUBLE PRECISION work array of
  C             length LRW which provides the code with needed storage
  C             space.
  C
  C      IWORK(*), LIW -- IWORK(*) is an INTEGER work array of length LIW
  C             which provides the code with needed storage space and an
  C             across call flag.
  C
  C      RPAR, IPAR -- These are DOUBLE PRECISION and INTEGER parameter
  C             arrays which you can use for communication between your
  C             calling program and the DF subroutine.
  C
  C  Quantities which are used as input items are
  C             NEQ, T, Y(*), TOUT, INFO(*),
  C             RTOL, ATOL, RWORK(1), LRW and LIW.
  C
  C  Quantities which may be altered by the code are
  C             T, Y(*), INFO(1), RTOL, ATOL,
  C             IDID, RWORK(*) and IWORK(*).
  C
  C **********************************************************************
  C * INPUT -- What To Do On The First Call To DDEABM *
  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      DF -- Provide a subroutine of the form
  C                               DF(X,U,UPRIME,RPAR,IPAR)
  C             to define the system of first order differential equations
  C             which is to be solved.  For the given values of X and the
  C             vector  U(*)=(U(1),U(2),...,U(NEQ)) , the subroutine must
  C             evaluate the NEQ components of the system of differential
  C             equations  DU/DX=DF(X,U)  and store the derivatives in the
  C             array UPRIME(*), that is,  UPRIME(I) = * DU(I)/DX *  for
  C             equations I=1,...,NEQ.
  C
  C             Subroutine DF must NOT alter X or U(*).  You must declare
  C             the name df in an external statement in your program that
  C             calls DDEABM.  You must dimension U and UPRIME in DF.
  C
  C             RPAR and IPAR are DOUBLE PRECISION and INTEGER parameter
  C             arrays which you can use for communication between your
  C             calling program and subroutine DF. They are not used or
  C             altered by DDEABM.  If you do not need RPAR or IPAR,
  C             ignore these parameters by treating them as dummy
  C             arguments. If you do choose to use them, dimension them in
  C             your calling program and in DF as arrays of appropriate
  C             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             You must use a program variable for T because the code
  C             changes its value.
  C
  C      Y(*) -- Set this vector to the initial values of the NEQ solution
  C             components at the initial point.  You must dimension Y at
  C             least NEQ in your calling program.
  C
  C      TOUT -- Set it to the first point at which a solution
  C             is desired.  You can take TOUT = T, in which case the code
  C             will evaluate the derivative of the solution at T and
  C             return. 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 following
  C             each intermediate step (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 permissible 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 to accommodate other members of
  C             DEPAC or possible future extensions, though DDEABM uses
  C             only the first four 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      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             program variables because the code may change them.  You
  C             have two choices --
  C                  Both RTOL and ATOL are scalars. (INFO(2)=0)
  C                  Both RTOL and ATOL are vectors. (INFO(2)=1)
  C             In either case all components must be non-negative.
  C
  C             The tolerances are used by the code in a local error 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 Euclidean norm is used to measure
  C             the size of vectors, and the error test uses the magnitude
  C             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.  Roughly speaking, they produce a solution Y(T)
  C             which satisfies the differential equations with a
  C             residual R(T),    DY(T)/DT = DF(T,Y(T)) + R(T)   ,
  C             and, almost always, R(T) is bounded by the error
  C             tolerances.  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.D0 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             Proper selection of the absolute error control parameters
  C             ATOL  requires you to have some idea of the scale of the
  C             solution components.  To acquire this information may mean
  C             that you will have to solve the problem more than once. In
  C             the absence of scale information, you should ask for some
  C             relative accuracy in all the components (by setting  RTOL
  C             values non-zero) and perhaps impose extremely small
  C             absolute error tolerances to protect against the danger of
  C             a solution component becoming zero.
  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 DOUBLE PRECISION work array of length
  C             LRW in your calling program.
  C
  C      RWORK(1) -- If you have set INFO(4)=0, you can ignore this
  C             optional input parameter.  Otherwise you must define a
  C             stopping point TSTOP by setting   RWORK(1) = TSTOP.
  C             (for some problems it may not be permissible 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.)
  C
  C      LRW -- Set it to the declared length of the RWORK array.
  C             You must have  LRW .GE. 130+21*NEQ
  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. 51
  C
  C      RPAR, IPAR -- These are parameter arrays, of DOUBLE PRECISION and
  C             INTEGER type, respectively.  You can use them for
  C             communication between your program that calls DDEABM and
  C             the  DF subroutine.  They are not used or altered by
  C             DDEABM.  If you do not need RPAR or IPAR, ignore these
  C             parameters by treating them as dummy arguments.  If you do
  C             choose to use them, dimension them in your calling program
  C             and in DF as arrays of appropriate length.
  C
  C **********************************************************************
  C * OUTPUT -- After Any Return From DDEABM *
  C **********************************************************************
  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             You may also be interested in the approximate derivative
  C             of the solution at T.  It is contained in
  C             RWORK(21),...,RWORK(20+NEQ).
  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
  C                         *** Task Interrupted ***
  C                   Reported by negative values of IDID
  C
  C             IDID = -1 -- A large amount of work has been expended.
  C                       (500 steps attempted)
  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 = -4 -- The problem appears to be stiff.
  C
  C             IDID = -5,-6,-7,..,-32  -- Not applicable for this code
  C                       but used by other members of DEPAC or possible
  C                       future extensions.
  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(11)--which contains the step size H to be
  C                        attempted on the next step.
  C
  C             RWORK(12)--if the tolerances have been increased by the
  C                        code (IDID = -2) , they were multiplied by the
  C                        value in RWORK(12).
  C
  C             RWORK(13)--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(20+I)--Which contains the approximate derivative
  C                        of the solution component Y(I).  In DDEABM, it
  C                        is obtained by calling subroutine DF to
  C                        evaluate the differential equation using T and
  C                        Y(*) when IDID=1 or 2, and by interpolation
  C                        when IDID=3.
  C
  C **********************************************************************
  C * INPUT -- What To Do To Continue The Integration *
  C *             (calls after the first)             *
  C **********************************************************************
  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(*), RWORK(*), IWORK(*) or
  C        the differential equation in subroutine DF. Any such alteration
  C        constitutes a new problem and must be treated as such, i.e.
  C        you must start afresh.
  C
  C        You cannot change from vector to scalar error control or vice
  C        versa (INFO(2)) but you can change the size of the entries of
  C        RTOL, ATOL.  Increasing a tolerance makes the equation easier
  C        to integrate.  Decreasing a tolerance will make the equation
  C        harder to integrate and should generally be avoided.
  C
  C        You can switch from the intermediate-output mode to the
  C        interval mode (INFO(3)) or vice versa at any time.
  C
  C        If it has been necessary to prevent the integration from going
  C        past a point TSTOP (INFO(4), RWORK(1)), keep in mind that the
  C        code will not integrate to any TOUT beyond the currently
  C        specified TSTOP.  Once TSTOP has been reached you must change
  C        the value of TSTOP or set INFO(4)=0.  You may change INFO(4)
  C        or TSTOP at any time but you must supply the value of TSTOP in
  C        RWORK(1) whenever you set INFO(4)=1.
  C
  C        The parameter INFO(1) is used by the code to indicate the
  C        beginning of a new problem and to indicate whether integration
  C        is to be continued.  You must input the value  INFO(1) = 0
  C        when starting a new problem.  You must input the value
  C        INFO(1) = 1  if you wish to continue after an interrupted task.
  C        Do not set  INFO(1) = 0  on a continuation call unless you
  C        want the code to restart at the current T.
  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 reset INFO(1) = 1
  C         If
  C             IDID = -1, the code has attempted 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             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, the problem appears to be stiff.  It is very
  C                     inefficient to solve such problems with DDEABM.
  C                     The code DDEBDF in DEPAC handles this task
  C                     efficiently.  If you are absolutely sure you want
  C                     to continue with DDEABM, set INFO(1)=1 and call
  C                     the code again.
  C
  C             IDID = -5,-6,-7,..,-32  --- cannot occur with this code
  C                     but used by other members of DEPAC or possible
  C                     future extensions.
  C
  C                         *** Following A Terminated Task ***
  C         If
  C             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 *Long Description:
  C
  C **********************************************************************
  C *             DEPAC Package Overview           *
  C **********************************************************************
  C
  C ....   You have a choice of three differential equation solvers from
  C ....   DEPAC. The following brief descriptions are meant to aid you in
  C ....   choosing the most appropriate code for your problem.
  C
  C ....   DDERKF is a fifth order Runge-Kutta code. It is the simplest of
  C ....   the three choices, both algorithmically and in the use of the
  C ....   code. DDERKF is primarily designed to solve non-stiff and
  C ....   mildly stiff differential equations when derivative evaluations
  C ....   are not expensive. It should generally not be used to get high
  C ....   accuracy results nor answers at a great many specific points.
  C ....   Because DDERKF has very low overhead costs, it will usually
  C ....   result in the least expensive integration when solving
  C ....   problems requiring a modest amount of accuracy and having
  C ....   equations that are not costly to evaluate. DDERKF attempts to
  C ....   discover when it is not suitable for the task posed.
  C
  C ....   DDEABM is a variable order (one through twelve) Adams code.
  C ....   Its complexity lies somewhere between that of DDERKF and
  C ....   DDEBDF.  DDEABM is primarily designed to solve non-stiff and
  C ....   mildly stiff differential equations when derivative evaluations
  C ....   are expensive, high accuracy results are needed or answers at
  C ....   many specific points are required. DDEABM attempts to discover
  C ....   when it is not suitable for the task posed.
  C
  C ....   DDEBDF is a variable order (one through five) backward
  C ....   differentiation formula code. it is the most complicated of
  C ....   the three choices. DDEBDF is primarily designed to solve stiff
  C ....   differential equations at crude to moderate tolerances.
  C ....   If the problem is very stiff at all, DDERKF and DDEABM will be
  C ....   quite inefficient compared to DDEBDF. However, DDEBDF will be
  C ....   inefficient compared to DDERKF and DDEABM on non-stiff problems
  C ....   because it uses much more storage, has a much larger overhead,
  C ....   and the low order formulas will not give high accuracies
  C ....   efficiently.
  C
  C ....   The concept of stiffness cannot be described in a few words.
  C ....   If you do not know the problem to be stiff, try either DDERKF
  C ....   or DDEABM. Both of these codes will inform you of stiffness
  C ....   when the cost of solving such problems becomes important.
  C
  C *********************************************************************
  C
  C***REFERENCES  L. F. Shampine and H. A. Watts, DEPAC - design of a user
  C                 oriented package of ODE solvers, Report SAND79-2374,
  C                 Sandia Laboratories, 1979.
  C***ROUTINES CALLED  DDES, XERMSG
  C***REVISION HISTORY  (YYMMDD)
  C   820301  DATE WRITTEN
  C   890531  Changed all specific intrinsics to generic.  (WRB)
  C   890831  Modified array declarations.  (WRB)
  C   891006  Cosmetic changes to prologue.  (WRB)
  C   891024  Changed references from DVNORM to DHVNRM.  (WRB)
  C   891024  REVISION DATE from Version 3.2
  C   891214  Prologue converted to Version 4.0 format.  (BAB)
  C   900510  Convert XERRWV calls to XERMSG calls.  (RWC)
  C   920501  Reformatted the REFERENCES section.  (WRB)
  C***END PROLOGUE  DDEABM
  C
        INTEGER IALPHA, IBETA, IDELSN, IDID, IFOURU, IG, IHOLD,
       1      INFO, IP, IPAR, IPHI, IPSI, ISIG, ITOLD, ITSTAR, ITWOU,
       2      IV, IW, IWORK, IWT, IYP, IYPOUT, IYY, LIW, LRW, NEQ
        DOUBLE PRECISION ATOL, RPAR, RTOL, RWORK, T, TOUT, Y
        LOGICAL START,PHASE1,NORND,STIFF,INTOUT
  C
        DIMENSION Y(*),INFO(15),RTOL(*),ATOL(*),RWORK(*),IWORK(*),
       1          RPAR(*),IPAR(*)
  C
        CHARACTER*8 XERN1
        CHARACTER*16 XERN3
  C
        EXTERNAL DF
  C
  C     CHECK FOR AN APPARENT INFINITE LOOP
  C
  C***FIRST EXECUTABLE STATEMENT  DDEABM
        IF ( INFO(1) .EQ. 0 ) IWORK(LIW) = 0
        IF (IWORK(LIW) .GE. 5) THEN
           IF (T .EQ. RWORK(21 + NEQ)) THEN
              WRITE (XERN3, '(1PE15.6)') T
              CALL XERMSG ('SLATEC', 'DDEABM',
       *         'AN APPARENT INFINITE LOOP HAS BEEN DETECTED.$$' //
       *         'YOU HAVE MADE REPEATED CALLS AT T = ' // XERN3 //
       *         ' AND THE INTEGRATION HAS NOT ADVANCED.  CHECK THE ' //
       *         'WAY YOU HAVE SET PARAMETERS FOR THE CALL TO THE ' //
       *         'CODE, PARTICULARLY INFO(1).', 13, 2)
              RETURN
           ENDIF
        ENDIF
  C
  C     CHECK LRW AND LIW FOR SUFFICIENT STORAGE ALLOCATION
  C
        IDID=0
        IF (LRW .LT. 130+21*NEQ) THEN
           WRITE (XERN1, '(I8)') LRW
           CALL XERMSG ('SLATEC', 'DDEABM', 'THE LENGTH OF THE RWORK ' //
       *      'ARRAY MUST BE AT LEAST 130 + 21*NEQ.$$' //
       *      'YOU HAVE CALLED THE CODE WITH LRW = ' // XERN1, 1, 1)
           IDID=-33
        ENDIF
  C
        IF (LIW .LT. 51) THEN
           WRITE (XERN1, '(I8)') LIW
           CALL XERMSG ('SLATEC', 'DDEABM', 'THE LENGTH OF THE IWORK ' //
       *      'ARRAY MUST BE AT LEAST 51.$$YOU HAVE CALLED THE CODE ' //
       *      'WITH LIW = ' // XERN1, 2, 1)
           IDID=-33
        ENDIF
  C
  C     COMPUTE THE INDICES FOR THE ARRAYS TO BE STORED IN THE WORK ARRAY
  C
        IYPOUT = 21
        ITSTAR = NEQ + 21
        IYP = 1 + ITSTAR
        IYY = NEQ + IYP
        IWT = NEQ + IYY
        IP = NEQ + IWT
        IPHI = NEQ + IP
        IALPHA = (NEQ*16) + IPHI
        IBETA = 12 + IALPHA
        IPSI = 12 + IBETA
        IV = 12 + IPSI
        IW = 12 + IV
        ISIG = 12 + IW
        IG = 13 + ISIG
        IGI = 13 + IG
        IXOLD = 11 + IGI
        IHOLD = 1 + IXOLD
        ITOLD = 1 + IHOLD
        IDELSN = 1 + ITOLD
        ITWOU = 1 + IDELSN
        IFOURU = 1 + ITWOU
  C
        RWORK(ITSTAR) = T
        IF (INFO(1) .EQ. 0) GO TO 50
        START = IWORK(21) .NE. (-1)
        PHASE1 = IWORK(22) .NE. (-1)
        NORND = IWORK(23) .NE. (-1)
        STIFF = IWORK(24) .NE. (-1)
        INTOUT = IWORK(25) .NE. (-1)
  C
   50   CALL DDES(DF,NEQ,T,Y,TOUT,INFO,RTOL,ATOL,IDID,RWORK(IYPOUT),
       1         RWORK(IYP),RWORK(IYY),RWORK(IWT),RWORK(IP),RWORK(IPHI),
       2         RWORK(IALPHA),RWORK(IBETA),RWORK(IPSI),RWORK(IV),
       3         RWORK(IW),RWORK(ISIG),RWORK(IG),RWORK(IGI),RWORK(11),
       4         RWORK(12),RWORK(13),RWORK(IXOLD),RWORK(IHOLD),
       5         RWORK(ITOLD),RWORK(IDELSN),RWORK(1),RWORK(ITWOU),
       5         RWORK(IFOURU),START,PHASE1,NORND,STIFF,INTOUT,IWORK(26),
       6         IWORK(27),IWORK(28),IWORK(29),IWORK(30),IWORK(31),
       7         IWORK(32),IWORK(33),IWORK(34),IWORK(35),IWORK(45),
       8         RPAR,IPAR)
  C
        IWORK(21) = -1
        IF (START) IWORK(21) = 1
        IWORK(22) = -1
        IF (PHASE1) IWORK(22) = 1
        IWORK(23) = -1
        IF (NORND) IWORK(23) = 1
        IWORK(24) = -1
        IF (STIFF) IWORK(24) = 1
        IWORK(25) = -1
        IF (INTOUT) IWORK(25) = 1
  C
        IF (IDID .NE. (-2)) IWORK(LIW) = IWORK(LIW) + 1
        IF (T .NE. RWORK(ITSTAR)) IWORK(LIW) = 0
  C
        RETURN
        END
  *DECK DDES
        SUBROUTINE DDES (DF, NEQ, T, Y, TOUT, INFO, RTOL, ATOL, IDID,
       +   YPOUT, YP, YY, WT, P, PHI, ALPHA, BETA, PSI, V, W, SIG, G, GI,
       +   H, EPS, X, XOLD, HOLD, TOLD, DELSGN, TSTOP, TWOU, FOURU, START,
       +   PHASE1, NORND, STIFF, INTOUT, NS, KORD, KOLD, INIT, KSTEPS,
       +   KLE4, IQUIT, KPREV, IVC, IV, KGI, RPAR, IPAR)
  C***BEGIN PROLOGUE  DDES
  C***SUBSIDIARY
  C***PURPOSE  Subsidiary to DDEABM
  C***LIBRARY   SLATEC
  C***TYPE      DOUBLE PRECISION (DES-S, DDES-D)
  C***AUTHOR  Watts, H. A., (SNLA)
  C***DESCRIPTION
  C
  C   DDEABM merely allocates storage for DDES to relieve the user of the
  C   inconvenience of a long call list.  Consequently  DDES  is used as
  C   described in the comments for  DDEABM .
  C
  C***SEE ALSO  DDEABM
  C***ROUTINES CALLED  D1MACH, DINTP, DSTEPS, XERMSG
  C***REVISION HISTORY  (YYMMDD)
  C   820301  DATE WRITTEN
  C   890531  Changed all specific intrinsics to generic.  (WRB)
  C   890831  Modified array declarations.  (WRB)
  C   891214  Prologue converted to Version 4.0 format.  (BAB)
  C   900328  Added TYPE section.  (WRB)
  C   900510  Convert XERRWV calls to XERMSG calls, cvt GOTOs to
  C           IF-THEN-ELSE.  (RWC)
  C   910722  Updated AUTHOR section.  (ALS)
  C***END PROLOGUE  DDES
  C
        INTEGER IDID, INFO, INIT, IPAR, IQUIT, IV, IVC, K, KGI, KLE4,
       1      KOLD, KORD, KPREV, KSTEPS, L, LTOL, MAXNUM, NATOLP, NEQ,
       2      NRTOLP, NS
        DOUBLE PRECISION A, ABSDEL, ALPHA, ATOL, BETA, D1MACH,
       1      DEL, DELSGN, DT, EPS, FOURU, G, GI, H,
       2      HA, HOLD, P, PHI, PSI, RPAR, RTOL, SIG, T, TOLD, TOUT,
       3      TSTOP, TWOU, U, V, W, WT, X, XOLD, Y, YP, YPOUT, YY
        LOGICAL STIFF,CRASH,START,PHASE1,NORND,INTOUT
  C
        DIMENSION Y(*),YY(*),WT(*),PHI(NEQ,16),P(*),YP(*),
       1  YPOUT(*),PSI(12),ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13),
       2  GI(11),IV(10),INFO(15),RTOL(*),ATOL(*),RPAR(*),IPAR(*)
        CHARACTER*8 XERN1
        CHARACTER*16 XERN3, XERN4
  C
        EXTERNAL DF
  C
  C.......................................................................
  C
  C  THE EXPENSE OF SOLVING THE PROBLEM IS MONITORED BY COUNTING THE
  C  NUMBER OF  STEPS ATTEMPTED. WHEN THIS EXCEEDS  MAXNUM, THE COUNTER
  C  IS RESET TO ZERO AND THE USER IS INFORMED ABOUT POSSIBLE EXCESSIVE
  C  WORK.
  C
        SAVE MAXNUM
        DATA MAXNUM/500/
  C
  C.......................................................................
  C
  C***FIRST EXECUTABLE STATEMENT  DDES
        IF (INFO(1) .EQ. 0) THEN
  C
  C ON THE FIRST CALL , PERFORM INITIALIZATION --
  C        DEFINE THE MACHINE UNIT ROUNDOFF QUANTITY  U  BY CALLING THE
  C        FUNCTION ROUTINE  D1MACH. THE USER MUST MAKE SURE THAT THE
  C        VALUES SET IN D1MACH ARE RELEVANT TO THE COMPUTER BEING USED.
  C
           U=D1MACH(4)
  C                       -- SET ASSOCIATED MACHINE DEPENDENT PARAMETERS
           TWOU=2.D0*U
           FOURU=4.D0*U
  C                       -- SET TERMINATION FLAG
           IQUIT=0
  C                       -- SET INITIALIZATION INDICATOR
           INIT=0
  C                       -- SET COUNTER FOR ATTEMPTED STEPS
           KSTEPS=0
  C                       -- SET INDICATOR FOR INTERMEDIATE-OUTPUT
           INTOUT= .FALSE.
  C                       -- SET INDICATOR FOR STIFFNESS DETECTION
           STIFF= .FALSE.
  C                       -- SET STEP COUNTER FOR STIFFNESS DETECTION
           KLE4=0
  C                       -- SET INDICATORS FOR STEPS CODE
           START= .TRUE.
           PHASE1= .TRUE.
           NORND= .TRUE.
  C                       -- RESET INFO(1) FOR SUBSEQUENT CALLS
           INFO(1)=1
        ENDIF
  C
  C.......................................................................
  C
  C      CHECK VALIDITY OF INPUT PARAMETERS ON EACH ENTRY
  C
        IF (INFO(1) .NE. 0  .AND.  INFO(1) .NE. 1) THEN
           WRITE (XERN1, '(I8)') INFO(1)
           CALL XERMSG ('SLATEC', 'DDES', 'IN DDEABM, INFO(1) MUST BE ' //
       *      'SET TO 0 FOR THE START OF A NEW PROBLEM, AND MUST BE ' //
       *      'SET TO 1 FOLLOWING AN INTERRUPTED TASK.  YOU ARE ' //
       *      'ATTEMPTING TO CONTINUE THE INTEGRATION ILLEGALLY BY ' //
       *      'CALLING THE CODE WITH INFO(1) = ' // XERN1, 3, 1)
           IDID=-33
        ENDIF
  C
        IF (INFO(2) .NE. 0  .AND.  INFO(2) .NE. 1) THEN
           WRITE (XERN1, '(I8)') INFO(2)
           CALL XERMSG ('SLATEC', 'DDES', 'IN DDEABM, INFO(2) MUST BE ' //
       *      '0 OR 1 INDICATING SCALAR AND VECTOR ERROR TOLERANCES, ' //
       *      'RESPECTIVELY.  YOU HAVE CALLED THE CODE WITH INFO(2) = ' //
       *      XERN1, 4, 1)
           IDID=-33
        ENDIF
  C
        IF (INFO(3) .NE. 0  .AND.  INFO(3) .NE. 1) THEN
           WRITE (XERN1, '(I8)') INFO(3)
           CALL XERMSG ('SLATEC', 'DDES', 'IN DDEABM, INFO(3) MUST BE ' //
       *      '0 OR 1 INDICATING THE INTERVAL OR INTERMEDIATE-OUTPUT ' //
       *      'MODE OF INTEGRATION, RESPECTIVELY.  YOU HAVE CALLED ' //
       *      'THE CODE WITH  INFO(3) = ' // XERN1, 5, 1)
           IDID=-33
        ENDIF
  C
        IF (INFO(4) .NE. 0  .AND.  INFO(4) .NE. 1) THEN
           WRITE (XERN1, '(I8)') INFO(4)
           CALL XERMSG ('SLATEC', 'DDES', 'IN DDEABM, INFO(4) MUST BE ' //
       *      '0 OR 1 INDICATING WHETHER OR NOT THE INTEGRATION ' //
       *      'INTERVAL IS TO BE RESTRICTED BY A POINT TSTOP.  YOU ' //
       *      'HAVE CALLED THE CODE WITH INFO(4) = ' // XERN1, 14, 1)
           IDID=-33
        ENDIF
  C
        IF (NEQ .LT. 1) THEN
           WRITE (XERN1, '(I8)') NEQ
           CALL XERMSG ('SLATEC', 'DDES', 'IN DDEABM,  THE NUMBER OF ' //
       *      'EQUATIONS NEQ MUST BE A POSITIVE INTEGER.  YOU HAVE ' //
       *      'CALLED THE CODE WITH  NEQ = ' // XERN1, 6, 1)
           IDID=-33
        ENDIF
  C
        NRTOLP = 0
        NATOLP = 0
        DO 90 K=1,NEQ
           IF (NRTOLP .EQ. 0 .AND. RTOL(K) .LT. 0.D0) THEN
              WRITE (XERN1, '(I8)') K
              WRITE (XERN3, '(1PE15.6)') RTOL(K)
              CALL XERMSG ('SLATEC', 'DDES', 'IN DDEABM, THE RELATIVE ' //
       *         'ERROR TOLERANCES RTOL MUST BE NON-NEGATIVE.  YOU ' //
       *         'HAVE CALLED THE CODE WITH  RTOL(' // XERN1 // ') = ' //
       *         XERN3 // '.  IN THE CASE OF VECTOR ERROR TOLERANCES, ' //
       *         'NO FURTHER CHECKING OF RTOL COMPONENTS IS DONE.', 7, 1)
              IDID = -33
              NRTOLP = 1
           ENDIF
  C
           IF (NATOLP .EQ. 0 .AND. ATOL(K) .LT. 0.D0) THEN
              WRITE (XERN1, '(I8)') K
              WRITE (XERN3, '(1PE15.6)') ATOL(K)
              CALL XERMSG ('SLATEC', 'DDES', 'IN DDEABM, THE ABSOLUTE ' //
       *         'ERROR TOLERANCES ATOL MUST BE NON-NEGATIVE.  YOU ' //
       *         'HAVE CALLED THE CODE WITH  ATOL(' // XERN1 // ') = ' //
       *         XERN3 // '.  IN THE CASE OF VECTOR ERROR TOLERANCES, ' //
       *         'NO FURTHER CHECKING OF ATOL COMPONENTS IS DONE.', 8, 1)
              IDID = -33
              NATOLP = 1
           ENDIF
  C
           IF (INFO(2) .EQ. 0) GO TO 100
           IF (NATOLP.GT.0 .AND. NRTOLP.GT.0) GO TO 100
     90 CONTINUE
  C
    100 IF (INFO(4) .EQ. 1) THEN
           IF (SIGN(1.D0,TOUT-T) .NE. SIGN(1.D0,TSTOP-T)
       1      .OR. ABS(TOUT-T) .GT. ABS(TSTOP-T)) THEN
              WRITE (XERN3, '(1PE15.6)') TOUT
              WRITE (XERN4, '(1PE15.6)') TSTOP
              CALL XERMSG ('SLATEC', 'DDES', 'IN DDEABM, YOU HAVE ' //
       *         'CALLED THE CODE WITH  TOUT = ' // XERN3 // ' BUT ' //
       *         'YOU HAVE ALSO TOLD THE CODE (INFO(4) = 1) NOT TO ' //
       *         'INTEGRATE PAST THE POINT TSTOP = ' // XERN4 //
       *         ' THESE INSTRUCTIONS CONFLICT.', 14, 1)
              IDID=-33
           ENDIF
        ENDIF
  C
  C     CHECK SOME CONTINUATION POSSIBILITIES
  C
        IF (INIT .NE. 0) THEN
           IF (T .EQ. TOUT) THEN
              WRITE (XERN3, '(1PE15.6)') T
              CALL XERMSG ('SLATEC', 'DDES', 'IN DDEABM, YOU HAVE ' //
       *         'CALLED THE CODE WITH  T = TOUT = ' // XERN3 //
       *         '$$THIS IS NOT ALLOWED ON CONTINUATION CALLS.', 9, 1)
              IDID=-33
           ENDIF
  C
           IF (T .NE. TOLD) THEN
              WRITE (XERN3, '(1PE15.6)') TOLD
              WRITE (XERN4, '(1PE15.6)') T
              CALL XERMSG ('SLATEC', 'DDES', 'IN DDEABM, YOU HAVE ' //
       *         'CHANGED THE VALUE OF T FROM ' // XERN3 // ' TO ' //
       *         XERN4 //'  THIS IS NOT ALLOWED ON CONTINUATION CALLS.',
       *         10, 1)
              IDID=-33
           ENDIF
  C
           IF (INIT .NE. 1) THEN
              IF (DELSGN*(TOUT-T) .LT. 0.D0) THEN
                 WRITE (XERN3, '(1PE15.6)') TOUT
                 CALL XERMSG ('SLATEC', 'DDES', 'IN DDEABM, BY ' //
       *            'CALLING THE CODE WITH TOUT = ' // XERN3 //
       *            ' YOU ARE ATTEMPTING TO CHANGE THE DIRECTION OF ' //
       *            'INTEGRATION.$$THIS IS NOT ALLOWED WITHOUT ' //
       *            'RESTARTING.', 11, 1)
                 IDID=-33
              ENDIF
           ENDIF
        ENDIF
  C
  C     INVALID INPUT DETECTED
  C
        IF (IDID .EQ. (-33)) THEN
           IF (IQUIT .NE. (-33)) THEN
              IQUIT = -33
              INFO(1) = -1
           ELSE
              CALL XERMSG ('SLATEC', 'DDES', 'IN DDEABM, INVALID ' //
       *         'INPUT WAS DETECTED ON SUCCESSIVE ENTRIES.  IT IS ' //
       *         'IMPOSSIBLE TO PROCEED BECAUSE YOU HAVE NOT ' //
       *         'CORRECTED THE PROBLEM, SO EXECUTION IS BEING ' //
       *         'TERMINATED.', 12, 2)
           ENDIF
           RETURN
        ENDIF
  C
  C.......................................................................
  C
  C     RTOL = ATOL = 0. IS ALLOWED AS VALID INPUT AND INTERPRETED AS
  C     ASKING FOR THE MOST ACCURATE SOLUTION POSSIBLE. IN THIS CASE,
  C     THE RELATIVE ERROR TOLERANCE RTOL IS RESET TO THE SMALLEST VALUE
  C     FOURU WHICH IS LIKELY TO BE REASONABLE FOR THIS METHOD AND MACHINE
  C
        DO 180 K=1,NEQ
          IF (RTOL(K)+ATOL(K) .GT. 0.D0) GO TO 170
          RTOL(K)=FOURU
          IDID=-2
    170   IF (INFO(2) .EQ. 0) GO TO 190
    180   CONTINUE
  C
    190 IF (IDID .NE. (-2)) GO TO 200
  C                       RTOL=ATOL=0 ON INPUT, SO RTOL IS CHANGED TO A
  C                                                SMALL POSITIVE VALUE
        INFO(1)=-1
        RETURN
  C
  C     BRANCH ON STATUS OF INITIALIZATION INDICATOR
  C            INIT=0 MEANS INITIAL DERIVATIVES AND NOMINAL STEP SIZE
  C                   AND DIRECTION NOT YET SET
  C            INIT=1 MEANS NOMINAL STEP SIZE AND DIRECTION NOT YET SET
  C            INIT=2 MEANS NO FURTHER INITIALIZATION REQUIRED
  C
    200 IF (INIT .EQ. 0) GO TO 210
        IF (INIT .EQ. 1) GO TO 220
        GO TO 240
  C
  C.......................................................................
  C
  C     MORE INITIALIZATION --
  C                         -- EVALUATE INITIAL DERIVATIVES
  C
    210 INIT=1
        A=T
        CALL DF(A,Y,YP,RPAR,IPAR)
        IF (T .NE. TOUT) GO TO 220
        IDID=2
        DO 215 L = 1,NEQ
    215    YPOUT(L) = YP(L)
        TOLD=T
        RETURN
  C
  C                         -- SET INDEPENDENT AND DEPENDENT VARIABLES
  C                                              X AND YY(*) FOR STEPS
  C                         -- SET SIGN OF INTEGRATION DIRECTION
  C                         -- INITIALIZE THE STEP SIZE
  C
    220 INIT = 2
        X = T
        DO 230 L = 1,NEQ
    230   YY(L) = Y(L)
        DELSGN = SIGN(1.0D0,TOUT-T)
        H = SIGN(MAX(FOURU*ABS(X),ABS(TOUT-X)),TOUT-X)
  C
  C.......................................................................
  C
  C   ON EACH CALL SET INFORMATION WHICH DETERMINES THE ALLOWED INTERVAL
  C   OF INTEGRATION BEFORE RETURNING WITH AN ANSWER AT TOUT
  C
    240 DEL = TOUT - T
        ABSDEL = ABS(DEL)
  C
  C.......................................................................
  C
  C   IF ALREADY PAST OUTPUT POINT, INTERPOLATE AND RETURN
  C
    250 IF(ABS(X-T) .LT. ABSDEL) GO TO 260
        CALL DINTP(X,YY,TOUT,Y,YPOUT,NEQ,KOLD,PHI,IVC,IV,KGI,GI,
       1                                        ALPHA,G,W,XOLD,P)
        IDID = 3
        IF (X .NE. TOUT) GO TO 255
        IDID = 2
        INTOUT = .FALSE.
    255 T = TOUT
        TOLD = T
        RETURN
  C
  C   IF CANNOT GO PAST TSTOP AND SUFFICIENTLY CLOSE,
  C   EXTRAPOLATE AND RETURN
  C
    260 IF (INFO(4) .NE. 1) GO TO 280
        IF (ABS(TSTOP-X) .GE. FOURU*ABS(X)) GO TO 280
        DT = TOUT - X
        DO 270 L = 1,NEQ
    270   Y(L) = YY(L) + DT*YP(L)
        CALL DF(TOUT,Y,YPOUT,RPAR,IPAR)
        IDID = 3
        T = TOUT
        TOLD = T
        RETURN
  C
    280 IF (INFO(3) .EQ. 0  .OR.  .NOT.INTOUT) GO TO 300
  C
  C   INTERMEDIATE-OUTPUT MODE
  C
        IDID = 1
        DO 290 L = 1,NEQ
          Y(L)=YY(L)
    290   YPOUT(L) = YP(L)
        T = X
        TOLD = T
        INTOUT = .FALSE.
        RETURN
  C
  C.......................................................................
  C
  C     MONITOR NUMBER OF STEPS ATTEMPTED
  C
    300 IF (KSTEPS .LE. MAXNUM) GO TO 330
  C
  C                       A SIGNIFICANT AMOUNT OF WORK HAS BEEN EXPENDED
        IDID=-1
        KSTEPS=0
        IF (.NOT. STIFF) GO TO 310
  C
  C                       PROBLEM APPEARS TO BE STIFF
        IDID=-4
        STIFF= .FALSE.
        KLE4=0
  C
    310 DO 320 L = 1,NEQ
          Y(L) = YY(L)
    320   YPOUT(L) = YP(L)
        T = X
        TOLD = T
        INFO(1) = -1
        INTOUT = .FALSE.
        RETURN
  C
  C.......................................................................
  C
  C   LIMIT STEP SIZE, SET WEIGHT VECTOR AND TAKE A STEP
  C
    330 HA = ABS(H)
        IF (INFO(4) .NE. 1) GO TO 340
        HA = MIN(HA,ABS(TSTOP-X))
    340 H = SIGN(HA,H)
        EPS = 1.0D0
        LTOL = 1
        DO 350 L = 1,NEQ
          IF (INFO(2) .EQ. 1) LTOL = L
          WT(L) = RTOL(LTOL)*ABS(YY(L)) + ATOL(LTOL)
          IF (WT(L) .LE. 0.0D0) GO TO 360
    350   CONTINUE
        GO TO 380
  C
  C                       RELATIVE ERROR CRITERION INAPPROPRIATE
    360 IDID = -3
        DO 370 L = 1,NEQ
          Y(L) = YY(L)
    370   YPOUT(L) = YP(L)
        T = X
        TOLD = T
        INFO(1) = -1
        INTOUT = .FALSE.
        RETURN
  C
    380 CALL DSTEPS(DF,NEQ,YY,X,H,EPS,WT,START,HOLD,KORD,KOLD,CRASH,PHI,P,
       1           YP,PSI,ALPHA,BETA,SIG,V,W,G,PHASE1,NS,NORND,KSTEPS,
       2           TWOU,FOURU,XOLD,KPREV,IVC,IV,KGI,GI,RPAR,IPAR)
  C
  C.......................................................................
  C
        IF(.NOT.CRASH) GO TO 420
  C
  C                       TOLERANCES TOO SMALL
        IDID = -2
        RTOL(1) = EPS*RTOL(1)
        ATOL(1) = EPS*ATOL(1)
        IF (INFO(2) .EQ. 0) GO TO 400
        DO 390 L = 2,NEQ
          RTOL(L) = EPS*RTOL(L)
    390   ATOL(L) = EPS*ATOL(L)
    400 DO 410 L = 1,NEQ
          Y(L) = YY(L)
    410   YPOUT(L) = YP(L)
        T = X
        TOLD = T
        INFO(1) = -1
        INTOUT = .FALSE.
        RETURN
  C
  C   (STIFFNESS TEST) COUNT NUMBER OF CONSECUTIVE STEPS TAKEN WITH THE
  C   ORDER OF THE METHOD BEING LESS OR EQUAL TO FOUR
  C
    420 KLE4 = KLE4 + 1
        IF(KOLD .GT. 4) KLE4 = 0
        IF(KLE4 .GE. 50) STIFF = .TRUE.
        INTOUT = .TRUE.
        GO TO 250
        END
  *DECK DINTP
        SUBROUTINE DINTP (X, Y, XOUT, YOUT, YPOUT, NEQN, KOLD, PHI, IVC,
       +   IV, KGI, GI, ALPHA, OG, OW, OX, OY)
  C***BEGIN PROLOGUE  DINTP
  C***PURPOSE  Approximate the solution at XOUT by evaluating the
  C            polynomial computed in DSTEPS at XOUT.  Must be used in
  C            conjunction with DSTEPS.
  C***LIBRARY   SLATEC (DEPAC)
  C***CATEGORY  I1A1B
  C***TYPE      DOUBLE PRECISION (SINTRP-S, DINTP-D)
  C***KEYWORDS  ADAMS METHOD, DEPAC, INITIAL VALUE PROBLEMS, ODE,
  C             ORDINARY DIFFERENTIAL EQUATIONS, PREDICTOR-CORRECTOR,
  C             SMOOTH INTERPOLANT
  C***AUTHOR  Watts, H. A., (SNLA)
  C***DESCRIPTION
  C
  C   The methods in subroutine  DSTEPS  approximate the solution near  X
  C   by a polynomial.  Subroutine  DINTP  approximates the solution at
  C   XOUT  by evaluating the polynomial there.  Information defining this
  C   polynomial is passed from  DSTEPS  so  DINTP  cannot be used alone.
  C
  C   Subroutine DSTEPS is completely explained and documented in the text
  C   "Computer Solution of Ordinary Differential Equations, the Initial
  C   Value Problem"  by L. F. Shampine and M. K. Gordon.
  C
  C   Input to DINTP --
  C
  C   The user provides storage in the calling program for the arrays in
  C   the call list
  C      DIMENSION Y(NEQN),YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16),OY(NEQN)
  C                AND ALPHA(12),OG(13),OW(12),GI(11),IV(10)
  C   and defines
  C      XOUT -- point at which solution is desired.
  C   The remaining parameters are defined in  DSTEPS  and passed to
  C   DINTP  from that subroutine
  C
  C   Output from  DINTP --
  C
  C      YOUT(*) -- solution at  XOUT
  C      YPOUT(*) -- derivative of solution at  XOUT
  C   The remaining parameters are returned unaltered from their input
  C   values.  Integration with  DSTEPS  may be continued.
  C
  C***REFERENCES  H. A. Watts, A smoother interpolant for DE/STEP, INTRP
  C                 II, Report SAND84-0293, Sandia Laboratories, 1984.
  C***ROUTINES CALLED  (NONE)
  C***REVISION HISTORY  (YYMMDD)
  C   840201  DATE WRITTEN
  C   890831  Modified array declarations.  (WRB)
  C   890831  REVISION DATE from Version 3.2
  C   891214  Prologue converted to Version 4.0 format.  (BAB)
  C   920501  Reformatted the REFERENCES section.  (WRB)
  C***END PROLOGUE  DINTP
  C
        INTEGER I, IQ, IV, IVC, IW, J, JQ, KGI, KOLD, KP1, KP2,
       1        L, M, NEQN
        DOUBLE PRECISION ALP, ALPHA, C, G, GDI, GDIF, GI, GAMMA, H, HI,
       1       HMU, OG, OW, OX, OY, PHI, RMU, SIGMA, TEMP1, TEMP2, TEMP3,
       2       W, X, XI, XIM1, XIQ, XOUT, Y, YOUT, YPOUT
  C
        DIMENSION Y(*),YOUT(*),YPOUT(*),PHI(NEQN,16),OY(*)
        DIMENSION G(13),C(13),W(13),OG(13),OW(12),ALPHA(12),GI(11),IV(10)
  C
  C***FIRST EXECUTABLE STATEMENT  DINTP
        KP1 = KOLD + 1
        KP2 = KOLD + 2
  C
        HI = XOUT - OX
        H = X - OX
        XI = HI/H
        XIM1 = XI - 1.D0
  C
  C   INITIALIZE W(*) FOR COMPUTING G(*)
  C
        XIQ = XI
        DO 10 IQ = 1,KP1
          XIQ = XI*XIQ
          TEMP1 = IQ*(IQ+1)
   10     W(IQ) = XIQ/TEMP1
  C
  C   COMPUTE THE DOUBLE INTEGRAL TERM GDI
  C
        IF (KOLD .LE. KGI) GO TO 50
        IF (IVC .GT. 0) GO TO 20
        GDI = 1.0D0/TEMP1
        M = 2
        GO TO 30
   20   IW = IV(IVC)
        GDI = OW(IW)
        M = KOLD - IW + 3
   30   IF (M .GT. KOLD) GO TO 60
        DO 40 I = M,KOLD
   40     GDI = OW(KP2-I) - ALPHA(I)*GDI
        GO TO 60
   50   GDI = GI(KOLD)
  C
  C   COMPUTE G(*) AND C(*)
  C
   60   G(1) = XI
        G(2) = 0.5D0*XI*XI
        C(1) = 1.0D0
        C(2) = XI
        IF (KOLD .LT. 2) GO TO 90
        DO 80 I = 2,KOLD
          ALP = ALPHA(I)
          GAMMA = 1.0D0 + XIM1*ALP
          L = KP2 - I
          DO 70 JQ = 1,L
   70       W(JQ) = GAMMA*W(JQ) - ALP*W(JQ+1)
          G(I+1) = W(1)
   80     C(I+1) = GAMMA*C(I)
  C
  C   DEFINE INTERPOLATION PARAMETERS
  C
   90   SIGMA = (W(2) - XIM1*W(1))/GDI
        RMU = XIM1*C(KP1)/GDI
        HMU = RMU/H
  C
  C   INTERPOLATE FOR THE SOLUTION -- YOUT
  C   AND FOR THE DERIVATIVE OF THE SOLUTION -- YPOUT
  C
        DO 100 L = 1,NEQN
          YOUT(L) = 0.0D0
   100    YPOUT(L) = 0.0D0
        DO 120 J = 1,KOLD
          I = KP2 - J
          GDIF = OG(I) - OG(I-1)
          TEMP2 = (G(I) - G(I-1)) - SIGMA*GDIF
          TEMP3 = (C(I) - C(I-1)) + RMU*GDIF
          DO 110 L = 1,NEQN
            YOUT(L) = YOUT(L) + TEMP2*PHI(L,I)
   110      YPOUT(L) = YPOUT(L) + TEMP3*PHI(L,I)
   120    CONTINUE
        DO 130 L = 1,NEQN
          YOUT(L) = ((1.0D0 - SIGMA)*OY(L) + SIGMA*Y(L)) +
       1             H*(YOUT(L) + (G(1) - SIGMA*OG(1))*PHI(L,1))
   130    YPOUT(L) = HMU*(OY(L) - Y(L)) +
       1                (YPOUT(L) + (C(1) + RMU*OG(1))*PHI(L,1))
  C
        RETURN
        END
  *DECK XERMSG
        SUBROUTINE XERMSG (LIBRAR, SUBROU, MESSG, NERR, LEVEL)
  C***BEGIN PROLOGUE  XERMSG
  C***PURPOSE  Process error messages for SLATEC and other libraries.
  C***LIBRARY   SLATEC (XERROR)
  C***CATEGORY  R3C
  C***TYPE      ALL (XERMSG-A)
  C***KEYWORDS  ERROR MESSAGE, XERROR
  C***AUTHOR  Fong, Kirby, (NMFECC at LLNL)
  C***DESCRIPTION
  C
  C   XERMSG processes a diagnostic message in a manner determined by the
  C   value of LEVEL and the current value of the library error control
  C   flag, KONTRL.  See subroutine XSETF for details.
  C
  C    LIBRAR   A character constant (or character variable) with the name
  C             of the library.  This will be 'SLATEC' for the SLATEC
  C             Common Math Library.  The error handling package is
  C             general enough to be used by many libraries
  C             simultaneously, so it is desirable for the routine that
  C             detects and reports an error to identify the library name
  C             as well as the routine name.
  C
  C    SUBROU   A character constant (or character variable) with the name
  C             of the routine that detected the error.  Usually it is the
  C             name of the routine that is calling XERMSG.  There are
  C             some instances where a user callable library routine calls
  C             lower level subsidiary routines where the error is
  C             detected.  In such cases it may be more informative to
  C             supply the name of the routine the user called rather than
  C             the name of the subsidiary routine that detected the
  C             error.
  C
  C    MESSG    A character constant (or character variable) with the text
  C             of the error or warning message.  In the example below,
  C             the message is a character constant that contains a
  C             generic message.
  C
  C                   CALL XERMSG ('SLATEC', 'MMPY',
  C                  *'THE ORDER OF THE MATRIX EXCEEDS THE ROW DIMENSION',
  C                  *3, 1)
  C
  C             It is possible (and is sometimes desirable) to generate a
  C             specific message--e.g., one that contains actual numeric
  C             values.  Specific numeric values can be converted into
  C             character strings using formatted WRITE statements into
  C             character variables.  This is called standard Fortran
  C             internal file I/O and is exemplified in the first three
  C             lines of the following example.  You can also catenate
  C             substrings of characters to construct the error message.
  C             Here is an example showing the use of both writing to
  C             an internal file and catenating character strings.
  C
  C                   CHARACTER*5 CHARN, CHARL
  C                   WRITE (CHARN,10) N
  C                   WRITE (CHARL,10) LDA
  C                10 FORMAT(I5)
  C                   CALL XERMSG ('SLATEC', 'MMPY', 'THE ORDER'//CHARN//
  C                  *   ' OF THE MATRIX EXCEEDS ITS ROW DIMENSION OF'//
  C                  *   CHARL, 3, 1)
  C
  C             There are two subtleties worth mentioning.  One is that
  C             the // for character catenation is used to construct the
  C             error message so that no single character constant is
  C             continued to the next line.  This avoids confusion as to
  C             whether there are trailing blanks at the end of the line.
  C             The second is that by catenating the parts of the message
  C             as an actual argument rather than encoding the entire
  C             message into one large character variable, we avoid
  C             having to know how long the message will be in order to
  C             declare an adequate length for that large character
  C             variable.  XERMSG calls XERPRN to print the message using
  C             multiple lines if necessary.  If the message is very long,
  C             XERPRN will break it into pieces of 72 characters (as
  C             requested by XERMSG) for printing on multiple lines.
  C             Also, XERMSG asks XERPRN to prefix each line with ' *  '
  C             so that the total line length could be 76 characters.
  C             Note also that XERPRN scans the error message backwards
  C             to ignore trailing blanks.  Another feature is that
  C             the substring '$$' is treated as a new line sentinel
  C             by XERPRN.  If you want to construct a multiline
  C             message without having to count out multiples of 72
  C             characters, just use '$$' as a separator.  '$$'
  C             obviously must occur within 72 characters of the
  C             start of each line to have its intended effect since
  C             XERPRN is asked to wrap around at 72 characters in
  C             addition to looking for '$$'.
  C
  C    NERR     An integer value that is chosen by the library routine's
  C             author.  It must be in the range -99 to 999 (three
  C             printable digits).  Each distinct error should have its
  C             own error number.  These error numbers should be described
  C             in the machine readable documentation for the routine.
  C             The error numbers need be unique only within each routine,
  C             so it is reasonable for each routine to start enumerating
  C             errors from 1 and proceeding to the next integer.
  C
  C    LEVEL    An integer value in the range 0 to 2 that indicates the
  C             level (severity) of the error.  Their meanings are
  C
  C            -1  A warning message.  This is used if it is not clear
  C                that there really is an error, but the user's attention
  C                may be needed.  An attempt is made to only print this
  C                message once.
  C
  C             0  A warning message.  This is used if it is not clear
  C                that there really is an error, but the user's attention
  C                may be needed.
  C
  C             1  A recoverable error.  This is used even if the error is
  C                so serious that the routine cannot return any useful
  C                answer.  If the user has told the error package to
  C                return after recoverable errors, then XERMSG will
  C                return to the Library routine which can then return to
  C                the user's routine.  The user may also permit the error
  C                package to terminate the program upon encountering a
  C                recoverable error.
  C
  C             2  A fatal error.  XERMSG will not return to its caller
  C                after it receives a fatal error.  This level should
  C                hardly ever be used; it is much better to allow the
  C                user a chance to recover.  An example of one of the few
  C                cases in which it is permissible to declare a level 2
  C                error is a reverse communication Library routine that
  C                is likely to be called repeatedly until it integrates
  C                across some interval.  If there is a serious error in
  C                the input such that another step cannot be taken and
  C                the Library routine is called again without the input
  C                error having been corrected by the caller, the Library
  C                routine will probably be called forever with improper
  C                input.  In this case, it is reasonable to declare the
  C                error to be fatal.
  C
  C    Each of the arguments to XERMSG is input; none will be modified by
  C    XERMSG.  A routine may make multiple calls to XERMSG with warning
  C    level messages; however, after a call to XERMSG with a recoverable
  C    error, the routine should return to the user.  Do not try to call
  C    XERMSG with a second recoverable error after the first recoverable
  C    error because the error package saves the error number.  The user
  C    can retrieve this error number by calling another entry point in
  C    the error handling package and then clear the error number when
  C    recovering from the error.  Calling XERMSG in succession causes the
  C    old error number to be overwritten by the latest error number.
  C    This is considered harmless for error numbers associated with
  C    warning messages but must not be done for error numbers of serious
  C    errors.  After a call to XERMSG with a recoverable error, the user
  C    must be given a chance to call NUMXER or XERCLR to retrieve or
  C    clear the error number.
  C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
  C                 Error-handling Package, SAND82-0800, Sandia
  C                 Laboratories, 1982.
  C***ROUTINES CALLED  FDUMP, J4SAVE, XERCNT, XERHLT, XERPRN, XERSVE
  C***REVISION HISTORY  (YYMMDD)
  C   880101  DATE WRITTEN
  C   880621  REVISED AS DIRECTED AT SLATEC CML MEETING OF FEBRUARY 1988.
  C           THERE ARE TWO BASIC CHANGES.
  C           1.  A NEW ROUTINE, XERPRN, IS USED INSTEAD OF XERPRT TO
  C               PRINT MESSAGES.  THIS ROUTINE WILL BREAK LONG MESSAGES
  C               INTO PIECES FOR PRINTING ON MULTIPLE LINES.  '$$' IS
  C               ACCEPTED AS A NEW LINE SENTINEL.  A PREFIX CAN BE
  C               ADDED TO EACH LINE TO BE PRINTED.  XERMSG USES EITHER
  C               ' ***' OR ' *  ' AND LONG MESSAGES ARE BROKEN EVERY
  C               72 CHARACTERS (AT MOST) SO THAT THE MAXIMUM LINE
  C               LENGTH OUTPUT CAN NOW BE AS GREAT AS 76.
  C           2.  THE TEXT OF ALL MESSAGES IS NOW IN UPPER CASE SINCE THE
  C               FORTRAN STANDARD DOCUMENT DOES NOT ADMIT THE EXISTENCE
  C               OF LOWER CASE.
  C   880708  REVISED AFTER THE SLATEC CML MEETING OF JUNE 29 AND 30.
  C           THE PRINCIPAL CHANGES ARE
  C           1.  CLARIFY COMMENTS IN THE PROLOGUES
  C           2.  RENAME XRPRNT TO XERPRN
  C           3.  REWORK HANDLING OF '$$' IN XERPRN TO HANDLE BLANK LINES
  C               SIMILAR TO THE WAY FORMAT STATEMENTS HANDLE THE /
  C               CHARACTER FOR NEW RECORDS.
  C   890706  REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO
  C           CLEAN UP THE CODING.
  C   890721  REVISED TO USE NEW FEATURE IN XERPRN TO COUNT CHARACTERS IN
  C           PREFIX.
  C   891013  REVISED TO CORRECT COMMENTS.
  C   891214  Prologue converted to Version 4.0 format.  (WRB)
  C   900510  Changed test on NERR to be -9999999 < NERR < 99999999, but
  C           NERR .ne. 0, and on LEVEL to be -2 < LEVEL < 3.  Added
  C           LEVEL=-1 logic, changed calls to XERSAV to XERSVE, and
  C           XERCTL to XERCNT.  (RWC)
  C   920501  Reformatted the REFERENCES section.  (WRB)
  C***END PROLOGUE  XERMSG
        CHARACTER*(*) LIBRAR, SUBROU, MESSG
        CHARACTER*8 XLIBR, XSUBR
        CHARACTER*72  TEMP
        CHARACTER*20  LFIRST
  C***FIRST EXECUTABLE STATEMENT  XERMSG
        LKNTRL = J4SAVE (2, 0, .FALSE.)
        MAXMES = J4SAVE (4, 0, .FALSE.)
  C
  C       LKNTRL IS A LOCAL COPY OF THE CONTROL FLAG KONTRL.
  C       MAXMES IS THE MAXIMUM NUMBER OF TIMES ANY PARTICULAR MESSAGE
  C          SHOULD BE PRINTED.
  C
  C       WE PRINT A FATAL ERROR MESSAGE AND TERMINATE FOR AN ERROR IN
  C          CALLING XERMSG.  THE ERROR NUMBER SHOULD BE POSITIVE,
  C          AND THE LEVEL SHOULD BE BETWEEN 0 AND 2.
  C
        IF (NERR.LT.-9999999 .OR. NERR.GT.99999999 .OR. NERR.EQ.0 .OR.
       *   LEVEL.LT.-1 .OR. LEVEL.GT.2) THEN
           CALL XERPRN (' ***', -1, 'FATAL ERROR IN...$$ ' //
       *      'XERMSG -- INVALID ERROR NUMBER OR LEVEL$$ '//
       *      'JOB ABORT DUE TO FATAL ERROR.', 72)
           CALL XERSVE (' ', ' ', ' ', 0, 0, 0, KDUMMY)
           CALL XERHLT (' ***XERMSG -- INVALID INPUT')
           RETURN
        ENDIF
  C
  C       RECORD THE MESSAGE.
  C
        I = J4SAVE (1, NERR, .TRUE.)
        CALL XERSVE (LIBRAR, SUBROU, MESSG, 1, NERR, LEVEL, KOUNT)
  C
  C       HANDLE PRINT-ONCE WARNING MESSAGES.
  C
        IF (LEVEL.EQ.-1 .AND. KOUNT.GT.1) RETURN
  C
  C       ALLOW TEMPORARY USER OVERRIDE OF THE CONTROL FLAG.
  C
        XLIBR  = LIBRAR
        XSUBR  = SUBROU
        LFIRST = MESSG
        LERR   = NERR
        LLEVEL = LEVEL
        CALL XERCNT (XLIBR, XSUBR, LFIRST, LERR, LLEVEL, LKNTRL)
  C
        LKNTRL = MAX(-2, MIN(2,LKNTRL))
        MKNTRL = ABS(LKNTRL)
  C
  C       SKIP PRINTING IF THE CONTROL FLAG VALUE AS RESET IN XERCNT IS
  C       ZERO AND THE ERROR IS NOT FATAL.
  C
        IF (LEVEL.LT.2 .AND. LKNTRL.EQ.0) GO TO 30
        IF (LEVEL.EQ.0 .AND. KOUNT.GT.MAXMES) GO TO 30
        IF (LEVEL.EQ.1 .AND. KOUNT.GT.MAXMES .AND. MKNTRL.EQ.1) GO TO 30
        IF (LEVEL.EQ.2 .AND. KOUNT.GT.MAX(1,MAXMES)) GO TO 30
  C
  C       ANNOUNCE THE NAMES OF THE LIBRARY AND SUBROUTINE BY BUILDING A
  C       MESSAGE IN CHARACTER VARIABLE TEMP (NOT EXCEEDING 66 CHARACTERS)
  C       AND SENDING IT OUT VIA XERPRN.  PRINT ONLY IF CONTROL FLAG
  C       IS NOT ZERO.
  C
        IF (LKNTRL .NE. 0) THEN
           TEMP(1:21) = 'MESSAGE FROM ROUTINE '
           I = MIN(LEN(SUBROU), 16)
           TEMP(22:21+I) = SUBROU(1:I)
           TEMP(22+I:33+I) = ' IN LIBRARY '
           LTEMP = 33 + I
           I = MIN(LEN(LIBRAR), 16)
           TEMP(LTEMP+1:LTEMP+I) = LIBRAR (1:I)
           TEMP(LTEMP+I+1:LTEMP+I+1) = '.'
           LTEMP = LTEMP + I + 1
           CALL XERPRN (' ***', -1, TEMP(1:LTEMP), 72)
        ENDIF
  C
  C       IF LKNTRL IS POSITIVE, PRINT AN INTRODUCTORY LINE BEFORE
  C       PRINTING THE MESSAGE.  THE INTRODUCTORY LINE TELLS THE CHOICE
  C       FROM EACH OF THE FOLLOWING THREE OPTIONS.
  C       1.  LEVEL OF THE MESSAGE
  C              'INFORMATIVE MESSAGE'
  C              'POTENTIALLY RECOVERABLE ERROR'
  C              'FATAL ERROR'
  C       2.  WHETHER CONTROL FLAG WILL ALLOW PROGRAM TO CONTINUE
  C              'PROG CONTINUES'
  C              'PROG ABORTED'
  C       3.  WHETHER OR NOT A TRACEBACK WAS REQUESTED.  (THE TRACEBACK
  C           MAY NOT BE IMPLEMENTED AT SOME SITES, SO THIS ONLY TELLS
  C           WHAT WAS REQUESTED, NOT WHAT WAS DELIVERED.)
  C              'TRACEBACK REQUESTED'
  C              'TRACEBACK NOT REQUESTED'
  C       NOTICE THAT THE LINE INCLUDING FOUR PREFIX CHARACTERS WILL NOT
  C       EXCEED 74 CHARACTERS.
  C       WE SKIP THE NEXT BLOCK IF THE INTRODUCTORY LINE IS NOT NEEDED.
  C
        IF (LKNTRL .GT. 0) THEN
  C
  C       THE FIRST PART OF THE MESSAGE TELLS ABOUT THE LEVEL.
  C
           IF (LEVEL .LE. 0) THEN
              TEMP(1:20) = 'INFORMATIVE MESSAGE,'
              LTEMP = 20
           ELSEIF (LEVEL .EQ. 1) THEN
              TEMP(1:30) = 'POTENTIALLY RECOVERABLE ERROR,'
              LTEMP = 30
           ELSE
              TEMP(1:12) = 'FATAL ERROR,'
              LTEMP = 12
           ENDIF
  C
  C       THEN WHETHER THE PROGRAM WILL CONTINUE.
  C
           IF ((MKNTRL.EQ.2 .AND. LEVEL.GE.1) .OR.
       *       (MKNTRL.EQ.1 .AND. LEVEL.EQ.2)) THEN
              TEMP(LTEMP+1:LTEMP+14) = ' PROG ABORTED,'
              LTEMP = LTEMP + 14
           ELSE
              TEMP(LTEMP+1:LTEMP+16) = ' PROG CONTINUES,'
              LTEMP = LTEMP + 16
           ENDIF
  C
  C       FINALLY TELL WHETHER THERE SHOULD BE A TRACEBACK.
  C
           IF (LKNTRL .GT. 0) THEN
              TEMP(LTEMP+1:LTEMP+20) = ' TRACEBACK REQUESTED'
              LTEMP = LTEMP + 20
           ELSE
              TEMP(LTEMP+1:LTEMP+24) = ' TRACEBACK NOT REQUESTED'
              LTEMP = LTEMP + 24
           ENDIF
           CALL XERPRN (' ***', -1, TEMP(1:LTEMP), 72)
        ENDIF
  C
  C       NOW SEND OUT THE MESSAGE.
  C
        CALL XERPRN (' *  ', -1, MESSG, 72)
  C
  C       IF LKNTRL IS POSITIVE, WRITE THE ERROR NUMBER AND REQUEST A
  C          TRACEBACK.
  C
        IF (LKNTRL .GT. 0) THEN
           WRITE (TEMP, '(''ERROR NUMBER = '', I8)') NERR
           DO 10 I=16,22
              IF (TEMP(I:I) .NE. ' ') GO TO 20
     10    CONTINUE
  C
     20    CALL XERPRN (' *  ', -1, TEMP(1:15) // TEMP(I:23), 72)
           CALL FDUMP
        ENDIF
  C
  C       IF LKNTRL IS NOT ZERO, PRINT A BLANK LINE AND AN END OF MESSAGE.
  C
        IF (LKNTRL .NE. 0) THEN
           CALL XERPRN (' *  ', -1, ' ', 72)
           CALL XERPRN (' ***', -1, 'END OF MESSAGE', 72)
           CALL XERPRN ('    ',  0, ' ', 72)
        ENDIF
  C
  C       IF THE ERROR IS NOT FATAL OR THE ERROR IS RECOVERABLE AND THE
  C       CONTROL FLAG IS SET FOR RECOVERY, THEN RETURN.
  C
     30 IF (LEVEL.LE.0 .OR. (LEVEL.EQ.1 .AND. MKNTRL.LE.1)) RETURN
  C
  C       THE PROGRAM WILL BE STOPPED DUE TO AN UNRECOVERED ERROR OR A
  C       FATAL ERROR.  PRINT THE REASON FOR THE ABORT AND THE ERROR
  C       SUMMARY IF THE CONTROL FLAG AND THE MAXIMUM ERROR COUNT PERMIT.
  C
        IF (LKNTRL.GT.0 .AND. KOUNT.LT.MAX(1,MAXMES)) THEN
           IF (LEVEL .EQ. 1) THEN
              CALL XERPRN
       *         (' ***', -1, 'JOB ABORT DUE TO UNRECOVERED ERROR.', 72)
           ELSE
              CALL XERPRN(' ***', -1, 'JOB ABORT DUE TO FATAL ERROR.', 72)
           ENDIF
           CALL XERSVE (' ', ' ', ' ', -1, 0, 0, KDUMMY)
           CALL XERHLT (' ')
        ELSE
           CALL XERHLT (MESSG)
        ENDIF
        RETURN
        END
  *DECK XERPRN
        SUBROUTINE XERPRN (PREFIX, NPREF, MESSG, NWRAP)
  C***BEGIN PROLOGUE  XERPRN
  C***SUBSIDIARY
  C***PURPOSE  Print error messages processed by XERMSG.
  C***LIBRARY   SLATEC (XERROR)
  C***CATEGORY  R3C
  C***TYPE      ALL (XERPRN-A)
  C***KEYWORDS  ERROR MESSAGES, PRINTING, XERROR
  C***AUTHOR  Fong, Kirby, (NMFECC at LLNL)
  C***DESCRIPTION
  C
  C This routine sends one or more lines to each of the (up to five)
  C logical units to which error messages are to be sent.  This routine
  C is called several times by XERMSG, sometimes with a single line to
  C print and sometimes with a (potentially very long) message that may
  C wrap around into multiple lines.
  C
  C PREFIX  Input argument of type CHARACTER.  This argument contains
  C         characters to be put at the beginning of each line before
  C         the body of the message.  No more than 16 characters of
  C         PREFIX will be used.
  C
  C NPREF   Input argument of type INTEGER.  This argument is the number
  C         of characters to use from PREFIX.  If it is negative, the
  C         intrinsic function LEN is used to determine its length.  If
  C         it is zero, PREFIX is not used.  If it exceeds 16 or if
  C         LEN(PREFIX) exceeds 16, only the first 16 characters will be
  C         used.  If NPREF is positive and the length of PREFIX is less
  C         than NPREF, a copy of PREFIX extended with blanks to length
  C         NPREF will be used.
  C
  C MESSG   Input argument of type CHARACTER.  This is the text of a
  C         message to be printed.  If it is a long message, it will be
  C         broken into pieces for printing on multiple lines.  Each line
  C         will start with the appropriate prefix and be followed by a
  C         piece of the message.  NWRAP is the number of characters per
  C         piece; that is, after each NWRAP characters, we break and
  C         start a new line.  In addition the characters '$$' embedded
  C         in MESSG are a sentinel for a new line.  The counting of
  C         characters up to NWRAP starts over for each new line.  The
  C         value of NWRAP typically used by XERMSG is 72 since many
  C         older error messages in the SLATEC Library are laid out to
  C         rely on wrap-around every 72 characters.
  C
  C NWRAP   Input argument of type INTEGER.  This gives the maximum size
  C         piece into which to break MESSG for printing on multiple
  C         lines.  An embedded '$$' ends a line, and the count restarts
  C         at the following character.  If a line break does not occur
  C         on a blank (it would split a word) that word is moved to the
  C         next line.  Values of NWRAP less than 16 will be treated as
  C         16.  Values of NWRAP greater than 132 will be treated as 132.
  C         The actual line length will be NPREF + NWRAP after NPREF has
  C         been adjusted to fall between 0 and 16 and NWRAP has been
  C         adjusted to fall between 16 and 132.
  C
  C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
  C                 Error-handling Package, SAND82-0800, Sandia
  C                 Laboratories, 1982.
  C***ROUTINES CALLED  I1MACH, XGETUA
  C***REVISION HISTORY  (YYMMDD)
  C   880621  DATE WRITTEN
  C   880708  REVISED AFTER THE SLATEC CML SUBCOMMITTEE MEETING OF
  C           JUNE 29 AND 30 TO CHANGE THE NAME TO XERPRN AND TO REWORK
  C           THE HANDLING OF THE NEW LINE SENTINEL TO BEHAVE LIKE THE
  C           SLASH CHARACTER IN FORMAT STATEMENTS.
  C   890706  REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO
  C           STREAMLINE THE CODING AND FIX A BUG THAT CAUSED EXTRA BLANK
  C           LINES TO BE PRINTED.
  C   890721  REVISED TO ADD A NEW FEATURE.  A NEGATIVE VALUE OF NPREF
  C           CAUSES LEN(PREFIX) TO BE USED AS THE LENGTH.
  C   891013  REVISED TO CORRECT ERROR IN CALCULATING PREFIX LENGTH.
  C   891214  Prologue converted to Version 4.0 format.  (WRB)
  C   900510  Added code to break messages between words.  (RWC)
  C   920501  Reformatted the REFERENCES section.  (WRB)
  C***END PROLOGUE  XERPRN
        CHARACTER*(*) PREFIX, MESSG
        INTEGER NPREF, NWRAP
        CHARACTER*148 CBUFF
        INTEGER IU(5), NUNIT
        CHARACTER*2 NEWLIN
        PARAMETER (NEWLIN = '$$')
  C***FIRST EXECUTABLE STATEMENT  XERPRN
        CALL XGETUA(IU,NUNIT)
  C
  C       A ZERO VALUE FOR A LOGICAL UNIT NUMBER MEANS TO USE THE STANDARD
  C       ERROR MESSAGE UNIT INSTEAD.  I1MACH(4) RETRIEVES THE STANDARD
  C       ERROR MESSAGE UNIT.
  C
        N = I1MACH(4)
        DO 10 I=1,NUNIT
           IF (IU(I) .EQ. 0) IU(I) = N
     10 CONTINUE
  C
  C       LPREF IS THE LENGTH OF THE PREFIX.  THE PREFIX IS PLACED AT THE
  C       BEGINNING OF CBUFF, THE CHARACTER BUFFER, AND KEPT THERE DURING
  C       THE REST OF THIS ROUTINE.
  C
        IF ( NPREF .LT. 0 ) THEN
           LPREF = LEN(PREFIX)
        ELSE
           LPREF = NPREF
        ENDIF
        LPREF = MIN(16, LPREF)
        IF (LPREF .NE. 0) CBUFF(1:LPREF) = PREFIX
  C
  C       LWRAP IS THE MAXIMUM NUMBER OF CHARACTERS WE WANT TO TAKE AT ONE
  C       TIME FROM MESSG TO PRINT ON ONE LINE.
  C
        LWRAP = MAX(16, MIN(132, NWRAP))
  C
  C       SET LENMSG TO THE LENGTH OF MESSG, IGNORE ANY TRAILING BLANKS.
  C
        LENMSG = LEN(MESSG)
        N = LENMSG
        DO 20 I=1,N
           IF (MESSG(LENMSG:LENMSG) .NE. ' ') GO TO 30
           LENMSG = LENMSG - 1
     20 CONTINUE
     30 CONTINUE
  C
  C       IF THE MESSAGE IS ALL BLANKS, THEN PRINT ONE BLANK LINE.
  C
        IF (LENMSG .EQ. 0) THEN
           CBUFF(LPREF+1:LPREF+1) = ' '
           DO 40 I=1,NUNIT
              WRITE(IU(I), '(A)') CBUFF(1:LPREF+1)
     40    CONTINUE
           RETURN
        ENDIF
  C
  C       SET NEXTC TO THE POSITION IN MESSG WHERE THE NEXT SUBSTRING
  C       STARTS.  FROM THIS POSITION WE SCAN FOR THE NEW LINE SENTINEL.
  C       WHEN NEXTC EXCEEDS LENMSG, THERE IS NO MORE TO PRINT.
  C       WE LOOP BACK TO LABEL 50 UNTIL ALL PIECES HAVE BEEN PRINTED.
  C
  C       WE LOOK FOR THE NEXT OCCURRENCE OF THE NEW LINE SENTINEL.  THE
  C       INDEX INTRINSIC FUNCTION RETURNS ZERO IF THERE IS NO OCCURRENCE
  C       OR IF THE LENGTH OF THE FIRST ARGUMENT IS LESS THAN THE LENGTH
  C       OF THE SECOND ARGUMENT.
  C
  C       THERE ARE SEVERAL CASES WHICH SHOULD BE CHECKED FOR IN THE
  C       FOLLOWING ORDER.  WE ARE ATTEMPTING TO SET LPIECE TO THE NUMBER
  C       OF CHARACTERS THAT SHOULD BE TAKEN FROM MESSG STARTING AT
  C       POSITION NEXTC.
  C
  C       LPIECE .EQ. 0   THE NEW LINE SENTINEL DOES NOT OCCUR IN THE
  C                       REMAINDER OF THE CHARACTER STRING.  LPIECE
  C                       SHOULD BE SET TO LWRAP OR LENMSG+1-NEXTC,
  C                       WHICHEVER IS LESS.
  C
  C       LPIECE .EQ. 1   THE NEW LINE SENTINEL STARTS AT MESSG(NEXTC:
  C                       NEXTC).  LPIECE IS EFFECTIVELY ZERO, AND WE
  C                       PRINT NOTHING TO AVOID PRODUCING UNNECESSARY
  C                       BLANK LINES.  THIS TAKES CARE OF THE SITUATION
  C                       WHERE THE LIBRARY ROUTINE HAS A MESSAGE OF
  C                       EXACTLY 72 CHARACTERS FOLLOWED BY A NEW LINE
  C                       SENTINEL FOLLOWED BY MORE CHARACTERS.  NEXTC
  C                       SHOULD BE INCREMENTED BY 2.
  C
  C       LPIECE .GT. LWRAP+1  REDUCE LPIECE TO LWRAP.
  C
  C       ELSE            THIS LAST CASE MEANS 2 .LE. LPIECE .LE. LWRAP+1
  C                       RESET LPIECE = LPIECE-1.  NOTE THAT THIS
  C                       PROPERLY HANDLES THE END CASE WHERE LPIECE .EQ.
  C                       LWRAP+1.  THAT IS, THE SENTINEL FALLS EXACTLY
  C                       AT THE END OF A LINE.
  C
        NEXTC = 1
     50 LPIECE = INDEX(MESSG(NEXTC:LENMSG), NEWLIN)
        IF (LPIECE .EQ. 0) THEN
  C
  C       THERE WAS NO NEW LINE SENTINEL FOUND.
  C
           IDELTA = 0
           LPIECE = MIN(LWRAP, LENMSG+1-NEXTC)
           IF (LPIECE .LT. LENMSG+1-NEXTC) THEN
              DO 52 I=LPIECE+1,2,-1
                 IF (MESSG(NEXTC+I-1:NEXTC+I-1) .EQ. ' ') THEN
                    LPIECE = I-1
                    IDELTA = 1
                    GOTO 54
                 ENDIF
     52       CONTINUE
           ENDIF
     54    CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1)
           NEXTC = NEXTC + LPIECE + IDELTA
        ELSEIF (LPIECE .EQ. 1) THEN
  C
  C       WE HAVE A NEW LINE SENTINEL AT MESSG(NEXTC:NEXTC+1).
  C       DON'T PRINT A BLANK LINE.
  C
           NEXTC = NEXTC + 2
           GO TO 50
        ELSEIF (LPIECE .GT. LWRAP+1) THEN
  C
  C       LPIECE SHOULD BE SET DOWN TO LWRAP.
  C
           IDELTA = 0
           LPIECE = LWRAP
           DO 56 I=LPIECE+1,2,-1
              IF (MESSG(NEXTC+I-1:NEXTC+I-1) .EQ. ' ') THEN
                 LPIECE = I-1
                 IDELTA = 1
                 GOTO 58
              ENDIF
     56    CONTINUE
     58    CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1)
           NEXTC = NEXTC + LPIECE + IDELTA
        ELSE
  C
  C       IF WE ARRIVE HERE, IT MEANS 2 .LE. LPIECE .LE. LWRAP+1.
  C       WE SHOULD DECREMENT LPIECE BY ONE.
  C
           LPIECE = LPIECE - 1
           CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1)
           NEXTC  = NEXTC + LPIECE + 2
        ENDIF
  C
  C       PRINT
  C
        DO 60 I=1,NUNIT
           WRITE(IU(I), '(A)') CBUFF(1:LPREF+LPIECE)
     60 CONTINUE
  C
        IF (NEXTC .LE. LENMSG) GO TO 50
        RETURN
        END
  *DECK XERSVE
        SUBROUTINE XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL,
       +   ICOUNT)
  C***BEGIN PROLOGUE  XERSVE
  C***SUBSIDIARY
  C***PURPOSE  Record that an error has occurred.
  C***LIBRARY   SLATEC (XERROR)
  C***CATEGORY  R3
  C***TYPE      ALL (XERSVE-A)
  C***KEYWORDS  ERROR, XERROR
  C***AUTHOR  Jones, R. E., (SNLA)
  C***DESCRIPTION
  C
  C *Usage:
  C
  C        INTEGER  KFLAG, NERR, LEVEL, ICOUNT
  C        CHARACTER * (len) LIBRAR, SUBROU, MESSG
  C
  C        CALL XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL, ICOUNT)
  C
  C *Arguments:
  C
  C        LIBRAR :IN    is the library that the message is from.
  C        SUBROU :IN    is the subroutine that the message is from.
  C        MESSG  :IN    is the message to be saved.
  C        KFLAG  :IN    indicates the action to be performed.
  C                      when KFLAG > 0, the message in MESSG is saved.
  C                      when KFLAG=0 the tables will be dumped and
  C                      cleared.
  C                      when KFLAG < 0, the tables will be dumped and
  C                      not cleared.
  C        NERR   :IN    is the error number.
  C        LEVEL  :IN    is the error severity.
  C        ICOUNT :OUT   the number of times this message has been seen,
  C                      or zero if the table has overflowed and does not
  C                      contain this message specifically.  When KFLAG=0,
  C                      ICOUNT will not be altered.
  C
  C *Description:
  C
  C   Record that this error occurred and possibly dump and clear the
  C   tables.
  C
  C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
  C                 Error-handling Package, SAND82-0800, Sandia
  C                 Laboratories, 1982.
  C***ROUTINES CALLED  I1MACH, XGETUA
  C***REVISION HISTORY  (YYMMDD)
  C   800319  DATE WRITTEN
  C   861211  REVISION DATE from Version 3.2
  C   891214  Prologue converted to Version 4.0 format.  (BAB)
  C   900413  Routine modified to remove reference to KFLAG.  (WRB)
  C   900510  Changed to add LIBRARY NAME and SUBROUTINE to calling
  C           sequence, use IF-THEN-ELSE, make number of saved entries
  C           easily changeable, changed routine name from XERSAV to
  C           XERSVE.  (RWC)
  C   910626  Added LIBTAB and SUBTAB to SAVE statement.  (BKS)
  C   920501  Reformatted the REFERENCES section.  (WRB)
  C***END PROLOGUE  XERSVE
        PARAMETER (LENTAB=10)
        INTEGER LUN(5)
        CHARACTER*(*) LIBRAR, SUBROU, MESSG
        CHARACTER*8  LIBTAB(LENTAB), SUBTAB(LENTAB), LIB, SUB
        CHARACTER*20 MESTAB(LENTAB), MES
        DIMENSION NERTAB(LENTAB), LEVTAB(LENTAB), KOUNT(LENTAB)
        SAVE LIBTAB, SUBTAB, MESTAB, NERTAB, LEVTAB, KOUNT, KOUNTX, NMSG
        DATA KOUNTX/0/, NMSG/0/
  C***FIRST EXECUTABLE STATEMENT  XERSVE
  C
        IF (KFLAG.LE.0) THEN
  C
  C        Dump the table.
  C
           IF (NMSG.EQ.0) RETURN
  C
  C        Print to each unit.
  C
           CALL XGETUA (LUN, NUNIT)
           DO 20 KUNIT = 1,NUNIT
              IUNIT = LUN(KUNIT)
              IF (IUNIT.EQ.0) IUNIT = I1MACH(4)
  C
  C           Print the table header.
  C
              WRITE (IUNIT,9000)
  C
  C           Print body of table.
  C
              DO 10 I = 1,NMSG
                 WRITE (IUNIT,9010) LIBTAB(I), SUBTAB(I), MESTAB(I),
       *            NERTAB(I),LEVTAB(I),KOUNT(I)
     10       CONTINUE
  C
  C           Print number of other errors.
  C
              IF (KOUNTX.NE.0) WRITE (IUNIT,9020) KOUNTX
              WRITE (IUNIT,9030)
     20    CONTINUE
  C
  C        Clear the error tables.
  C
           IF (KFLAG.EQ.0) THEN
              NMSG = 0
              KOUNTX = 0
           ENDIF
        ELSE
  C
  C        PROCESS A MESSAGE...
  C        SEARCH FOR THIS MESSG, OR ELSE AN EMPTY SLOT FOR THIS MESSG,
  C        OR ELSE DETERMINE THAT THE ERROR TABLE IS FULL.
  C
           LIB = LIBRAR
           SUB = SUBROU
           MES = MESSG
           DO 30 I = 1,NMSG
              IF (LIB.EQ.LIBTAB(I) .AND. SUB.EQ.SUBTAB(I) .AND.
       *         MES.EQ.MESTAB(I) .AND. NERR.EQ.NERTAB(I) .AND.
       *         LEVEL.EQ.LEVTAB(I)) THEN
                    KOUNT(I) = KOUNT(I) + 1
                    ICOUNT = KOUNT(I)
                    RETURN
              ENDIF
     30    CONTINUE
  C
           IF (NMSG.LT.LENTAB) THEN
  C
  C           Empty slot found for new message.
  C
              NMSG = NMSG + 1
              LIBTAB(I) = LIB
              SUBTAB(I) = SUB
              MESTAB(I) = MES
              NERTAB(I) = NERR
              LEVTAB(I) = LEVEL
              KOUNT (I) = 1
              ICOUNT    = 1
           ELSE
  C
  C           Table is full.
  C
              KOUNTX = KOUNTX+1
              ICOUNT = 0
           ENDIF
        ENDIF
        RETURN
  C
  C     Formats.
  C
   9000 FORMAT ('0          ERROR MESSAGE SUMMARY' /
       +   ' LIBRARY    SUBROUTINE MESSAGE START             NERR',
       +   '     LEVEL     COUNT')
   9010 FORMAT (1X,A,3X,A,3X,A,3I10)
   9020 FORMAT ('0OTHER ERRORS NOT INDIVIDUALLY TABULATED = ', I10)
   9030 FORMAT (1X)
        END
  *DECK D1MACH
        DOUBLE PRECISION FUNCTION D1MACH (I)
  C***BEGIN PROLOGUE  D1MACH
  C***PURPOSE  Return floating point machine dependent constants.
  C***LIBRARY   SLATEC
  C***CATEGORY  R1
  C***TYPE      DOUBLE PRECISION (R1MACH-S, D1MACH-D)
  C***KEYWORDS  MACHINE CONSTANTS
  C***AUTHOR  Fox, P. A., (Bell Labs)
  C           Hall, A. D., (Bell Labs)
  C           Schryer, N. L., (Bell Labs)
  C***DESCRIPTION
  C
  C   D1MACH can be used to obtain machine-dependent parameters for the
  C   local machine environment.  It is a function subprogram with one
  C   (input) argument, and can be referenced as follows:
  C
  C        D = D1MACH(I)
  C
  C   where I=1,...,5.  The (output) value of D above is determined by
  C   the (input) value of I.  The results for various values of I are
  C   discussed below.
  C
  C   D1MACH( 1) = B**(EMIN-1), the smallest positive magnitude.
  C   D1MACH( 2) = B**EMAX*(1 - B**(-T)), the largest magnitude.
  C   D1MACH( 3) = B**(-T), the smallest relative spacing.
  C   D1MACH( 4) = B**(1-T), the largest relative spacing.
  C   D1MACH( 5) = LOG10(B)
  C
  C   Assume double precision numbers are represented in the T-digit,
  C   base-B form
  C
  C              sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
  C
  C   where 0 .LE. X(I) .LT. B for I=1,...,T, 0 .LT. X(1), and
  C   EMIN .LE. E .LE. EMAX.
  C
  C   The values of B, T, EMIN and EMAX are provided in I1MACH as
  C   follows:
  C   I1MACH(10) = B, the base.
  C   I1MACH(14) = T, the number of base-B digits.
  C   I1MACH(15) = EMIN, the smallest exponent E.
  C   I1MACH(16) = EMAX, the largest exponent E.
  C
  C   To alter this function for a particular environment, the desired
  C   set of DATA statements should be activated by removing the C from
  C   column 1.  Also, the values of D1MACH(1) - D1MACH(4) should be
  C   checked for consistency with the local operating system.
  C
  C***REFERENCES  P. A. Fox, A. D. Hall and N. L. Schryer, Framework for
  C                 a portable library, ACM Transactions on Mathematical
  C                 Software 4, 2 (June 1978), pp. 177-188.
  C***ROUTINES CALLED  XERMSG
  C***REVISION HISTORY  (YYMMDD)
  C   750101  DATE WRITTEN
  C   890213  REVISION DATE from Version 3.2
  C   891214  Prologue converted to Version 4.0 format.  (BAB)
  C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
  C   900618  Added DEC RISC constants.  (WRB)
  C   900723  Added IBM RS 6000 constants.  (WRB)
  C   900911  Added SUN 386i constants.  (WRB)
  C   910710  Added HP 730 constants.  (SMR)
  C   911114  Added Convex IEEE constants.  (WRB)
  C   920121  Added SUN -r8 compiler option constants.  (WRB)
  C   920229  Added Touchstone Delta i860 constants.  (WRB)
  C   920501  Reformatted the REFERENCES section.  (WRB)
  C   920625  Added CONVEX -p8 and -pd8 compiler option constants.
  C           (BKS, WRB)
  C   930201  Added DEC Alpha and SGI constants.  (RWC and WRB)
  C   010817  Elevated IEEE to highest importance; see next set of
  C           comments below.  (DWL)
  C***END PROLOGUE  D1MACH
  C
        INTEGER SMALL(4)
        INTEGER LARGE(4)
        INTEGER RIGHT(4)
        INTEGER DIVER(4)
        INTEGER LOG10(4)
  C
  C Initial data here correspond to the IEEE standard.  The values for
  C DMACH(1), DMACH(3) and DMACH(4) are slight upper bounds.  The value
  C for DMACH(2) is a slight lower bound.  The value for DMACH(5) is
  C a 20-digit approximation.  If one of the sets of initial data below
  C is preferred, do the necessary commenting and uncommenting. (DWL)
        DOUBLE PRECISION DMACH(5)
        DATA DMACH / 2.23D-308, 1.79D+308, 1.111D-16, 2.222D-16,
       1             0.30102999566398119521D0 /
        SAVE DMACH
  C
  cc      EQUIVALENCE (DMACH(1),SMALL(1))
  cc      EQUIVALENCE (DMACH(2),LARGE(1))
  cc      EQUIVALENCE (DMACH(3),RIGHT(1))
  cc      EQUIVALENCE (DMACH(4),DIVER(1))
  cc      EQUIVALENCE (DMACH(5),LOG10(1))
  C
  C     MACHINE CONSTANTS FOR THE AMIGA
  C     ABSOFT FORTRAN COMPILER USING THE 68020/68881 COMPILER OPTION
  C
  C     DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' /
  C     DATA LARGE(1), LARGE(2) / Z'7FEFFFFF', Z'FFFFFFFF' /
  C     DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' /
  C     DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' /
  C     DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' /
  C
  C     MACHINE CONSTANTS FOR THE AMIGA
  C     ABSOFT FORTRAN COMPILER USING SOFTWARE FLOATING POINT
  C
  C     DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' /
  C     DATA LARGE(1), LARGE(2) / Z'7FDFFFFF', Z'FFFFFFFF' /
  C     DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' /
  C     DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' /
  C     DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' /
  C
  C     MACHINE CONSTANTS FOR THE APOLLO
  C
  C     DATA SMALL(1), SMALL(2) / 16#00100000, 16#00000000 /
  C     DATA LARGE(1), LARGE(2) / 16#7FFFFFFF, 16#FFFFFFFF /
  C     DATA RIGHT(1), RIGHT(2) / 16#3CA00000, 16#00000000 /
  C     DATA DIVER(1), DIVER(2) / 16#3CB00000, 16#00000000 /
  C     DATA LOG10(1), LOG10(2) / 16#3FD34413, 16#509F79FF /
  C
  C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM
  C
  C     DATA SMALL(1) / ZC00800000 /
  C     DATA SMALL(2) / Z000000000 /
  C     DATA LARGE(1) / ZDFFFFFFFF /
  C     DATA LARGE(2) / ZFFFFFFFFF /
  C     DATA RIGHT(1) / ZCC5800000 /
  C     DATA RIGHT(2) / Z000000000 /
  C     DATA DIVER(1) / ZCC6800000 /
  C     DATA DIVER(2) / Z000000000 /
  C     DATA LOG10(1) / ZD00E730E7 /
  C     DATA LOG10(2) / ZC77800DC0 /
  C
  C     MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM
  C
  C     DATA SMALL(1) / O1771000000000000 /
  C     DATA SMALL(2) / O0000000000000000 /
  C     DATA LARGE(1) / O0777777777777777 /
  C     DATA LARGE(2) / O0007777777777777 /
  C     DATA RIGHT(1) / O1461000000000000 /
  C     DATA RIGHT(2) / O0000000000000000 /
  C     DATA DIVER(1) / O1451000000000000 /
  C     DATA DIVER(2) / O0000000000000000 /
  C     DATA LOG10(1) / O1157163034761674 /
  C     DATA LOG10(2) / O0006677466732724 /
  C
  C     MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS
  C
  C     DATA SMALL(1) / O1771000000000000 /
  C     DATA SMALL(2) / O7770000000000000 /
  C     DATA LARGE(1) / O0777777777777777 /
  C     DATA LARGE(2) / O7777777777777777 /
  C     DATA RIGHT(1) / O1461000000000000 /
  C     DATA RIGHT(2) / O0000000000000000 /
  C     DATA DIVER(1) / O1451000000000000 /
  C     DATA DIVER(2) / O0000000000000000 /
  C     DATA LOG10(1) / O1157163034761674 /
  C     DATA LOG10(2) / O0006677466732724 /
  C
  C     MACHINE CONSTANTS FOR THE CDC 170/180 SERIES USING NOS/VE
  C
  C     DATA SMALL(1) / Z"3001800000000000" /
  C     DATA SMALL(2) / Z"3001000000000000" /
  C     DATA LARGE(1) / Z"4FFEFFFFFFFFFFFE" /
  C     DATA LARGE(2) / Z"4FFE000000000000" /
  C     DATA RIGHT(1) / Z"3FD2800000000000" /
  C     DATA RIGHT(2) / Z"3FD2000000000000" /
  C     DATA DIVER(1) / Z"3FD3800000000000" /
  C     DATA DIVER(2) / Z"3FD3000000000000" /
  C     DATA LOG10(1) / Z"3FFF9A209A84FBCF" /
  C     DATA LOG10(2) / Z"3FFFF7988F8959AC" /
  C
  C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES
  C
  C     DATA SMALL(1) / 00564000000000000000B /
  C     DATA SMALL(2) / 00000000000000000000B /
  C     DATA LARGE(1) / 37757777777777777777B /
  C     DATA LARGE(2) / 37157777777777777777B /
  C     DATA RIGHT(1) / 15624000000000000000B /
  C     DATA RIGHT(2) / 00000000000000000000B /
  C     DATA DIVER(1) / 15634000000000000000B /
  C     DATA DIVER(2) / 00000000000000000000B /
  C     DATA LOG10(1) / 17164642023241175717B /
  C     DATA LOG10(2) / 16367571421742254654B /
  C
  C     MACHINE CONSTANTS FOR THE CELERITY C1260
  C
  C     DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' /
  C     DATA LARGE(1), LARGE(2) / Z'7FEFFFFF', Z'FFFFFFFF' /
  C     DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' /
  C     DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' /
  C     DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' /
  C
  C     MACHINE CONSTANTS FOR THE CONVEX
  C     USING THE -fn OR -pd8 COMPILER OPTION
  C
  C     DATA DMACH(1) / Z'0010000000000000' /
  C     DATA DMACH(2) / Z'7FFFFFFFFFFFFFFF' /
  C     DATA DMACH(3) / Z'3CC0000000000000' /
  C     DATA DMACH(4) / Z'3CD0000000000000' /
  C     DATA DMACH(5) / Z'3FF34413509F79FF' /
  C
  C     MACHINE CONSTANTS FOR THE CONVEX
  C     USING THE -fi COMPILER OPTION
  C
  C     DATA DMACH(1) / Z'0010000000000000' /
  C     DATA DMACH(2) / Z'7FEFFFFFFFFFFFFF' /
  C     DATA DMACH(3) / Z'3CA0000000000000' /
  C     DATA DMACH(4) / Z'3CB0000000000000' /
  C     DATA DMACH(5) / Z'3FD34413509F79FF' /
  C
  C     MACHINE CONSTANTS FOR THE CONVEX
  C     USING THE -p8 COMPILER OPTION
  C
  C     DATA DMACH(1) / Z'00010000000000000000000000000000' /
  C     DATA DMACH(2) / Z'7FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF' /
  C     DATA DMACH(3) / Z'3F900000000000000000000000000000' /
  C     DATA DMACH(4) / Z'3F910000000000000000000000000000' /
  C     DATA DMACH(5) / Z'3FFF34413509F79FEF311F12B35816F9' /
  C
  C     MACHINE CONSTANTS FOR THE CRAY
  C
  C     DATA SMALL(1) / 201354000000000000000B /
  C     DATA SMALL(2) / 000000000000000000000B /
  C     DATA LARGE(1) / 577767777777777777777B /
  C     DATA LARGE(2) / 000007777777777777774B /
  C     DATA RIGHT(1) / 376434000000000000000B /
  C     DATA RIGHT(2) / 000000000000000000000B /
  C     DATA DIVER(1) / 376444000000000000000B /
  C     DATA DIVER(2) / 000000000000000000000B /
  C     DATA LOG10(1) / 377774642023241175717B /
  C     DATA LOG10(2) / 000007571421742254654B /
  C
  C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
  C     NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD -
  C     STATIC DMACH(5)
  C
  C     DATA SMALL /    20K, 3*0 /
  C     DATA LARGE / 77777K, 3*177777K /
  C     DATA RIGHT / 31420K, 3*0 /
  C     DATA DIVER / 32020K, 3*0 /
  C     DATA LOG10 / 40423K, 42023K, 50237K, 74776K /
  C
  C     MACHINE CONSTANTS FOR THE DEC ALPHA
  C     USING G_FLOAT
  C
  C     DATA DMACH(1) / '0000000000000010'X /
  C     DATA DMACH(2) / 'FFFFFFFFFFFF7FFF'X /
  C     DATA DMACH(3) / '0000000000003CC0'X /
  C     DATA DMACH(4) / '0000000000003CD0'X /
  C     DATA DMACH(5) / '79FF509F44133FF3'X /
  C
  C     MACHINE CONSTANTS FOR THE DEC ALPHA
  C     USING IEEE_FORMAT
  C
  C     DATA DMACH(1) / '0010000000000000'X /
  C     DATA DMACH(2) / '7FEFFFFFFFFFFFFF'X /
  C     DATA DMACH(3) / '3CA0000000000000'X /
  C     DATA DMACH(4) / '3CB0000000000000'X /
  C     DATA DMACH(5) / '3FD34413509F79FF'X /
  C
  C     MACHINE CONSTANTS FOR THE DEC RISC
  C
  C     DATA SMALL(1), SMALL(2) / Z'00000000', Z'00100000'/
  C     DATA LARGE(1), LARGE(2) / Z'FFFFFFFF', Z'7FEFFFFF'/
  C     DATA RIGHT(1), RIGHT(2) / Z'00000000', Z'3CA00000'/
  C     DATA DIVER(1), DIVER(2) / Z'00000000', Z'3CB00000'/
  C     DATA LOG10(1), LOG10(2) / Z'509F79FF', Z'3FD34413'/
  C
  C     MACHINE CONSTANTS FOR THE DEC VAX
  C     USING D_FLOATING
  C     (EXPRESSED IN INTEGER AND HEXADECIMAL)
  C     THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS
  C     THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS
  C
  C     DATA SMALL(1), SMALL(2) /        128,           0 /
  C     DATA LARGE(1), LARGE(2) /     -32769,          -1 /
  C     DATA RIGHT(1), RIGHT(2) /       9344,           0 /
  C     DATA DIVER(1), DIVER(2) /       9472,           0 /
  C     DATA LOG10(1), LOG10(2) /  546979738,  -805796613 /
  C
  C     DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 /
  C     DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /
  C     DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 /
  C     DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 /
  C     DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFF884FB /
  C
  C     MACHINE CONSTANTS FOR THE DEC VAX
  C     USING G_FLOATING
  C     (EXPRESSED IN INTEGER AND HEXADECIMAL)
  C     THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS
  C     THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS
  C
  C     DATA SMALL(1), SMALL(2) /         16,           0 /
  C     DATA LARGE(1), LARGE(2) /     -32769,          -1 /
  C     DATA RIGHT(1), RIGHT(2) /      15552,           0 /
  C     DATA DIVER(1), DIVER(2) /      15568,           0 /
  C     DATA LOG10(1), LOG10(2) /  1142112243, 2046775455 /
  C
  C     DATA SMALL(1), SMALL(2) / Z00000010, Z00000000 /
  C     DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF /
  C     DATA RIGHT(1), RIGHT(2) / Z00003CC0, Z00000000 /
  C     DATA DIVER(1), DIVER(2) / Z00003CD0, Z00000000 /
  C     DATA LOG10(1), LOG10(2) / Z44133FF3, Z79FF509F /
  C
  C     MACHINE CONSTANTS FOR THE ELXSI 6400
  C     (ASSUMING REAL*8 IS THE DEFAULT DOUBLE PRECISION)
  C
  C     DATA SMALL(1), SMALL(2) / '00100000'X,'00000000'X /
  C     DATA LARGE(1), LARGE(2) / '7FEFFFFF'X,'FFFFFFFF'X /
  C     DATA RIGHT(1), RIGHT(2) / '3CB00000'X,'00000000'X /
  C     DATA DIVER(1), DIVER(2) / '3CC00000'X,'00000000'X /
  C     DATA LOG10(1), LOG10(2) / '3FD34413'X,'509F79FF'X /
  C
  C     MACHINE CONSTANTS FOR THE HARRIS 220
  C
  C     DATA SMALL(1), SMALL(2) / '20000000, '00000201 /
  C     DATA LARGE(1), LARGE(2) / '37777777, '37777577 /
  C     DATA RIGHT(1), RIGHT(2) / '20000000, '00000333 /
  C     DATA DIVER(1), DIVER(2) / '20000000, '00000334 /
  C     DATA LOG10(1), LOG10(2) / '23210115, '10237777 /
  C
  C     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES
  C
  C     DATA SMALL(1), SMALL(2) / O402400000000, O000000000000 /
  C     DATA LARGE(1), LARGE(2) / O376777777777, O777777777777 /
  C     DATA RIGHT(1), RIGHT(2) / O604400000000, O000000000000 /
  C     DATA DIVER(1), DIVER(2) / O606400000000, O000000000000 /
  C     DATA LOG10(1), LOG10(2) / O776464202324, O117571775714 /
  C
  C     MACHINE CONSTANTS FOR THE HP 730
  C
  C     DATA DMACH(1) / Z'0010000000000000' /
  C     DATA DMACH(2) / Z'7FEFFFFFFFFFFFFF' /
  C     DATA DMACH(3) / Z'3CA0000000000000' /
  C     DATA DMACH(4) / Z'3CB0000000000000' /
  C     DATA DMACH(5) / Z'3FD34413509F79FF' /
  C
  C     MACHINE CONSTANTS FOR THE HP 2100
  C     THREE WORD DOUBLE PRECISION OPTION WITH FTN4
  C
  C     DATA SMALL(1), SMALL(2), SMALL(3) / 40000B,       0,       1 /
  C     DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B /
  C     DATA RIGHT(1), RIGHT(2), RIGHT(3) / 40000B,       0,    265B /
  C     DATA DIVER(1), DIVER(2), DIVER(3) / 40000B,       0,    276B /
  C     DATA LOG10(1), LOG10(2), LOG10(3) / 46420B,  46502B,  77777B /
  C
  C     MACHINE CONSTANTS FOR THE HP 2100
  C     FOUR WORD DOUBLE PRECISION OPTION WITH FTN4
  C
  C     DATA SMALL(1), SMALL(2) /  40000B,       0 /
  C     DATA SMALL(3), SMALL(4) /       0,       1 /
  C     DATA LARGE(1), LARGE(2) /  77777B, 177777B /
  C     DATA LARGE(3), LARGE(4) / 177777B, 177776B /
  C     DATA RIGHT(1), RIGHT(2) /  40000B,       0 /
  C     DATA RIGHT(3), RIGHT(4) /       0,    225B /
  C     DATA DIVER(1), DIVER(2) /  40000B,       0 /
  C     DATA DIVER(3), DIVER(4) /       0,    227B /
  C     DATA LOG10(1), LOG10(2) /  46420B,  46502B /
  C     DATA LOG10(3), LOG10(4) /  76747B, 176377B /
  C
  C     MACHINE CONSTANTS FOR THE HP 9000
  C
  C     DATA SMALL(1), SMALL(2) / 00040000000B, 00000000000B /
  C     DATA LARGE(1), LARGE(2) / 17737777777B, 37777777777B /
  C     DATA RIGHT(1), RIGHT(2) / 07454000000B, 00000000000B /
  C     DATA DIVER(1), DIVER(2) / 07460000000B, 00000000000B /
  C     DATA LOG10(1), LOG10(2) / 07764642023B, 12047674777B /
  C
  C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
  C     THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND
  C     THE PERKIN ELMER (INTERDATA) 7/32.
  C
  C     DATA SMALL(1), SMALL(2) / Z00100000, Z00000000 /
  C     DATA LARGE(1), LARGE(2) / Z7FFFFFFF, ZFFFFFFFF /
  C     DATA RIGHT(1), RIGHT(2) / Z33100000, Z00000000 /
  C     DATA DIVER(1), DIVER(2) / Z34100000, Z00000000 /
  C     DATA LOG10(1), LOG10(2) / Z41134413, Z509F79FF /
  C
  C     MACHINE CONSTANTS FOR THE IBM PC
  C     ASSUMES THAT ALL ARITHMETIC IS DONE IN DOUBLE PRECISION
  C     ON 8088, I.E., NOT IN 80 BIT FORM FOR THE 8087.
  C
  C     DATA SMALL(1) / 2.23D-308  /
  C     DATA LARGE(1) / 1.79D+308  /
  C     DATA RIGHT(1) / 1.11D-16   /
  C     DATA DIVER(1) / 2.22D-16   /
  C     DATA LOG10(1) / 0.301029995663981195D0 /
  C
  C     MACHINE CONSTANTS FOR THE IBM RS 6000
  C
  C     DATA DMACH(1) / Z'0010000000000000' /
  C     DATA DMACH(2) / Z'7FEFFFFFFFFFFFFF' /
  C     DATA DMACH(3) / Z'3CA0000000000000' /
  C     DATA DMACH(4) / Z'3CB0000000000000' /
  C     DATA DMACH(5) / Z'3FD34413509F79FF' /
  C
  C     MACHINE CONSTANTS FOR THE INTEL i860
  C
  C     DATA DMACH(1) / Z'0010000000000000' /
  C     DATA DMACH(2) / Z'7FEFFFFFFFFFFFFF' /
  C     DATA DMACH(3) / Z'3CA0000000000000' /
  C     DATA DMACH(4) / Z'3CB0000000000000' /
  C     DATA DMACH(5) / Z'3FD34413509F79FF' /
  C
  C     MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR)
  C
  C     DATA SMALL(1), SMALL(2) / "033400000000, "000000000000 /
  C     DATA LARGE(1), LARGE(2) / "377777777777, "344777777777 /
  C     DATA RIGHT(1), RIGHT(2) / "113400000000, "000000000000 /
  C     DATA DIVER(1), DIVER(2) / "114400000000, "000000000000 /
  C     DATA LOG10(1), LOG10(2) / "177464202324, "144117571776 /
  C
  C     MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR)
  C
  C     DATA SMALL(1), SMALL(2) / "000400000000, "000000000000 /
  C     DATA LARGE(1), LARGE(2) / "377777777777, "377777777777 /
  C     DATA RIGHT(1), RIGHT(2) / "103400000000, "000000000000 /
  C     DATA DIVER(1), DIVER(2) / "104400000000, "000000000000 /
  C     DATA LOG10(1), LOG10(2) / "177464202324, "476747767461 /
  C
  C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
  C     32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
  C
  C     DATA SMALL(1), SMALL(2) /    8388608,           0 /
  C     DATA LARGE(1), LARGE(2) / 2147483647,          -1 /
  C     DATA RIGHT(1), RIGHT(2) /  612368384,           0 /
  C     DATA DIVER(1), DIVER(2) /  620756992,           0 /
  C     DATA LOG10(1), LOG10(2) / 1067065498, -2063872008 /
  C
  C     DATA SMALL(1), SMALL(2) / O00040000000, O00000000000 /
  C     DATA LARGE(1), LARGE(2) / O17777777777, O37777777777 /
  C     DATA RIGHT(1), RIGHT(2) / O04440000000, O00000000000 /
  C     DATA DIVER(1), DIVER(2) / O04500000000, O00000000000 /
  C     DATA LOG10(1), LOG10(2) / O07746420232, O20476747770 /
  C
  C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
  C     16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL).
  C
  C     DATA SMALL(1), SMALL(2) /    128,      0 /
  C     DATA SMALL(3), SMALL(4) /      0,      0 /
  C     DATA LARGE(1), LARGE(2) /  32767,     -1 /
  C     DATA LARGE(3), LARGE(4) /     -1,     -1 /
  C     DATA RIGHT(1), RIGHT(2) /   9344,      0 /
  C     DATA RIGHT(3), RIGHT(4) /      0,      0 /
  C     DATA DIVER(1), DIVER(2) /   9472,      0 /
  C     DATA DIVER(3), DIVER(4) /      0,      0 /
  C     DATA LOG10(1), LOG10(2) /  16282,   8346 /
  C     DATA LOG10(3), LOG10(4) / -31493, -12296 /
  C
  C     DATA SMALL(1), SMALL(2) / O000200, O000000 /
  C     DATA SMALL(3), SMALL(4) / O000000, O000000 /
  C     DATA LARGE(1), LARGE(2) / O077777, O177777 /
  C     DATA LARGE(3), LARGE(4) / O177777, O177777 /
  C     DATA RIGHT(1), RIGHT(2) / O022200, O000000 /
  C     DATA RIGHT(3), RIGHT(4) / O000000, O000000 /
  C     DATA DIVER(1), DIVER(2) / O022400, O000000 /
  C     DATA DIVER(3), DIVER(4) / O000000, O000000 /
  C     DATA LOG10(1), LOG10(2) / O037632, O020232 /
  C     DATA LOG10(3), LOG10(4) / O102373, O147770 /
  C
  C     MACHINE CONSTANTS FOR THE SILICON GRAPHICS
  C
  C     DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' /
  C     DATA LARGE(1), LARGE(2) / Z'7FEFFFFF', Z'FFFFFFFF' /
  C     DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' /
  C     DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' /
  C     DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' /
  C
  C     MACHINE CONSTANTS FOR THE SUN
  C
  C     DATA DMACH(1) / Z'0010000000000000' /
  C     DATA DMACH(2) / Z'7FEFFFFFFFFFFFFF' /
  C     DATA DMACH(3) / Z'3CA0000000000000' /
  C     DATA DMACH(4) / Z'3CB0000000000000' /
  C     DATA DMACH(5) / Z'3FD34413509F79FF' /
  C
  C     MACHINE CONSTANTS FOR THE SUN
  C     USING THE -r8 COMPILER OPTION
  C
  C     DATA DMACH(1) / Z'00010000000000000000000000000000' /
  C     DATA DMACH(2) / Z'7FFEFFFFFFFFFFFFFFFFFFFFFFFFFFFF' /
  C     DATA DMACH(3) / Z'3F8E0000000000000000000000000000' /
  C     DATA DMACH(4) / Z'3F8F0000000000000000000000000000' /
  C     DATA DMACH(5) / Z'3FFD34413509F79FEF311F12B35816F9' /
  C
  C     MACHINE CONSTANTS FOR THE SUN 386i
  C
  C     DATA SMALL(1), SMALL(2) / Z'FFFFFFFD', Z'000FFFFF' /
  C     DATA LARGE(1), LARGE(2) / Z'FFFFFFB0', Z'7FEFFFFF' /
  C     DATA RIGHT(1), RIGHT(2) / Z'000000B0', Z'3CA00000' /
  C     DATA DIVER(1), DIVER(2) / Z'FFFFFFCB', Z'3CAFFFFF'
  C     DATA LOG10(1), LOG10(2) / Z'509F79E9', Z'3FD34413' /
  C
  C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES FTN COMPILER
  C
  C     DATA SMALL(1), SMALL(2) / O000040000000, O000000000000 /
  C     DATA LARGE(1), LARGE(2) / O377777777777, O777777777777 /
  C     DATA RIGHT(1), RIGHT(2) / O170540000000, O000000000000 /
  C     DATA DIVER(1), DIVER(2) / O170640000000, O000000000000 /
  C     DATA LOG10(1), LOG10(2) / O177746420232, O411757177572 /
  C
  C***FIRST EXECUTABLE STATEMENT  D1MACH
        IF (I .LT. 1 .OR. I .GT. 5) CALL XERMSG ('SLATEC', 'D1MACH',
       +   'I OUT OF BOUNDS', 1, 2)
  C
        D1MACH = DMACH(I)
        RETURN
  C
        END
  *DECK XGETUA
        SUBROUTINE XGETUA (IUNITA, N)
  C***BEGIN PROLOGUE  XGETUA
  C***PURPOSE  Return unit number(s) to which error messages are being
  C            sent.
  C***LIBRARY   SLATEC (XERROR)
  C***CATEGORY  R3C
  C***TYPE      ALL (XGETUA-A)
  C***KEYWORDS  ERROR, XERROR
  C***AUTHOR  Jones, R. E., (SNLA)
  C***DESCRIPTION
  C
  C     Abstract
  C        XGETUA may be called to determine the unit number or numbers
  C        to which error messages are being sent.
  C        These unit numbers may have been set by a call to XSETUN,
  C        or a call to XSETUA, or may be a default value.
  C
  C     Description of Parameters
  C      --Output--
  C        IUNIT - an array of one to five unit numbers, depending
  C                on the value of N.  A value of zero refers to the
  C                default unit, as defined by the I1MACH machine
  C                constant routine.  Only IUNIT(1),...,IUNIT(N) are
  C                defined by XGETUA.  The values of IUNIT(N+1),...,
  C                IUNIT(5) are not defined (for N .LT. 5) or altered
  C                in any way by XGETUA.
  C        N     - the number of units to which copies of the
  C                error messages are being sent.  N will be in the
  C                range from 1 to 5.
  C
  C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
  C                 Error-handling Package, SAND82-0800, Sandia
  C                 Laboratories, 1982.
  C***ROUTINES CALLED  J4SAVE
  C***REVISION HISTORY  (YYMMDD)
  C   790801  DATE WRITTEN
  C   861211  REVISION DATE from Version 3.2
  C   891214  Prologue converted to Version 4.0 format.  (BAB)
  C   920501  Reformatted the REFERENCES section.  (WRB)
  C***END PROLOGUE  XGETUA
        DIMENSION IUNITA(5)
  C***FIRST EXECUTABLE STATEMENT  XGETUA
        N = J4SAVE(5,0,.FALSE.)
        DO 30 I=1,N
           INDEX = I+4
           IF (I.EQ.1) INDEX = 3
           IUNITA(I) = J4SAVE(INDEX,0,.FALSE.)
     30 CONTINUE
        RETURN
        END
  *DECK DSTEPS
        SUBROUTINE DSTEPS (DF, NEQN, Y, X, H, EPS, WT, START, HOLD, K,
       +   KOLD, CRASH, PHI, P, YP, PSI, ALPHA, BETA, SIG, V, W, G,
       +   PHASE1, NS, NORND, KSTEPS, TWOU, FOURU, XOLD, KPREV, IVC, IV,
       +   KGI, GI, RPAR, IPAR)
  C***BEGIN PROLOGUE  DSTEPS
  C***PURPOSE  Integrate a system of first order ordinary differential
  C            equations one step.
  C***LIBRARY   SLATEC (DEPAC)
  C***CATEGORY  I1A1B
  C***TYPE      DOUBLE PRECISION (STEPS-S, DSTEPS-D)
  C***KEYWORDS  ADAMS METHOD, DEPAC, INITIAL VALUE PROBLEMS, ODE,
  C             ORDINARY DIFFERENTIAL EQUATIONS, PREDICTOR-CORRECTOR
  C***AUTHOR  Shampine, L. F., (SNLA)
  C           Gordon, M. K., (SNLA)
  C             MODIFIED BY H.A. WATTS
  C***DESCRIPTION
  C
  C   Written by L. F. Shampine and M. K. Gordon
  C
  C   Abstract
  C
  C   Subroutine  DSTEPS  is normally used indirectly through subroutine
  C   DDEABM .  Because  DDEABM  suffices for most problems and is much
  C   easier to use, using it should be considered before using  DSTEPS
  C   alone.
  C
  C   Subroutine DSTEPS integrates a system of  NEQN  first order ordinary
  C   differential equations one step, normally from X to X+H, using a
  C   modified divided difference form of the Adams Pece formulas.  Local
  C   extrapolation is used to improve absolute stability and accuracy.
  C   The code adjusts its order and step size to control the local error
  C   per unit step in a generalized sense.  Special devices are included
  C   to control roundoff error and to detect when the user is requesting
  C   too much accuracy.
  C
  C   This code is completely explained and documented in the text,
  C   Computer Solution of Ordinary Differential Equations, The Initial
  C   Value Problem  by L. F. Shampine and M. K. Gordon.
  C   Further details on use of this code are available in "Solving
  C   Ordinary Differential Equations with ODE, STEP, and INTRP",
  C   by L. F. Shampine and M. K. Gordon, SLA-73-1060.
  C
  C
  C   The parameters represent --
  C      DF -- subroutine to evaluate derivatives
  C      NEQN -- number of equations to be integrated
  C      Y(*) -- solution vector at X
  C      X -- independent variable
  C      H -- appropriate step size for next step.  Normally determined by
  C           code
  C      EPS -- local error tolerance
  C      WT(*) -- vector of weights for error criterion
  C      START -- logical variable set .TRUE. for first step,  .FALSE.
  C           otherwise
  C      HOLD -- step size used for last successful step
  C      K -- appropriate order for next step (determined by code)
  C      KOLD -- order used for last successful step
  C      CRASH -- logical variable set .TRUE. when no step can be taken,
  C           .FALSE. otherwise.
  C      YP(*) -- derivative of solution vector at  X  after successful
  C           step
  C      KSTEPS -- counter on attempted steps
  C      TWOU -- 2.*U where U is machine unit roundoff quantity
  C      FOURU -- 4.*U where U is machine unit roundoff quantity
  C      RPAR,IPAR -- parameter arrays which you may choose to use
  C            for communication between your program and subroutine F.
  C            They are not altered or used by DSTEPS.
  C   The variables X,XOLD,KOLD,KGI and IVC and the arrays Y,PHI,ALPHA,G,
  C   W,P,IV and GI are required for the interpolation subroutine SINTRP.
  C   The remaining variables and arrays are included in the call list
  C   only to eliminate local retention of variables between calls.
  C
  C   Input to DSTEPS
  C
  C      First call --
  C
  C   The user must provide storage in his calling program for all arrays
  C   in the call list, namely
  C
  C     DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),PSI(12),
  C    1  ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13),GI(11),IV(10),
  C    2  RPAR(*),IPAR(*)
  C
  C    **Note**
  C
  C   The user must also declare  START ,  CRASH ,  PHASE1  and  NORND
  C   logical variables and  DF  an EXTERNAL subroutine, supply the
  C   subroutine  DF(X,Y,YP)  to evaluate
  C      DY(I)/DX = YP(I) = DF(X,Y(1),Y(2),...,Y(NEQN))
  C   and initialize only the following parameters.
  C      NEQN -- number of equations to be integrated
  C      Y(*) -- vector of initial values of dependent variables
  C      X -- initial value of the independent variable
  C      H -- nominal step size indicating direction of integration
  C           and maximum size of step.  Must be variable
  C      EPS -- local error tolerance per step.  Must be variable
  C      WT(*) -- vector of non-zero weights for error criterion
  C      START -- .TRUE.
  C      YP(*) -- vector of initial derivative values
  C      KSTEPS -- set KSTEPS to zero
  C      TWOU -- 2.*U where U is machine unit roundoff quantity
  C      FOURU -- 4.*U where U is machine unit roundoff quantity
  C   Define U to be the machine unit roundoff quantity by calling
  C   the function routine  D1MACH,  U = D1MACH(4), or by
  C   computing U so that U is the smallest positive number such
  C   that 1.0+U .GT. 1.0.
  C
  C   DSTEPS  requires that the L2 norm of the vector with components
  C   LOCAL ERROR(L)/WT(L)  be less than  EPS  for a successful step.  The
  C   array  WT  allows the user to specify an error test appropriate
  C   for his problem.  For example,
  C      WT(L) = 1.0  specifies absolute error,
  C            = ABS(Y(L))  error relative to the most recent value of the
  C                 L-th component of the solution,
  C            = ABS(YP(L))  error relative to the most recent value of
  C                 the L-th component of the derivative,
  C            = MAX(WT(L),ABS(Y(L)))  error relative to the largest
  C                 magnitude of L-th component obtained so far,
  C            = ABS(Y(L))*RELERR/EPS + ABSERR/EPS  specifies a mixed
  C                 relative-absolute test where  RELERR  is relative
  C                 error,  ABSERR  is absolute error and  EPS =
  C                 MAX(RELERR,ABSERR) .
  C
  C      Subsequent calls --
  C
  C   Subroutine  DSTEPS  is designed so that all information needed to
  C   continue the integration, including the step size  H  and the order
  C   K , is returned with each step.  With the exception of the step
  C   size, the error tolerance, and the weights, none of the parameters
  C   should be altered.  The array  WT  must be updated after each step
  C   to maintain relative error tests like those above.  Normally the
  C   integration is continued just beyond the desired endpoint and the
  C   solution interpolated there with subroutine  SINTRP .  If it is
  C   impossible to integrate beyond the endpoint, the step size may be
  C   reduced to hit the endpoint since the code will not take a step
  C   larger than the  H  input.  Changing the direction of integration,
  C   i.e., the sign of  H , requires the user set  START = .TRUE. before
  C   calling  DSTEPS  again.  This is the only situation in which  START
  C   should be altered.
  C
  C   Output from DSTEPS
  C
  C      Successful Step --
  C
  C   The subroutine returns after each successful step with  START  and
  C   CRASH  set .FALSE. .  X  represents the independent variable
  C   advanced one step of length  HOLD  from its value on input and  Y
  C   the solution vector at the new value of  X .  All other parameters
  C   represent information corresponding to the new  X  needed to
  C   continue the integration.
  C
  C      Unsuccessful Step --
  C
  C   When the error tolerance is too small for the machine precision,
  C   the subroutine returns without taking a step and  CRASH = .TRUE. .
  C   An appropriate step size and error tolerance for continuing are
  C   estimated and all other information is restored as upon input
  C   before returning.  To continue with the larger tolerance, the user
  C   just calls the code again.  A restart is neither required nor
  C   desirable.
  C
  C***REFERENCES  L. F. Shampine and M. K. Gordon, Solving ordinary
  C                 differential equations with ODE, STEP, and INTRP,
  C                 Report SLA-73-1060, Sandia Laboratories, 1973.
  C***ROUTINES CALLED  D1MACH, DHSTRT
  C***REVISION HISTORY  (YYMMDD)
  C   740101  DATE WRITTEN
  C   890531  Changed all specific intrinsics to generic.  (WRB)
  C   890831  Modified array declarations.  (WRB)
  C   890831  REVISION DATE from Version 3.2
  C   891214  Prologue converted to Version 4.0 format.  (BAB)
  C   920501  Reformatted the REFERENCES section.  (WRB)
  C***END PROLOGUE  DSTEPS
  C
        INTEGER I, IFAIL, IM1, IP1, IPAR, IQ, J, K, KM1, KM2, KNEW,
       1      KOLD, KP1, KP2, KSTEPS, L, LIMIT1, LIMIT2, NEQN, NS, NSM2,
       2      NSP1, NSP2
        DOUBLE PRECISION ABSH, ALPHA, BETA, BIG, D1MACH,
       1      EPS, ERK, ERKM1, ERKM2, ERKP1, ERR,
       2      FOURU, G, GI, GSTR, H, HNEW, HOLD, P, P5EPS, PHI, PSI, R,
       3      REALI, REALNS, RHO, ROUND, RPAR, SIG, TAU, TEMP1,
       4      TEMP2, TEMP3, TEMP4, TEMP5, TEMP6, TWO, TWOU, U, V, W, WT,
       5      X, XOLD, Y, YP
        LOGICAL START,CRASH,PHASE1,NORND
        DIMENSION Y(*),WT(*),PHI(NEQN,16),P(*),YP(*),PSI(12),
       1  ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13),GI(11),IV(10),
       2  RPAR(*),IPAR(*)
        DIMENSION TWO(13),GSTR(13)
        EXTERNAL DF
        SAVE TWO, GSTR
  C
        DATA TWO(1),TWO(2),TWO(3),TWO(4),TWO(5),TWO(6),TWO(7),TWO(8),
       1     TWO(9),TWO(10),TWO(11),TWO(12),TWO(13)
       2     /2.0D0,4.0D0,8.0D0,16.0D0,32.0D0,64.0D0,128.0D0,256.0D0,
       3      512.0D0,1024.0D0,2048.0D0,4096.0D0,8192.0D0/
        DATA GSTR(1),GSTR(2),GSTR(3),GSTR(4),GSTR(5),GSTR(6),GSTR(7),
       1     GSTR(8),GSTR(9),GSTR(10),GSTR(11),GSTR(12),GSTR(13)
       2     /0.5D0,0.0833D0,0.0417D0,0.0264D0,0.0188D0,0.0143D0,0.0114D0,
       3      0.00936D0,0.00789D0,0.00679D0,0.00592D0,0.00524D0,0.00468D0/
  C
  C       ***     BEGIN BLOCK 0     ***
  C   CHECK IF STEP SIZE OR ERROR TOLERANCE IS TOO SMALL FOR MACHINE
  C   PRECISION.  IF FIRST STEP, INITIALIZE PHI ARRAY AND ESTIMATE A
  C   STARTING STEP SIZE.
  C                   ***
  C
  C   IF STEP SIZE IS TOO SMALL, DETERMINE AN ACCEPTABLE ONE
  C
  C***FIRST EXECUTABLE STATEMENT  DSTEPS
        CRASH = .TRUE.
        IF(ABS(H) .GE. FOURU*ABS(X)) GO TO 5
        H = SIGN(FOURU*ABS(X),H)
        RETURN
   5    P5EPS = 0.5D0*EPS
  C
  C   IF ERROR TOLERANCE IS TOO SMALL, INCREASE IT TO AN ACCEPTABLE VALUE
  C
        ROUND = 0.0D0
        DO 10 L = 1,NEQN
   10     ROUND = ROUND + (Y(L)/WT(L))**2
        ROUND = TWOU*SQRT(ROUND)
        IF(P5EPS .GE. ROUND) GO TO 15
        EPS = 2.0D0*ROUND*(1.0D0 + FOURU)
        RETURN
   15   CRASH = .FALSE.
        G(1) = 1.0D0
        G(2) = 0.5D0
        SIG(1) = 1.0D0
        IF(.NOT.START) GO TO 99
  C
  C   INITIALIZE.  COMPUTE APPROPRIATE STEP SIZE FOR FIRST STEP
  C
  C     CALL DF(X,Y,YP,RPAR,IPAR)
  C     SUM = 0.0
        DO 20 L = 1,NEQN
          PHI(L,1) = YP(L)
     20   PHI(L,2) = 0.0D0
  C20     SUM = SUM + (YP(L)/WT(L))**2
  C     SUM = SQRT(SUM)
  C     ABSH = ABS(H)
  C     IF(EPS .LT. 16.0*SUM*H*H) ABSH = 0.25*SQRT(EPS/SUM)
  C     H = SIGN(MAX(ABSH,FOURU*ABS(X)),H)
  C
        U = D1MACH(4)
        BIG = SQRT(D1MACH(2))
        CALL DHSTRT(DF,NEQN,X,X+H,Y,YP,WT,1,U,BIG,
       1             PHI(1,3),PHI(1,4),PHI(1,5),PHI(1,6),RPAR,IPAR,H)
  C
        HOLD = 0.0D0
        K = 1
        KOLD = 0
        KPREV = 0
        START = .FALSE.
        PHASE1 = .TRUE.
        NORND = .TRUE.
        IF(P5EPS .GT. 100.0D0*ROUND) GO TO 99
        NORND = .FALSE.
        DO 25 L = 1,NEQN
   25     PHI(L,15) = 0.0D0
   99   IFAIL = 0
  C       ***     END BLOCK 0     ***
  C
  C       ***     BEGIN BLOCK 1     ***
  C   COMPUTE COEFFICIENTS OF FORMULAS FOR THIS STEP.  AVOID COMPUTING
  C   THOSE QUANTITIES NOT CHANGED WHEN STEP SIZE IS NOT CHANGED.
  C                   ***
  C
   100  KP1 = K+1
        KP2 = K+2
        KM1 = K-1
        KM2 = K-2
  C
  C   NS IS THE NUMBER OF DSTEPS TAKEN WITH SIZE H, INCLUDING THE CURRENT
  C   ONE.  WHEN K.LT.NS, NO COEFFICIENTS CHANGE
  C
        IF(H .NE. HOLD) NS = 0
        IF (NS.LE.KOLD) NS = NS+1
        NSP1 = NS+1
        IF (K .LT. NS) GO TO 199
  C
  C   COMPUTE THOSE COMPONENTS OF ALPHA(*),BETA(*),PSI(*),SIG(*) WHICH
  C   ARE CHANGED
  C
        BETA(NS) = 1.0D0
        REALNS = NS
        ALPHA(NS) = 1.0D0/REALNS
        TEMP1 = H*REALNS
        SIG(NSP1) = 1.0D0
        IF(K .LT. NSP1) GO TO 110
        DO 105 I = NSP1,K
          IM1 = I-1
          TEMP2 = PSI(IM1)
          PSI(IM1) = TEMP1
          BETA(I) = BETA(IM1)*PSI(IM1)/TEMP2
          TEMP1 = TEMP2 + H
          ALPHA(I) = H/TEMP1
          REALI = I
   105    SIG(I+1) = REALI*ALPHA(I)*SIG(I)
   110  PSI(K) = TEMP1
  C
  C   COMPUTE COEFFICIENTS G(*)
  C
  C   INITIALIZE V(*) AND SET W(*).
  C
        IF(NS .GT. 1) GO TO 120
        DO 115 IQ = 1,K
          TEMP3 = IQ*(IQ+1)
          V(IQ) = 1.0D0/TEMP3
   115    W(IQ) = V(IQ)
        IVC = 0
        KGI = 0
        IF (K .EQ. 1) GO TO 140
        KGI = 1
        GI(1) = W(2)
        GO TO 140
  C
  C   IF ORDER WAS RAISED, UPDATE DIAGONAL PART OF V(*)
  C
   120  IF(K .LE. KPREV) GO TO 130
        IF (IVC .EQ. 0) GO TO 122
        JV = KP1 - IV(IVC)
        IVC = IVC - 1
        GO TO 123
   122  JV = 1
        TEMP4 = K*KP1
        V(K) = 1.0D0/TEMP4
        W(K) = V(K)
        IF (K .NE. 2) GO TO 123
        KGI = 1
        GI(1) = W(2)
   123  NSM2 = NS-2
        IF(NSM2 .LT. JV) GO TO 130
        DO 125 J = JV,NSM2
          I = K-J
          V(I) = V(I) - ALPHA(J+1)*V(I+1)
   125    W(I) = V(I)
        IF (I .NE. 2) GO TO 130
        KGI = NS - 1
        GI(KGI) = W(2)
  C
  C   UPDATE V(*) AND SET W(*)
  C
   130  LIMIT1 = KP1 - NS
        TEMP5 = ALPHA(NS)
        DO 135 IQ = 1,LIMIT1
          V(IQ) = V(IQ) - TEMP5*V(IQ+1)
   135    W(IQ) = V(IQ)
        G(NSP1) = W(1)
        IF (LIMIT1 .EQ. 1) GO TO 137
        KGI = NS
        GI(KGI) = W(2)
   137  W(LIMIT1+1) = V(LIMIT1+1)
        IF (K .GE. KOLD) GO TO 140
        IVC = IVC + 1
        IV(IVC) = LIMIT1 + 2
  C
  C   COMPUTE THE G(*) IN THE WORK VECTOR W(*)
  C
   140  NSP2 = NS + 2
        KPREV = K
        IF(KP1 .LT. NSP2) GO TO 199
        DO 150 I = NSP2,KP1
          LIMIT2 = KP2 - I
          TEMP6 = ALPHA(I-1)
          DO 145 IQ = 1,LIMIT2
   145      W(IQ) = W(IQ) - TEMP6*W(IQ+1)
   150    G(I) = W(1)
   199    CONTINUE
  C       ***     END BLOCK 1     ***
  C
  C       ***     BEGIN BLOCK 2     ***
  C   PREDICT A SOLUTION P(*), EVALUATE DERIVATIVES USING PREDICTED
  C   SOLUTION, ESTIMATE LOCAL ERROR AT ORDER K AND ERRORS AT ORDERS K,
  C   K-1, K-2 AS IF CONSTANT STEP SIZE WERE USED.
  C                   ***
  C
  C   INCREMENT COUNTER ON ATTEMPTED DSTEPS
  C
        KSTEPS = KSTEPS + 1
  C
  C   CHANGE PHI TO PHI STAR
  C
        IF(K .LT. NSP1) GO TO 215
        DO 210 I = NSP1,K
          TEMP1 = BETA(I)
          DO 205 L = 1,NEQN
   205      PHI(L,I) = TEMP1*PHI(L,I)
   210    CONTINUE
  C
  C   PREDICT SOLUTION AND DIFFERENCES
  C
   215  DO 220 L = 1,NEQN
          PHI(L,KP2) = PHI(L,KP1)
          PHI(L,KP1) = 0.0D0
   220    P(L) = 0.0D0
        DO 230 J = 1,K
          I = KP1 - J
          IP1 = I+1
          TEMP2 = G(I)
          DO 225 L = 1,NEQN
            P(L) = P(L) + TEMP2*PHI(L,I)
   225      PHI(L,I) = PHI(L,I) + PHI(L,IP1)
   230    CONTINUE
        IF(NORND) GO TO 240
        DO 235 L = 1,NEQN
          TAU = H*P(L) - PHI(L,15)
          P(L) = Y(L) + TAU
   235    PHI(L,16) = (P(L) - Y(L)) - TAU
        GO TO 250
   240  DO 245 L = 1,NEQN
   245    P(L) = Y(L) + H*P(L)
   250  XOLD = X
        X = X + H
        ABSH = ABS(H)
        CALL DF(X,P,YP,RPAR,IPAR)
  C
  C   ESTIMATE ERRORS AT ORDERS K,K-1,K-2
  C
        ERKM2 = 0.0D0
        ERKM1 = 0.0D0
        ERK = 0.0D0
        DO 265 L = 1,NEQN
          TEMP3 = 1.0D0/WT(L)
          TEMP4 = YP(L) - PHI(L,1)
          IF(KM2)265,260,255
   255    ERKM2 = ERKM2 + ((PHI(L,KM1)+TEMP4)*TEMP3)**2
   260    ERKM1 = ERKM1 + ((PHI(L,K)+TEMP4)*TEMP3)**2
   265    ERK = ERK + (TEMP4*TEMP3)**2
        IF(KM2)280,275,270
   270  ERKM2 = ABSH*SIG(KM1)*GSTR(KM2)*SQRT(ERKM2)
   275  ERKM1 = ABSH*SIG(K)*GSTR(KM1)*SQRT(ERKM1)
   280  TEMP5 = ABSH*SQRT(ERK)
        ERR = TEMP5*(G(K)-G(KP1))
        ERK = TEMP5*SIG(KP1)*GSTR(K)
        KNEW = K
  C
  C   TEST IF ORDER SHOULD BE LOWERED
  C
        IF(KM2)299,290,285
   285  IF(MAX(ERKM1,ERKM2) .LE. ERK) KNEW = KM1
        GO TO 299
   290  IF(ERKM1 .LE. 0.5D0*ERK) KNEW = KM1
  C
  C   TEST IF STEP SUCCESSFUL
  C
   299  IF(ERR .LE. EPS) GO TO 400
  C       ***     END BLOCK 2     ***
  C
  C       ***     BEGIN BLOCK 3     ***
  C   THE STEP IS UNSUCCESSFUL.  RESTORE  X, PHI(*,*), PSI(*) .
  C   IF THIRD CONSECUTIVE FAILURE, SET ORDER TO ONE.  IF STEP FAILS MORE
  C   THAN THREE TIMES, CONSIDER AN OPTIMAL STEP SIZE.  DOUBLE ERROR
  C   TOLERANCE AND RETURN IF ESTIMATED STEP SIZE IS TOO SMALL FOR MACHINE
  C   PRECISION.
  C                   ***
  C
  C   RESTORE X, PHI(*,*) AND PSI(*)
  C
        PHASE1 = .FALSE.
        X = XOLD
        DO 310 I = 1,K
          TEMP1 = 1.0D0/BETA(I)
          IP1 = I+1
          DO 305 L = 1,NEQN
   305      PHI(L,I) = TEMP1*(PHI(L,I) - PHI(L,IP1))
   310    CONTINUE
        IF(K .LT. 2) GO TO 320
        DO 315 I = 2,K
   315    PSI(I-1) = PSI(I) - H
  C
  C   ON THIRD FAILURE, SET ORDER TO ONE.  THEREAFTER, USE OPTIMAL STEP
  C   SIZE
  C
   320  IFAIL = IFAIL + 1
        TEMP2 = 0.5D0
        IF(IFAIL - 3) 335,330,325
   325  IF(P5EPS .LT. 0.25D0*ERK) TEMP2 = SQRT(P5EPS/ERK)
   330  KNEW = 1
   335  H = TEMP2*H
        K = KNEW
        NS = 0
        IF(ABS(H) .GE. FOURU*ABS(X)) GO TO 340
        CRASH = .TRUE.
        H = SIGN(FOURU*ABS(X),H)
        EPS = EPS + EPS
        RETURN
   340  GO TO 100
  C       ***     END BLOCK 3     ***
  C
  C       ***     BEGIN BLOCK 4     ***
  C   THE STEP IS SUCCESSFUL.  CORRECT THE PREDICTED SOLUTION, EVALUATE
  C   THE DERIVATIVES USING THE CORRECTED SOLUTION AND UPDATE THE
  C   DIFFERENCES.  DETERMINE BEST ORDER AND STEP SIZE FOR NEXT STEP.
  C                   ***
   400  KOLD = K
        HOLD = H
  C
  C   CORRECT AND EVALUATE
  C
        TEMP1 = H*G(KP1)
        IF(NORND) GO TO 410
        DO 405 L = 1,NEQN
          TEMP3 = Y(L)
          RHO = TEMP1*(YP(L) - PHI(L,1)) - PHI(L,16)
          Y(L) = P(L) + RHO
          PHI(L,15) = (Y(L) - P(L)) - RHO
   405    P(L) = TEMP3
        GO TO 420
   410  DO 415 L = 1,NEQN
          TEMP3 = Y(L)
          Y(L) = P(L) + TEMP1*(YP(L) - PHI(L,1))
   415    P(L) = TEMP3
   420  CALL DF(X,Y,YP,RPAR,IPAR)
  C
  C   UPDATE DIFFERENCES FOR NEXT STEP
  C
        DO 425 L = 1,NEQN
          PHI(L,KP1) = YP(L) - PHI(L,1)
   425    PHI(L,KP2) = PHI(L,KP1) - PHI(L,KP2)
        DO 435 I = 1,K
          DO 430 L = 1,NEQN
   430      PHI(L,I) = PHI(L,I) + PHI(L,KP1)
   435    CONTINUE
  C
  C   ESTIMATE ERROR AT ORDER K+1 UNLESS:
  C     IN FIRST PHASE WHEN ALWAYS RAISE ORDER,
  C     ALREADY DECIDED TO LOWER ORDER,
  C     STEP SIZE NOT CONSTANT SO ESTIMATE UNRELIABLE
  C
        ERKP1 = 0.0D0
        IF(KNEW .EQ. KM1  .OR.  K .EQ. 12) PHASE1 = .FALSE.
        IF(PHASE1) GO TO 450
        IF(KNEW .EQ. KM1) GO TO 455
        IF(KP1 .GT. NS) GO TO 460
        DO 440 L = 1,NEQN
   440    ERKP1 = ERKP1 + (PHI(L,KP2)/WT(L))**2
        ERKP1 = ABSH*GSTR(KP1)*SQRT(ERKP1)
  C
  C   USING ESTIMATED ERROR AT ORDER K+1, DETERMINE APPROPRIATE ORDER
  C   FOR NEXT STEP
  C
        IF(K .GT. 1) GO TO 445
        IF(ERKP1 .GE. 0.5D0*ERK) GO TO 460
        GO TO 450
   445  IF(ERKM1 .LE. MIN(ERK,ERKP1)) GO TO 455
        IF(ERKP1 .GE. ERK  .OR.  K .EQ. 12) GO TO 460
  C
  C   HERE ERKP1 .LT. ERK .LT. MAX(ERKM1,ERKM2) ELSE ORDER WOULD HAVE
  C   BEEN LOWERED IN BLOCK 2.  THUS ORDER IS TO BE RAISED
  C
  C   RAISE ORDER
  C
   450  K = KP1
        ERK = ERKP1
        GO TO 460
  C
  C   LOWER ORDER
  C
   455  K = KM1
        ERK = ERKM1
  C
  C   WITH NEW ORDER DETERMINE APPROPRIATE STEP SIZE FOR NEXT STEP
  C
   460  HNEW = H + H
        IF(PHASE1) GO TO 465
        IF(P5EPS .GE. ERK*TWO(K+1)) GO TO 465
        HNEW = H
        IF(P5EPS .GE. ERK) GO TO 465
        TEMP2 = K+1
        R = (P5EPS/ERK)**(1.0D0/TEMP2)
        HNEW = ABSH*MAX(0.5D0,MIN(0.9D0,R))
        HNEW = SIGN(MAX(HNEW,FOURU*ABS(X)),H)
   465  H = HNEW
        RETURN
  C       ***     END BLOCK 4     ***
        END
  *DECK FDUMP
        SUBROUTINE FDUMP
  C***BEGIN PROLOGUE  FDUMP
  C***PURPOSE  Symbolic dump (should be locally written).
  C***LIBRARY   SLATEC (XERROR)
  C***CATEGORY  R3
  C***TYPE      ALL (FDUMP-A)
  C***KEYWORDS  ERROR, XERMSG
  C***AUTHOR  Jones, R. E., (SNLA)
  C***DESCRIPTION
  C
  C        ***Note*** Machine Dependent Routine
  C        FDUMP is intended to be replaced by a locally written
  C        version which produces a symbolic dump.  Failing this,
  C        it should be replaced by a version which prints the
  C        subprogram nesting list.  Note that this dump must be
  C        printed on each of up to five files, as indicated by the
  C        XGETUA routine.  See XSETUA and XGETUA for details.
  C
  C     Written by Ron Jones, with SLATEC Common Math Library Subcommittee
  C
  C***REFERENCES  (NONE)
  C***ROUTINES CALLED  (NONE)
  C***REVISION HISTORY  (YYMMDD)
  C   790801  DATE WRITTEN
  C   861211  REVISION DATE from Version 3.2
  C   891214  Prologue converted to Version 4.0 format.  (BAB)
  C***END PROLOGUE  FDUMP
  C***FIRST EXECUTABLE STATEMENT  FDUMP
        RETURN
        END
  *DECK I1MACH
        INTEGER FUNCTION I1MACH (I)
  C***BEGIN PROLOGUE  I1MACH
  C***PURPOSE  Return integer machine dependent constants.
  C***LIBRARY   SLATEC
  C***CATEGORY  R1
  C***TYPE      INTEGER (I1MACH-I)
  C***KEYWORDS  MACHINE CONSTANTS
  C***AUTHOR  Fox, P. A., (Bell Labs)
  C           Hall, A. D., (Bell Labs)
  C           Schryer, N. L., (Bell Labs)
  C***DESCRIPTION
  C
  C   I1MACH can be used to obtain machine-dependent parameters for the
  C   local machine environment.  It is a function subprogram with one
  C   (input) argument and can be referenced as follows:
  C
  C        K = I1MACH(I)
  C
  C   where I=1,...,16.  The (output) value of K above is determined by
  C   the (input) value of I.  The results for various values of I are
  C   discussed below.
  C
  C   I/O unit numbers:
  C     I1MACH( 1) = the standard input unit.
  C     I1MACH( 2) = the standard output unit.
  C     I1MACH( 3) = the standard punch unit.
  C     I1MACH( 4) = the standard error message unit.
  C
  C   Words:
  C     I1MACH( 5) = the number of bits per integer storage unit.
  C     I1MACH( 6) = the number of characters per integer storage unit.
  C
  C   Integers:
  C     assume integers are represented in the S-digit, base-A form
  C
  C                sign ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) )
  C
  C                where 0 .LE. X(I) .LT. A for I=0,...,S-1.
  C     I1MACH( 7) = A, the base.
  C     I1MACH( 8) = S, the number of base-A digits.
  C     I1MACH( 9) = A**S - 1, the largest magnitude.
  C
  C   Floating-Point Numbers:
  C     Assume floating-point numbers are represented in the T-digit,
  C     base-B form
  C                sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
  C
  C                where 0 .LE. X(I) .LT. B for I=1,...,T,
  C                0 .LT. X(1), and EMIN .LE. E .LE. EMAX.
  C     I1MACH(10) = B, the base.
  C
  C   Single-Precision:
  C     I1MACH(11) = T, the number of base-B digits.
  C     I1MACH(12) = EMIN, the smallest exponent E.
  C     I1MACH(13) = EMAX, the largest exponent E.
  C
  C   Double-Precision:
  C     I1MACH(14) = T, the number of base-B digits.
  C     I1MACH(15) = EMIN, the smallest exponent E.
  C     I1MACH(16) = EMAX, the largest exponent E.
  C
  C   To alter this function for a particular environment, the desired
  C   set of DATA statements should be activated by removing the C from
  C   column 1.  Also, the values of I1MACH(1) - I1MACH(4) should be
  C   checked for consistency with the local operating system.
  C
  C***REFERENCES  P. A. Fox, A. D. Hall and N. L. Schryer, Framework for
  C                 a portable library, ACM Transactions on Mathematical
  C                 Software 4, 2 (June 1978), pp. 177-188.
  C***ROUTINES CALLED  (NONE)
  C***REVISION HISTORY  (YYMMDD)
  C   750101  DATE WRITTEN
  C   891012  Added VAX G-floating constants.  (WRB)
  C   891012  REVISION DATE from Version 3.2
  C   891214  Prologue converted to Version 4.0 format.  (BAB)
  C   900618  Added DEC RISC constants.  (WRB)
  C   900723  Added IBM RS 6000 constants.  (WRB)
  C   901009  Correct I1MACH(7) for IBM Mainframes. Should be 2 not 16.
  C           (RWC)
  C   910710  Added HP 730 constants.  (SMR)
  C   911114  Added Convex IEEE constants.  (WRB)
  C   920121  Added SUN -r8 compiler option constants.  (WRB)
  C   920229  Added Touchstone Delta i860 constants.  (WRB)
  C   920501  Reformatted the REFERENCES section.  (WRB)
  C   920625  Added Convex -p8 and -pd8 compiler option constants.
  C           (BKS, WRB)
  C   930201  Added DEC Alpha and SGI constants.  (RWC and WRB)
  C   930618  Corrected I1MACH(5) for Convex -p8 and -pd8 compiler
  C           options.  (DWL, RWC and WRB).
  C   010817  Elevated IEEE to highest importance; see next set of
  C           comments below.  (DWL)
  C***END PROLOGUE  I1MACH
  C
  C Initial data here correspond to the IEEE standard.  If one of the
  C sets of initial data below is preferred, do the necessary commenting
  C and uncommenting. (DWL)
        INTEGER IMACH(16),OUTPUT
        DATA IMACH( 1) /          5 /
        DATA IMACH( 2) /          6 /
        DATA IMACH( 3) /          6 /
        DATA IMACH( 4) /          6 /
        DATA IMACH( 5) /         32 /
        DATA IMACH( 6) /          4 /
        DATA IMACH( 7) /          2 /
        DATA IMACH( 8) /         31 /
        DATA IMACH( 9) / 2147483647 /
        DATA IMACH(10) /          2 /
        DATA IMACH(11) /         24 /
        DATA IMACH(12) /       -126 /
        DATA IMACH(13) /        127 /
        DATA IMACH(14) /         53 /
        DATA IMACH(15) /      -1022 /
        DATA IMACH(16) /       1023 /
        SAVE IMACH
  cc      EQUIVALENCE (IMACH(4),OUTPUT)
  C
  C     MACHINE CONSTANTS FOR THE AMIGA
  C     ABSOFT COMPILER
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          5 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         32 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         31 /
  C     DATA IMACH( 9) / 2147483647 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -126 /
  C     DATA IMACH(13) /        127 /
  C     DATA IMACH(14) /         53 /
  C     DATA IMACH(15) /      -1022 /
  C     DATA IMACH(16) /       1023 /
  C
  C     MACHINE CONSTANTS FOR THE APOLLO
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          6 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         32 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         31 /
  C     DATA IMACH( 9) / 2147483647 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -125 /
  C     DATA IMACH(13) /        129 /
  C     DATA IMACH(14) /         53 /
  C     DATA IMACH(15) /      -1021 /
  C     DATA IMACH(16) /       1025 /
  C
  C     MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM
  C
  C     DATA IMACH( 1) /          7 /
  C     DATA IMACH( 2) /          2 /
  C     DATA IMACH( 3) /          2 /
  C     DATA IMACH( 4) /          2 /
  C     DATA IMACH( 5) /         36 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         33 /
  C     DATA IMACH( 9) / Z1FFFFFFFF /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -256 /
  C     DATA IMACH(13) /        255 /
  C     DATA IMACH(14) /         60 /
  C     DATA IMACH(15) /       -256 /
  C     DATA IMACH(16) /        255 /
  C
  C     MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          7 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         48 /
  C     DATA IMACH( 6) /          6 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         39 /
  C     DATA IMACH( 9) / O0007777777777777 /
  C     DATA IMACH(10) /          8 /
  C     DATA IMACH(11) /         13 /
  C     DATA IMACH(12) /        -50 /
  C     DATA IMACH(13) /         76 /
  C     DATA IMACH(14) /         26 /
  C     DATA IMACH(15) /        -50 /
  C     DATA IMACH(16) /         76 /
  C
  C     MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          7 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         48 /
  C     DATA IMACH( 6) /          6 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         39 /
  C     DATA IMACH( 9) / O0007777777777777 /
  C     DATA IMACH(10) /          8 /
  C     DATA IMACH(11) /         13 /
  C     DATA IMACH(12) /        -50 /
  C     DATA IMACH(13) /         76 /
  C     DATA IMACH(14) /         26 /
  C     DATA IMACH(15) /     -32754 /
  C     DATA IMACH(16) /      32780 /
  C
  C     MACHINE CONSTANTS FOR THE CDC 170/180 SERIES USING NOS/VE
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          7 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         64 /
  C     DATA IMACH( 6) /          8 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         63 /
  C     DATA IMACH( 9) / 9223372036854775807 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         47 /
  C     DATA IMACH(12) /      -4095 /
  C     DATA IMACH(13) /       4094 /
  C     DATA IMACH(14) /         94 /
  C     DATA IMACH(15) /      -4095 /
  C     DATA IMACH(16) /       4094 /
  C
  C     MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          7 /
  C     DATA IMACH( 4) /    6LOUTPUT/
  C     DATA IMACH( 5) /         60 /
  C     DATA IMACH( 6) /         10 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         48 /
  C     DATA IMACH( 9) / 00007777777777777777B /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         47 /
  C     DATA IMACH(12) /       -929 /
  C     DATA IMACH(13) /       1070 /
  C     DATA IMACH(14) /         94 /
  C     DATA IMACH(15) /       -929 /
  C     DATA IMACH(16) /       1069 /
  C
  C     MACHINE CONSTANTS FOR THE CELERITY C1260
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          6 /
  C     DATA IMACH( 4) /          0 /
  C     DATA IMACH( 5) /         32 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         31 /
  C     DATA IMACH( 9) / Z'7FFFFFFF' /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -126 /
  C     DATA IMACH(13) /        127 /
  C     DATA IMACH(14) /         53 /
  C     DATA IMACH(15) /      -1022 /
  C     DATA IMACH(16) /       1023 /
  C
  C     MACHINE CONSTANTS FOR THE CONVEX
  C     USING THE -fn COMPILER OPTION
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          7 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         32 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         31 /
  C     DATA IMACH( 9) / 2147483647 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -127 /
  C     DATA IMACH(13) /        127 /
  C     DATA IMACH(14) /         53 /
  C     DATA IMACH(15) /      -1023 /
  C     DATA IMACH(16) /       1023 /
  C
  C     MACHINE CONSTANTS FOR THE CONVEX
  C     USING THE -fi COMPILER OPTION
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          7 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         32 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         31 /
  C     DATA IMACH( 9) / 2147483647 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -125 /
  C     DATA IMACH(13) /        128 /
  C     DATA IMACH(14) /         53 /
  C     DATA IMACH(15) /      -1021 /
  C     DATA IMACH(16) /       1024 /
  C
  C     MACHINE CONSTANTS FOR THE CONVEX
  C     USING THE -p8 COMPILER OPTION
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          7 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         64 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         63 /
  C     DATA IMACH( 9) / 9223372036854775807 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         53 /
  C     DATA IMACH(12) /      -1023 /
  C     DATA IMACH(13) /       1023 /
  C     DATA IMACH(14) /        113 /
  C     DATA IMACH(15) /     -16383 /
  C     DATA IMACH(16) /      16383 /
  C
  C     MACHINE CONSTANTS FOR THE CONVEX
  C     USING THE -pd8 COMPILER OPTION
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          7 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         64 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         63 /
  C     DATA IMACH( 9) / 9223372036854775807 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         53 /
  C     DATA IMACH(12) /      -1023 /
  C     DATA IMACH(13) /       1023 /
  C     DATA IMACH(14) /         53 /
  C     DATA IMACH(15) /      -1023 /
  C     DATA IMACH(16) /       1023 /
  C
  C     MACHINE CONSTANTS FOR THE CRAY
  C     USING THE 46 BIT INTEGER COMPILER OPTION
  C
  C     DATA IMACH( 1) /        100 /
  C     DATA IMACH( 2) /        101 /
  C     DATA IMACH( 3) /        102 /
  C     DATA IMACH( 4) /        101 /
  C     DATA IMACH( 5) /         64 /
  C     DATA IMACH( 6) /          8 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         46 /
  C     DATA IMACH( 9) / 1777777777777777B /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         47 /
  C     DATA IMACH(12) /      -8189 /
  C     DATA IMACH(13) /       8190 /
  C     DATA IMACH(14) /         94 /
  C     DATA IMACH(15) /      -8099 /
  C     DATA IMACH(16) /       8190 /
  C
  C     MACHINE CONSTANTS FOR THE CRAY
  C     USING THE 64 BIT INTEGER COMPILER OPTION
  C
  C     DATA IMACH( 1) /        100 /
  C     DATA IMACH( 2) /        101 /
  C     DATA IMACH( 3) /        102 /
  C     DATA IMACH( 4) /        101 /
  C     DATA IMACH( 5) /         64 /
  C     DATA IMACH( 6) /          8 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         63 /
  C     DATA IMACH( 9) / 777777777777777777777B /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         47 /
  C     DATA IMACH(12) /      -8189 /
  C     DATA IMACH(13) /       8190 /
  C     DATA IMACH(14) /         94 /
  C     DATA IMACH(15) /      -8099 /
  C     DATA IMACH(16) /       8190 /
  C
  C     MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200
  C
  C     DATA IMACH( 1) /         11 /
  C     DATA IMACH( 2) /         12 /
  C     DATA IMACH( 3) /          8 /
  C     DATA IMACH( 4) /         10 /
  C     DATA IMACH( 5) /         16 /
  C     DATA IMACH( 6) /          2 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         15 /
  C     DATA IMACH( 9) /      32767 /
  C     DATA IMACH(10) /         16 /
  C     DATA IMACH(11) /          6 /
  C     DATA IMACH(12) /        -64 /
  C     DATA IMACH(13) /         63 /
  C     DATA IMACH(14) /         14 /
  C     DATA IMACH(15) /        -64 /
  C     DATA IMACH(16) /         63 /
  C
  C     MACHINE CONSTANTS FOR THE DEC ALPHA
  C     USING G_FLOAT
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          5 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         32 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         31 /
  C     DATA IMACH( 9) / 2147483647 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -127 /
  C     DATA IMACH(13) /        127 /
  C     DATA IMACH(14) /         53 /
  C     DATA IMACH(15) /      -1023 /
  C     DATA IMACH(16) /       1023 /
  C
  C     MACHINE CONSTANTS FOR THE DEC ALPHA
  C     USING IEEE_FLOAT
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          6 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         32 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         31 /
  C     DATA IMACH( 9) / 2147483647 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -125 /
  C     DATA IMACH(13) /        128 /
  C     DATA IMACH(14) /         53 /
  C     DATA IMACH(15) /      -1021 /
  C     DATA IMACH(16) /       1024 /
  C
  C     MACHINE CONSTANTS FOR THE DEC RISC
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          6 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         32 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         31 /
  C     DATA IMACH( 9) / 2147483647 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -125 /
  C     DATA IMACH(13) /        128 /
  C     DATA IMACH(14) /         53 /
  C     DATA IMACH(15) /      -1021 /
  C     DATA IMACH(16) /       1024 /
  C
  C     MACHINE CONSTANTS FOR THE DEC VAX
  C     USING D_FLOATING
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          5 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         32 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         31 /
  C     DATA IMACH( 9) / 2147483647 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -127 /
  C     DATA IMACH(13) /        127 /
  C     DATA IMACH(14) /         56 /
  C     DATA IMACH(15) /       -127 /
  C     DATA IMACH(16) /        127 /
  C
  C     MACHINE CONSTANTS FOR THE DEC VAX
  C     USING G_FLOATING
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          5 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         32 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         31 /
  C     DATA IMACH( 9) / 2147483647 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -127 /
  C     DATA IMACH(13) /        127 /
  C     DATA IMACH(14) /         53 /
  C     DATA IMACH(15) /      -1023 /
  C     DATA IMACH(16) /       1023 /
  C
  C     MACHINE CONSTANTS FOR THE ELXSI 6400
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          6 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         32 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         32 /
  C     DATA IMACH( 9) / 2147483647 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -126 /
  C     DATA IMACH(13) /        127 /
  C     DATA IMACH(14) /         53 /
  C     DATA IMACH(15) /      -1022 /
  C     DATA IMACH(16) /       1023 /
  C
  C     MACHINE CONSTANTS FOR THE HARRIS 220
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          0 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         24 /
  C     DATA IMACH( 6) /          3 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         23 /
  C     DATA IMACH( 9) /    8388607 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         23 /
  C     DATA IMACH(12) /       -127 /
  C     DATA IMACH(13) /        127 /
  C     DATA IMACH(14) /         38 /
  C     DATA IMACH(15) /       -127 /
  C     DATA IMACH(16) /        127 /
  C
  C     MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /         43 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         36 /
  C     DATA IMACH( 6) /          6 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         35 /
  C     DATA IMACH( 9) / O377777777777 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         27 /
  C     DATA IMACH(12) /       -127 /
  C     DATA IMACH(13) /        127 /
  C     DATA IMACH(14) /         63 /
  C     DATA IMACH(15) /       -127 /
  C     DATA IMACH(16) /        127 /
  C
  C     MACHINE CONSTANTS FOR THE HP 730
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          6 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         32 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         31 /
  C     DATA IMACH( 9) / 2147483647 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -125 /
  C     DATA IMACH(13) /        128 /
  C     DATA IMACH(14) /         53 /
  C     DATA IMACH(15) /      -1021 /
  C     DATA IMACH(16) /       1024 /
  C
  C     MACHINE CONSTANTS FOR THE HP 2100
  C     3 WORD DOUBLE PRECISION OPTION WITH FTN4
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          4 /
  C     DATA IMACH( 4) /          1 /
  C     DATA IMACH( 5) /         16 /
  C     DATA IMACH( 6) /          2 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         15 /
  C     DATA IMACH( 9) /      32767 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         23 /
  C     DATA IMACH(12) /       -128 /
  C     DATA IMACH(13) /        127 /
  C     DATA IMACH(14) /         39 /
  C     DATA IMACH(15) /       -128 /
  C     DATA IMACH(16) /        127 /
  C
  C     MACHINE CONSTANTS FOR THE HP 2100
  C     4 WORD DOUBLE PRECISION OPTION WITH FTN4
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          4 /
  C     DATA IMACH( 4) /          1 /
  C     DATA IMACH( 5) /         16 /
  C     DATA IMACH( 6) /          2 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         15 /
  C     DATA IMACH( 9) /      32767 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         23 /
  C     DATA IMACH(12) /       -128 /
  C     DATA IMACH(13) /        127 /
  C     DATA IMACH(14) /         55 /
  C     DATA IMACH(15) /       -128 /
  C     DATA IMACH(16) /        127 /
  C
  C     MACHINE CONSTANTS FOR THE HP 9000
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          6 /
  C     DATA IMACH( 4) /          7 /
  C     DATA IMACH( 5) /         32 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         32 /
  C     DATA IMACH( 9) / 2147483647 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -126 /
  C     DATA IMACH(13) /        127 /
  C     DATA IMACH(14) /         53 /
  C     DATA IMACH(15) /      -1015 /
  C     DATA IMACH(16) /       1017 /
  C
  C     MACHINE CONSTANTS FOR THE IBM 360/370 SERIES,
  C     THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND
  C     THE PERKIN ELMER (INTERDATA) 7/32.
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          7 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         32 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         31 /
  C     DATA IMACH( 9) /  Z7FFFFFFF /
  C     DATA IMACH(10) /         16 /
  C     DATA IMACH(11) /          6 /
  C     DATA IMACH(12) /        -64 /
  C     DATA IMACH(13) /         63 /
  C     DATA IMACH(14) /         14 /
  C     DATA IMACH(15) /        -64 /
  C     DATA IMACH(16) /         63 /
  C
  C     MACHINE CONSTANTS FOR THE IBM PC
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          0 /
  C     DATA IMACH( 4) /          0 /
  C     DATA IMACH( 5) /         32 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         31 /
  C     DATA IMACH( 9) / 2147483647 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -125 /
  C     DATA IMACH(13) /        127 /
  C     DATA IMACH(14) /         53 /
  C     DATA IMACH(15) /      -1021 /
  C     DATA IMACH(16) /       1023 /
  C
  C     MACHINE CONSTANTS FOR THE IBM RS 6000
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          6 /
  C     DATA IMACH( 4) /          0 /
  C     DATA IMACH( 5) /         32 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         31 /
  C     DATA IMACH( 9) / 2147483647 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -125 /
  C     DATA IMACH(13) /        128 /
  C     DATA IMACH(14) /         53 /
  C     DATA IMACH(15) /      -1021 /
  C     DATA IMACH(16) /       1024 /
  C
  C     MACHINE CONSTANTS FOR THE INTEL i860
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          6 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         32 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         31 /
  C     DATA IMACH( 9) / 2147483647 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -125 /
  C     DATA IMACH(13) /        128 /
  C     DATA IMACH(14) /         53 /
  C     DATA IMACH(15) /      -1021 /
  C     DATA IMACH(16) /       1024 /
  C
  C     MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR)
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          5 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         36 /
  C     DATA IMACH( 6) /          5 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         35 /
  C     DATA IMACH( 9) / "377777777777 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         27 /
  C     DATA IMACH(12) /       -128 /
  C     DATA IMACH(13) /        127 /
  C     DATA IMACH(14) /         54 /
  C     DATA IMACH(15) /       -101 /
  C     DATA IMACH(16) /        127 /
  C
  C     MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR)
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          5 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         36 /
  C     DATA IMACH( 6) /          5 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         35 /
  C     DATA IMACH( 9) / "377777777777 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         27 /
  C     DATA IMACH(12) /       -128 /
  C     DATA IMACH(13) /        127 /
  C     DATA IMACH(14) /         62 /
  C     DATA IMACH(15) /       -128 /
  C     DATA IMACH(16) /        127 /
  C
  C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
  C     32-BIT INTEGER ARITHMETIC.
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          5 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         32 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         31 /
  C     DATA IMACH( 9) / 2147483647 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -127 /
  C     DATA IMACH(13) /        127 /
  C     DATA IMACH(14) /         56 /
  C     DATA IMACH(15) /       -127 /
  C     DATA IMACH(16) /        127 /
  C
  C     MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING
  C     16-BIT INTEGER ARITHMETIC.
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          5 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         16 /
  C     DATA IMACH( 6) /          2 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         15 /
  C     DATA IMACH( 9) /      32767 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -127 /
  C     DATA IMACH(13) /        127 /
  C     DATA IMACH(14) /         56 /
  C     DATA IMACH(15) /       -127 /
  C     DATA IMACH(16) /        127 /
  C
  C     MACHINE CONSTANTS FOR THE SILICON GRAPHICS
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          6 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         32 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         31 /
  C     DATA IMACH( 9) / 2147483647 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -125 /
  C     DATA IMACH(13) /        128 /
  C     DATA IMACH(14) /         53 /
  C     DATA IMACH(15) /      -1021 /
  C     DATA IMACH(16) /       1024 /
  C
  C     MACHINE CONSTANTS FOR THE SUN
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          6 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         32 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         31 /
  C     DATA IMACH( 9) / 2147483647 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -125 /
  C     DATA IMACH(13) /        128 /
  C     DATA IMACH(14) /         53 /
  C     DATA IMACH(15) /      -1021 /
  C     DATA IMACH(16) /       1024 /
  C
  C     MACHINE CONSTANTS FOR THE SUN
  C     USING THE -r8 COMPILER OPTION
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          6 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         32 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         31 /
  C     DATA IMACH( 9) / 2147483647 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         53 /
  C     DATA IMACH(12) /      -1021 /
  C     DATA IMACH(13) /       1024 /
  C     DATA IMACH(14) /        113 /
  C     DATA IMACH(15) /     -16381 /
  C     DATA IMACH(16) /      16384 /
  C
  C     MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES FTN COMPILER
  C
  C     DATA IMACH( 1) /          5 /
  C     DATA IMACH( 2) /          6 /
  C     DATA IMACH( 3) /          1 /
  C     DATA IMACH( 4) /          6 /
  C     DATA IMACH( 5) /         36 /
  C     DATA IMACH( 6) /          4 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         35 /
  C     DATA IMACH( 9) / O377777777777 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         27 /
  C     DATA IMACH(12) /       -128 /
  C     DATA IMACH(13) /        127 /
  C     DATA IMACH(14) /         60 /
  C     DATA IMACH(15) /      -1024 /
  C     DATA IMACH(16) /       1023 /
  C
  C     MACHINE CONSTANTS FOR THE Z80 MICROPROCESSOR
  C
  C     DATA IMACH( 1) /          1 /
  C     DATA IMACH( 2) /          1 /
  C     DATA IMACH( 3) /          0 /
  C     DATA IMACH( 4) /          1 /
  C     DATA IMACH( 5) /         16 /
  C     DATA IMACH( 6) /          2 /
  C     DATA IMACH( 7) /          2 /
  C     DATA IMACH( 8) /         15 /
  C     DATA IMACH( 9) /      32767 /
  C     DATA IMACH(10) /          2 /
  C     DATA IMACH(11) /         24 /
  C     DATA IMACH(12) /       -127 /
  C     DATA IMACH(13) /        127 /
  C     DATA IMACH(14) /         56 /
  C     DATA IMACH(15) /       -127 /
  C     DATA IMACH(16) /        127 /
  C
  C***FIRST EXECUTABLE STATEMENT  I1MACH
        IF (I .LT. 1  .OR.  I .GT. 16) GO TO 10
  C
        I1MACH = IMACH(I)
        RETURN
  C
     10 CONTINUE
        WRITE (UNIT = OUTPUT, FMT = 9000)
   9000 FORMAT ('1ERROR    1 IN I1MACH - I OUT OF BOUNDS')
  C
  C     CALL FDUMP
  C
        STOP
        END
  *DECK DHSTRT
        SUBROUTINE DHSTRT (DF, NEQ, A, B, Y, YPRIME, ETOL, MORDER, SMALL,
       +   BIG, SPY, PV, YP, SF, RPAR, IPAR, H)
  C***BEGIN PROLOGUE  DHSTRT
  C***SUBSIDIARY
  C***PURPOSE  Subsidiary to DDEABM, DDEBDF and DDERKF
  C***LIBRARY   SLATEC
  C***TYPE      DOUBLE PRECISION (HSTART-S, DHSTRT-D)
  C***AUTHOR  Watts, H. A., (SNLA)
  C***DESCRIPTION
  C
  C   DHSTRT computes a starting step size to be used in solving initial
  C   value problems in ordinary differential equations.
  C
  C **********************************************************************
  C  ABSTRACT
  C
  C     Subroutine DHSTRT computes a starting step size to be used by an
  C     initial value method in solving ordinary differential equations.
  C     It is based on an estimate of the local Lipschitz constant for the
  C     differential equation   (lower bound on a norm of the Jacobian) ,
  C     a bound on the differential equation  (first derivative) , and
  C     a bound on the partial derivative of the equation with respect to
  C     the independent variable.
  C     (all approximated near the initial point A)
  C
  C     Subroutine DHSTRT uses a function subprogram DHVNRM for computing
  C     a vector norm. The maximum norm is presently utilized though it
  C     can easily be replaced by any other vector norm. It is presumed
  C     that any replacement norm routine would be carefully coded to
  C     prevent unnecessary underflows or overflows from occurring, and
  C     also, would not alter the vector or number of components.
  C
  C **********************************************************************
  C  On input you must provide the following
  C
  C      DF -- This is a subroutine of the form
  C                               DF(X,U,UPRIME,RPAR,IPAR)
  C             which defines the system of first order differential
  C             equations to be solved. For the given values of X and the
  C             vector  U(*)=(U(1),U(2),...,U(NEQ)) , the subroutine must
  C             evaluate the NEQ components of the system of differential
  C             equations  DU/DX=DF(X,U)  and store the derivatives in the
  C             array UPRIME(*), that is,  UPRIME(I) = * DU(I)/DX *  for
  C             equations I=1,...,NEQ.
  C
  C             Subroutine DF must not alter X or U(*). You must declare
  C             the name DF in an external statement in your program that
  C             calls DHSTRT. You must dimension U and UPRIME in DF.
  C
  C             RPAR and IPAR are DOUBLE PRECISION and INTEGER parameter
  C             arrays which you can use for communication between your
  C             program and subroutine DF. They are not used or altered by
  C             DHSTRT. If you do not need RPAR or IPAR, ignore these
  C             parameters by treating them as dummy arguments. If you do
  C             choose to use them, dimension them in your program and in
  C             DF as arrays of appropriate length.
  C
  C      NEQ -- This is the number of (first order) differential equations
  C             to be integrated.
  C
  C      A -- This is the initial point of integration.
  C
  C      B -- This is a value of the independent variable used to define
  C             the direction of integration. A reasonable choice is to
  C             set  B  to the first point at which a solution is desired.
  C             You can also use  B, if necessary, to restrict the length
  C             of the first integration step because the algorithm will
  C             not compute a starting step length which is bigger than
  C             ABS(B-A), unless  B  has been chosen too close to  A.
  C             (it is presumed that DHSTRT has been called with  B
  C             different from  A  on the machine being used. Also see the
  C             discussion about the parameter  SMALL.)
  C
  C      Y(*) -- This is the vector of initial values of the NEQ solution
  C             components at the initial point  A.
  C
  C      YPRIME(*) -- This is the vector of derivatives of the NEQ
  C             solution components at the initial point  A.
  C             (defined by the differential equations in subroutine DF)
  C
  C      ETOL -- This is the vector of error tolerances corresponding to
  C             the NEQ solution components. It is assumed that all
  C             elements are positive. Following the first integration
  C             step, the tolerances are expected to be used by the
  C             integrator in an error test which roughly requires that
  C                        ABS(LOCAL ERROR)  .LE.  ETOL
  C             for each vector component.
  C
  C      MORDER -- This is the order of the formula which will be used by
  C             the initial value method for taking the first integration
  C             step.
  C
  C      SMALL -- This is a small positive machine dependent constant
  C             which is used for protecting against computations with
  C             numbers which are too small relative to the precision of
  C             floating point arithmetic.  SMALL  should be set to
  C             (approximately) the smallest positive DOUBLE PRECISION
  C             number such that  (1.+SMALL) .GT. 1.  on the machine being
  C             used. The quantity  SMALL**(3/8)  is used in computing
  C             increments of variables for approximating derivatives by
  C             differences.  Also the algorithm will not compute a
  C             starting step length which is smaller than
  C             100*SMALL*ABS(A).
  C
  C      BIG -- This is a large positive machine dependent constant which
  C             is used for preventing machine overflows. A reasonable
  C             choice is to set big to (approximately) the square root of
  C             the largest DOUBLE PRECISION number which can be held in
  C             the machine.
  C
  C      SPY(*),PV(*),YP(*),SF(*) -- These are DOUBLE PRECISION work
  C             arrays of length NEQ which provide the routine with needed
  C             storage space.
  C
  C      RPAR,IPAR -- These are parameter arrays, of DOUBLE PRECISION and
  C             INTEGER type, respectively, which can be used for
  C             communication between your program and the DF subroutine.
  C             They are not used or altered by DHSTRT.
  C
  C **********************************************************************
  C  On Output  (after the return from DHSTRT),
  C
  C      H -- is an appropriate starting step size to be attempted by the
  C             differential equation method.
  C
  C           All parameters in the call list remain unchanged except for
  C           the working arrays SPY(*),PV(*),YP(*), and SF(*).
  C
  C **********************************************************************
  C
  C***SEE ALSO  DDEABM, DDEBDF, DDERKF
  C***ROUTINES CALLED  DHVNRM
  C***REVISION HISTORY  (YYMMDD)
  C   820301  DATE WRITTEN
  C   890531  Changed all specific intrinsics to generic.  (WRB)
  C   890831  Modified array declarations.  (WRB)
  C   890911  Removed unnecessary intrinsics.  (WRB)
  C   891024  Changed references from DVNORM to DHVNRM.  (WRB)
  C   891214  Prologue converted to Version 4.0 format.  (BAB)
  C   900328  Added TYPE section.  (WRB)
  C   910722  Updated AUTHOR section.  (ALS)
  C***END PROLOGUE  DHSTRT
  C
        INTEGER IPAR, J, K, LK, MORDER, NEQ
        DOUBLE PRECISION A, ABSDX, B, BIG, DA, DELF, DELY,
       1      DFDUB, DFDXB, DHVNRM,
       2      DX, DY, ETOL, FBND, H, PV, RELPER, RPAR, SF, SMALL, SPY,
       3      SRYDPB, TOLEXP, TOLMIN, TOLP, TOLSUM, Y, YDPB, YP, YPRIME
        DIMENSION Y(*),YPRIME(*),ETOL(*),SPY(*),PV(*),YP(*),
       1          SF(*),RPAR(*),IPAR(*)
        EXTERNAL DF
  C
  C     ..................................................................
  C
  C     BEGIN BLOCK PERMITTING ...EXITS TO 160
  C***FIRST EXECUTABLE STATEMENT  DHSTRT
           DX = B - A
           ABSDX = ABS(DX)
           RELPER = SMALL**0.375D0
  C
  C        ...............................................................
  C
  C             COMPUTE AN APPROXIMATE BOUND (DFDXB) ON THE PARTIAL
  C             DERIVATIVE OF THE EQUATION WITH RESPECT TO THE
  C             INDEPENDENT VARIABLE. PROTECT AGAINST AN OVERFLOW.
  C             ALSO COMPUTE A BOUND (FBND) ON THE FIRST DERIVATIVE
  C             LOCALLY.
  C
           DA = SIGN(MAX(MIN(RELPER*ABS(A),ABSDX),
       1                    100.0D0*SMALL*ABS(A)),DX)
           IF (DA .EQ. 0.0D0) DA = RELPER*DX
           CALL DF(A+DA,Y,SF,RPAR,IPAR)
           DO 10 J = 1, NEQ
              YP(J) = SF(J) - YPRIME(J)
     10    CONTINUE
           DELF = DHVNRM(YP,NEQ)
           DFDXB = BIG
           IF (DELF .LT. BIG*ABS(DA)) DFDXB = DELF/ABS(DA)
           FBND = DHVNRM(SF,NEQ)
  C
  C        ...............................................................
  C
  C             COMPUTE AN ESTIMATE (DFDUB) OF THE LOCAL LIPSCHITZ
  C             CONSTANT FOR THE SYSTEM OF DIFFERENTIAL EQUATIONS. THIS
  C             ALSO REPRESENTS AN ESTIMATE OF THE NORM OF THE JACOBIAN
  C             LOCALLY.  THREE ITERATIONS (TWO WHEN NEQ=1) ARE USED TO
  C             ESTIMATE THE LIPSCHITZ CONSTANT BY NUMERICAL DIFFERENCES.
  C             THE FIRST PERTURBATION VECTOR IS BASED ON THE INITIAL
  C             DERIVATIVES AND DIRECTION OF INTEGRATION. THE SECOND
  C             PERTURBATION VECTOR IS FORMED USING ANOTHER EVALUATION OF
  C             THE DIFFERENTIAL EQUATION.  THE THIRD PERTURBATION VECTOR
  C             IS FORMED USING PERTURBATIONS BASED ONLY ON THE INITIAL
  C             VALUES. COMPONENTS THAT ARE ZERO ARE ALWAYS CHANGED TO
  C             NON-ZERO VALUES (EXCEPT ON THE FIRST ITERATION). WHEN
  C             INFORMATION IS AVAILABLE, CARE IS TAKEN TO ENSURE THAT
  C             COMPONENTS OF THE PERTURBATION VECTOR HAVE SIGNS WHICH ARE
  C             CONSISTENT WITH THE SLOPES OF LOCAL SOLUTION CURVES.
  C             ALSO CHOOSE THE LARGEST BOUND (FBND) FOR THE FIRST
  C             DERIVATIVE.
  C
  C                               PERTURBATION VECTOR SIZE IS HELD
  C                               CONSTANT FOR ALL ITERATIONS. COMPUTE
  C                               THIS CHANGE FROM THE
  C                                       SIZE OF THE VECTOR OF INITIAL
  C                                       VALUES.
           DELY = RELPER*DHVNRM(Y,NEQ)
           IF (DELY .EQ. 0.0D0) DELY = RELPER
           DELY = SIGN(DELY,DX)
           DELF = DHVNRM(YPRIME,NEQ)
           FBND = MAX(FBND,DELF)
           IF (DELF .EQ. 0.0D0) GO TO 30
  C           USE INITIAL DERIVATIVES FOR FIRST PERTURBATION
              DO 20 J = 1, NEQ
                 SPY(J) = YPRIME(J)
                 YP(J) = YPRIME(J)
     20       CONTINUE
           GO TO 50
     30    CONTINUE
  C           CANNOT HAVE A NULL PERTURBATION VECTOR
              DO 40 J = 1, NEQ
                 SPY(J) = 0.0D0
                 YP(J) = 1.0D0
     40       CONTINUE
              DELF = DHVNRM(YP,NEQ)
     50    CONTINUE
  C
           DFDUB = 0.0D0
           LK = MIN(NEQ+1,3)
           DO 140 K = 1, LK
  C           DEFINE PERTURBED VECTOR OF INITIAL VALUES
              DO 60 J = 1, NEQ
                 PV(J) = Y(J) + DELY*(YP(J)/DELF)
     60       CONTINUE
              IF (K .EQ. 2) GO TO 80
  C              EVALUATE DERIVATIVES ASSOCIATED WITH PERTURBED
  C              VECTOR  AND  COMPUTE CORRESPONDING DIFFERENCES
                 CALL DF(A,PV,YP,RPAR,IPAR)
                 DO 70 J = 1, NEQ
                    PV(J) = YP(J) - YPRIME(J)
     70          CONTINUE
              GO TO 100
     80       CONTINUE
  C              USE A SHIFTED VALUE OF THE INDEPENDENT VARIABLE
  C                                    IN COMPUTING ONE ESTIMATE
                 CALL DF(A+DA,PV,YP,RPAR,IPAR)
                 DO 90 J = 1, NEQ
                    PV(J) = YP(J) - SF(J)
     90          CONTINUE
    100       CONTINUE
  C           CHOOSE LARGEST BOUNDS ON THE FIRST DERIVATIVE
  C                          AND A LOCAL LIPSCHITZ CONSTANT
              FBND = MAX(FBND,DHVNRM(YP,NEQ))
              DELF = DHVNRM(PV,NEQ)
  C        ...EXIT
              IF (DELF .GE. BIG*ABS(DELY)) GO TO 150
              DFDUB = MAX(DFDUB,DELF/ABS(DELY))
  C     ......EXIT
              IF (K .EQ. LK) GO TO 160
  C           CHOOSE NEXT PERTURBATION VECTOR
              IF (DELF .EQ. 0.0D0) DELF = 1.0D0
              DO 130 J = 1, NEQ
                 IF (K .EQ. 2) GO TO 110
                    DY = ABS(PV(J))
                    IF (DY .EQ. 0.0D0) DY = DELF
                 GO TO 120
    110          CONTINUE
                    DY = Y(J)
                    IF (DY .EQ. 0.0D0) DY = DELY/RELPER
    120          CONTINUE
                 IF (SPY(J) .EQ. 0.0D0) SPY(J) = YP(J)
                 IF (SPY(J) .NE. 0.0D0) DY = SIGN(DY,SPY(J))
                 YP(J) = DY
    130       CONTINUE
              DELF = DHVNRM(YP,NEQ)
    140    CONTINUE
    150    CONTINUE
  C
  C        PROTECT AGAINST AN OVERFLOW
           DFDUB = BIG
    160 CONTINUE
  C
  C     ..................................................................
  C
  C          COMPUTE A BOUND (YDPB) ON THE NORM OF THE SECOND DERIVATIVE
  C
        YDPB = DFDXB + DFDUB*FBND
  C
  C     ..................................................................
  C
  C          DEFINE THE TOLERANCE PARAMETER UPON WHICH THE STARTING STEP
  C          SIZE IS TO BE BASED.  A VALUE IN THE MIDDLE OF THE ERROR
  C          TOLERANCE RANGE IS SELECTED.
  C
        TOLMIN = BIG
        TOLSUM = 0.0D0
        DO 170 K = 1, NEQ
           TOLEXP = LOG10(ETOL(K))
           TOLMIN = MIN(TOLMIN,TOLEXP)
           TOLSUM = TOLSUM + TOLEXP
    170 CONTINUE
        TOLP = 10.0D0**(0.5D0*(TOLSUM/NEQ + TOLMIN)/(MORDER+1))
  C
  C     ..................................................................
  C
  C          COMPUTE A STARTING STEP SIZE BASED ON THE ABOVE FIRST AND
  C          SECOND DERIVATIVE INFORMATION
  C
  C                            RESTRICT THE STEP LENGTH TO BE NOT BIGGER
  C                            THAN ABS(B-A).   (UNLESS  B  IS TOO CLOSE
  C                            TO  A)
        H = ABSDX
  C
        IF (YDPB .NE. 0.0D0 .OR. FBND .NE. 0.0D0) GO TO 180
  C
  C        BOTH FIRST DERIVATIVE TERM (FBND) AND SECOND
  C                     DERIVATIVE TERM (YDPB) ARE ZERO
           IF (TOLP .LT. 1.0D0) H = ABSDX*TOLP
        GO TO 200
    180 CONTINUE
  C
        IF (YDPB .NE. 0.0D0) GO TO 190
  C
  C        ONLY SECOND DERIVATIVE TERM (YDPB) IS ZERO
           IF (TOLP .LT. FBND*ABSDX) H = TOLP/FBND
        GO TO 200
    190 CONTINUE
  C
  C        SECOND DERIVATIVE TERM (YDPB) IS NON-ZERO
           SRYDPB = SQRT(0.5D0*YDPB)
           IF (TOLP .LT. SRYDPB*ABSDX) H = TOLP/SRYDPB
    200 CONTINUE
  C
  C     FURTHER RESTRICT THE STEP LENGTH TO BE NOT
  C                               BIGGER THAN  1/DFDUB
        IF (H*DFDUB .GT. 1.0D0) H = 1.0D0/DFDUB
  C
  C     FINALLY, RESTRICT THE STEP LENGTH TO BE NOT
  C     SMALLER THAN  100*SMALL*ABS(A).  HOWEVER, IF
  C     A=0. AND THE COMPUTED H UNDERFLOWED TO ZERO,
  C     THE ALGORITHM RETURNS  SMALL*ABS(B)  FOR THE
  C                                     STEP LENGTH.
        H = MAX(H,100.0D0*SMALL*ABS(A))
        IF (H .EQ. 0.0D0) H = SMALL*ABS(B)
  C
  C     NOW SET DIRECTION OF INTEGRATION
        H = SIGN(H,DX)
  C
        RETURN
        END
  *DECK DHVNRM
        DOUBLE PRECISION FUNCTION DHVNRM (V, NCOMP)
  C***BEGIN PROLOGUE  DHVNRM
  C***SUBSIDIARY
  C***PURPOSE  Subsidiary to DDEABM, DDEBDF and DDERKF
  C***LIBRARY   SLATEC
  C***TYPE      DOUBLE PRECISION (HVNRM-S, DHVNRM-D)
  C***AUTHOR  Watts, H. A., (SNLA)
  C***DESCRIPTION
  C
  C     Compute the maximum norm of the vector V(*) of length NCOMP and
  C     return the result as DHVNRM
  C
  C***SEE ALSO  DDEABM, DDEBDF, DDERKF
  C***ROUTINES CALLED  (NONE)
  C***REVISION HISTORY  (YYMMDD)
  C   820301  DATE WRITTEN
  C   890531  Changed all specific intrinsics to generic.  (WRB)
  C   890831  Modified array declarations.  (WRB)
  C   891024  Changed references from DVNORM to DHVNRM.  (WRB)
  C   891024  Changed routine name from DVNORM to DHVNRM.  (WRB)
  C   891214  Prologue converted to Version 4.0 format.  (BAB)
  C   900328  Added TYPE section.  (WRB)
  C   910722  Updated AUTHOR section.  (ALS)
  C***END PROLOGUE  DHVNRM
  C
        INTEGER K, NCOMP
        DOUBLE PRECISION V
        DIMENSION V(*)
  C***FIRST EXECUTABLE STATEMENT  DHVNRM
        DHVNRM = 0.0D0
        DO 10 K = 1, NCOMP
           DHVNRM = MAX(DHVNRM,ABS(V(K)))
     10 CONTINUE
        RETURN
        END
  *DECK J4SAVE
        FUNCTION J4SAVE (IWHICH, IVALUE, ISET)
  C***BEGIN PROLOGUE  J4SAVE
  C***SUBSIDIARY
  C***PURPOSE  Save or recall global variables needed by error
  C            handling routines.
  C***LIBRARY   SLATEC (XERROR)
  C***TYPE      INTEGER (J4SAVE-I)
  C***KEYWORDS  ERROR MESSAGES, ERROR NUMBER, RECALL, SAVE, XERROR
  C***AUTHOR  Jones, R. E., (SNLA)
  C***DESCRIPTION
  C
  C     Abstract
  C        J4SAVE saves and recalls several global variables needed
  C        by the library error handling routines.
  C
  C     Description of Parameters
  C      --Input--
  C        IWHICH - Index of item desired.
  C                = 1 Refers to current error number.
  C                = 2 Refers to current error control flag.
  C                = 3 Refers to current unit number to which error
  C                    messages are to be sent.  (0 means use standard.)
  C                = 4 Refers to the maximum number of times any
  C                     message is to be printed (as set by XERMAX).
  C                = 5 Refers to the total number of units to which
  C                     each error message is to be written.
  C                = 6 Refers to the 2nd unit for error messages
  C                = 7 Refers to the 3rd unit for error messages
  C                = 8 Refers to the 4th unit for error messages
  C                = 9 Refers to the 5th unit for error messages
  C        IVALUE - The value to be set for the IWHICH-th parameter,
  C                 if ISET is .TRUE. .
  C        ISET   - If ISET=.TRUE., the IWHICH-th parameter will BE
  C                 given the value, IVALUE.  If ISET=.FALSE., the
  C                 IWHICH-th parameter will be unchanged, and IVALUE
  C                 is a dummy parameter.
  C      --Output--
  C        The (old) value of the IWHICH-th parameter will be returned
  C        in the function value, J4SAVE.
  C
  C***SEE ALSO  XERMSG
  C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
  C                 Error-handling Package, SAND82-0800, Sandia
  C                 Laboratories, 1982.
  C***ROUTINES CALLED  (NONE)
  C***REVISION HISTORY  (YYMMDD)
  C   790801  DATE WRITTEN
  C   891214  Prologue converted to Version 4.0 format.  (BAB)
  C   900205  Minor modifications to prologue.  (WRB)
  C   900402  Added TYPE section.  (WRB)
  C   910411  Added KEYWORDS section.  (WRB)
  C   920501  Reformatted the REFERENCES section.  (WRB)
  C***END PROLOGUE  J4SAVE
        LOGICAL ISET
        INTEGER IPARAM(9)
        SAVE IPARAM
        DATA IPARAM(1),IPARAM(2),IPARAM(3),IPARAM(4)/0,2,0,10/
        DATA IPARAM(5)/1/
        DATA IPARAM(6),IPARAM(7),IPARAM(8),IPARAM(9)/0,0,0,0/
  C***FIRST EXECUTABLE STATEMENT  J4SAVE
        J4SAVE = IPARAM(IWHICH)
        IF (ISET) IPARAM(IWHICH) = IVALUE
        RETURN
        END
  *DECK XERCNT
        SUBROUTINE XERCNT (LIBRAR, SUBROU, MESSG, NERR, LEVEL, KONTRL)
  C***BEGIN PROLOGUE  XERCNT
  C***SUBSIDIARY
  C***PURPOSE  Allow user control over handling of errors.
  C***LIBRARY   SLATEC (XERROR)
  C***CATEGORY  R3C
  C***TYPE      ALL (XERCNT-A)
  C***KEYWORDS  ERROR, XERROR
  C***AUTHOR  Jones, R. E., (SNLA)
  C***DESCRIPTION
  C
  C     Abstract
  C        Allows user control over handling of individual errors.
  C        Just after each message is recorded, but before it is
  C        processed any further (i.e., before it is printed or
  C        a decision to abort is made), a call is made to XERCNT.
  C        If the user has provided his own version of XERCNT, he
  C        can then override the value of KONTROL used in processing
  C        this message by redefining its value.
  C        KONTRL may be set to any value from -2 to 2.
  C        The meanings for KONTRL are the same as in XSETF, except
  C        that the value of KONTRL changes only for this message.
  C        If KONTRL is set to a value outside the range from -2 to 2,
  C        it will be moved back into that range.
  C
  C     Description of Parameters
  C
  C      --Input--
  C        LIBRAR - the library that the routine is in.
  C        SUBROU - the subroutine that XERMSG is being called from
  C        MESSG  - the first 20 characters of the error message.
  C        NERR   - same as in the call to XERMSG.
  C        LEVEL  - same as in the call to XERMSG.
  C        KONTRL - the current value of the control flag as set
  C                 by a call to XSETF.
  C
  C      --Output--
  C        KONTRL - the new value of KONTRL.  If KONTRL is not
  C                 defined, it will remain at its original value.
  C                 This changed value of control affects only
  C                 the current occurrence of the current message.
  C
  C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
  C                 Error-handling Package, SAND82-0800, Sandia
  C                 Laboratories, 1982.
  C***ROUTINES CALLED  (NONE)
  C***REVISION HISTORY  (YYMMDD)
  C   790801  DATE WRITTEN
  C   861211  REVISION DATE from Version 3.2
  C   891214  Prologue converted to Version 4.0 format.  (BAB)
  C   900206  Routine changed from user-callable to subsidiary.  (WRB)
  C   900510  Changed calling sequence to include LIBRARY and SUBROUTINE
  C           names, changed routine name from XERCTL to XERCNT.  (RWC)
  C   920501  Reformatted the REFERENCES section.  (WRB)
  C***END PROLOGUE  XERCNT
        CHARACTER*(*) LIBRAR, SUBROU, MESSG
  C***FIRST EXECUTABLE STATEMENT  XERCNT
        RETURN
        END
  *DECK XERHLT
        SUBROUTINE XERHLT (MESSG)
  C***BEGIN PROLOGUE  XERHLT
  C***SUBSIDIARY
  C***PURPOSE  Abort program execution and print error message.
  C***LIBRARY   SLATEC (XERROR)
  C***CATEGORY  R3C
  C***TYPE      ALL (XERHLT-A)
  C***KEYWORDS  ABORT PROGRAM EXECUTION, ERROR, XERROR
  C***AUTHOR  Jones, R. E., (SNLA)
  C***DESCRIPTION
  C
  C     Abstract
  C        ***Note*** machine dependent routine
  C        XERHLT aborts the execution of the program.
  C        The error message causing the abort is given in the calling
  C        sequence, in case one needs it for printing on a dayfile,
  C        for example.
  C
  C     Description of Parameters
  C        MESSG is as in XERMSG.
  C
  C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
  C                 Error-handling Package, SAND82-0800, Sandia
  C                 Laboratories, 1982.
  C***ROUTINES CALLED  (NONE)
  C***REVISION HISTORY  (YYMMDD)
  C   790801  DATE WRITTEN
  C   861211  REVISION DATE from Version 3.2
  C   891214  Prologue converted to Version 4.0 format.  (BAB)
  C   900206  Routine changed from user-callable to subsidiary.  (WRB)
  C   900510  Changed calling sequence to delete length of character
  C           and changed routine name from XERABT to XERHLT.  (RWC)
  C   920501  Reformatted the REFERENCES section.  (WRB)
  C***END PROLOGUE  XERHLT
        CHARACTER*(*) MESSG
  C***FIRST EXECUTABLE STATEMENT  XERHLT
        STOP
        END
  