fix import/export of ierode common on Windows
[scilab.git] / scilab / modules / differential_equations / src / fortran / ddasrt.f
1       SUBROUTINE DDASRT (RES,NEQ,T,Y,YPRIME,TOUT,
2      *  INFO,RTOL,ATOL,IDID,RWORK,LRW,IWORK,LIW,RPAR,IPAR,JAC,
3      *  G,NG,JROOT)
4 C***MODIF
5 C   WHEN A ROOT IS FOUND YPRIME WAS NOT UPDATED. see c*SS* modifications
6 C
7 C***BEGIN PROLOGUE  DDASRT
8 C***DATE WRITTEN   821001   (YYMMDD)
9 C***REVISION DATE  910624   (YYMMDD)
10 C***KEYWORDS  DIFFERENTIAL/ALGEBRAIC,BACKWARD DIFFERENTIATION FORMULAS
11 C             IMPLICIT DIFFERENTIAL SYSTEMS
12 C***AUTHOR  PETZOLD,LINDA R.,COMPUTING AND MATHEMATICS RESEARCH DIVISION
13 C             LAWRENCE LIVERMORE NATIONAL LABORATORY
14 C             L - 316, P.O. Box 808,
15 C             LIVERMORE, CA.    94550
16 C***PURPOSE  This code solves a system of differential/algebraic
17 C            equations of the form F(T,Y,YPRIME) = 0.
18 C***DESCRIPTION
19 C
20 C *Usage:
21 C
22 C      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
23 C      EXTERNAL RES, JAC, G
24 C      INTEGER NEQ, INFO(N), IDID, LRW, LIW, IWORK(LIW), IPAR, NG,
25 C     *   JROOT(NG)
26 C      DOUBLE PRECISION T, Y(NEQ), YPRIME(NEQ), TOUT, RTOL, ATOL,
27 C     *   RWORK(LRW), RPAR
28 C
29 C      CALL DDASRT (RES, NEQ, T, Y, YPRIME, TOUT, INFO, RTOL, ATOL,
30 C     *   IDID, RWORK, LRW, IWORK, LIW, RPAR, IPAR, JAC)
31 C
32 C
33 C
34 C *Arguments:
35 C
36 C  RES:EXT  This is a subroutine which you provide to define the
37 C           differential/algebraic system.
38 C
39 C  NEQ:IN  This is the number of equations to be solved.
40 C
41 C  T:INOUT  This is the current value of the independent variable.
42 C
43 C  Y(*):INOUT  This array contains the solution components at T.
44 C
45 C  YPRIME(*):INOUT  This array contains the derivatives of the solution
46 C                   components at T.
47 C
48 C  TOUT:IN  This is a point at which a solution is desired.
49 C
50 C  INFO(N):IN  The basic task of the code is to solve the system from T
51 C              to TOUT and return an answer at TOUT.  INFO is an integer
52 C              array which is used to communicate exactly how you want
53 C              this task to be carried out.  N must be greater than or
54 C              equal to 15.
55 C
56 C  RTOL,ATOL:INOUT  These quantities represent absolute and relative
57 C                   error tolerances which you provide to indicate how
58 C                   accurately you wish the solution to be computed.
59 C                   You may choose them to be both scalars or else
60 C                   both vectors.
61 C
62 C  IDID:OUT  This scalar quantity is an indicator reporting what the
63 C            code did.  You must monitor this integer variable to decide
64 C            what action to take next.
65 C
66 C  RWORK:WORK  A real work array of length LRW which provides the
67 C               code with needed storage space.
68 C
69 C  LRW:IN  The length of RWORK.
70 C
71 C  IWORK:WORK  An integer work array of length LIW which probides the
72 C               code with needed storage space.
73 C
74 C  LIW:IN  The length of IWORK.
75 C
76 C  RPAR,IPAR:IN  These are real and integer parameter arrays which
77 C                you can use for communication between your calling
78 C                program and the RES subroutine (and the JAC subroutine)
79 C
80 C  JAC:EXT  This is the name of a subroutine which you may choose to
81 C           provide for defining a matrix of partial derivatives
82 C           described below.
83 C
84 C  G  This is the name of the subroutine for defining
85 C     constraint functions, G(T,Y), whose roots are desired
86 C     during the integration.  This name must be declared
87 C     external in the calling program.
88 C
89 C  NG  This is the number of constraint functions G(I).
90 C      If there are none, set NG=0, and pass a dummy name
91 C      for G.
92 C
93 C  JROOT  This is an integer array of length NG for output
94 C         of root information.
95 C
96 C
97 C *Description
98 C
99 C  QUANTITIES WHICH MAY BE ALTERED BY THE CODE ARE
100 C     T,Y(*),YPRIME(*),INFO(1),RTOL,ATOL,
101 C     IDID,RWORK(*) AND IWORK(*).
102 C
103 C  Subroutine DDASRT uses the backward differentiation formulas of
104 C  orders one through five to solve a system of the above form for Y and
105 C  YPRIME.  Values for Y and YPRIME at the initial time must be given as
106 C  input.  These values must be consistent, (that is, if T,Y,YPRIME are
107 C  the given initial values, they must satisfy F(T,Y,YPRIME) = 0.).  The
108 C  subroutine solves the system from T to TOUT.
109 C  It is easy to continue the solution to get results at additional
110 C  TOUT.  This is the interval mode of operation.  Intermediate results
111 C  can also be obtained easily by using the intermediate-output
112 C  capability.  If DDASRT detects a sign-change in G(T,Y), then
113 C  it will return the intermediate value of T and Y for which
114 C  G(T,Y) = 0.
115 C
116 C  ---------INPUT-WHAT TO DO ON THE FIRST CALL TO DDASRT---------------
117 C
118 C
119 C  The first call of the code is defined to be the start of each new
120 C  problem. Read through the descriptions of all the following items,
121 C  provide sufficient storage space for designated arrays, set
122 C  appropriate variables for the initialization of the problem, and
123 C  give information about how you want the problem to be solved.
124 C
125 C
126 C  RES -- Provide a subroutine of the form
127 C             SUBROUTINE RES(T,Y,YPRIME,DELTA,IRES,RPAR,IPAR)
128 C         to define the system of differential/algebraic
129 C         equations which is to be solved. For the given values
130 C         of T,Y and YPRIME, the subroutine should
131 C         return the residual of the differential/algebraic
132 C         system
133 C             DELTA = F(T,Y,YPRIME)
134 C         (DELTA(*) is a vector of length NEQ which is
135 C         output for RES.)
136 C
137 C         Subroutine RES must not alter T,Y or YPRIME.
138 C         You must declare the name RES in an external
139 C         statement in your program that calls DDASRT.
140 C         You must dimension Y,YPRIME and DELTA in RES.
141 C
142 C         IRES is an integer flag which is always equal to
143 C         zero on input. Subroutine RES should alter IRES
144 C         only if it encounters an illegal value of Y or
145 C         a stop condition. Set IRES = -1 if an input value
146 C         is illegal, and DDASRT will try to solve the problem
147 C         without getting IRES = -1. If IRES = -2, DDASRT
148 C         will return control to the calling program
149 C         with IDID = -11.
150 C
151 C         RPAR and IPAR are real and integer parameter arrays which
152 C         you can use for communication between your calling program
153 C         and subroutine RES. They are not altered by DDASRT. If you
154 C         do not need RPAR or IPAR, ignore these parameters by treat-
155 C         ing them as dummy arguments. If you do choose to use them,
156 C         dimension them in your calling program and in RES as arrays
157 C         of appropriate length.
158 C
159 C  NEQ -- Set it to the number of differential equations.
160 C         (NEQ .GE. 1)
161 C
162 C  T -- Set it to the initial point of the integration.
163 C       T must be defined as a variable.
164 C
165 C  Y(*) -- Set this vector to the initial values of the NEQ solution
166 C          components at the initial point. You must dimension Y of
167 C          length at least NEQ in your calling program.
168 C
169 C  YPRIME(*) -- Set this vector to the initial values of
170 C               the NEQ first derivatives of the solution
171 C               components at the initial point. You
172 C               must dimension YPRIME at least NEQ
173 C               in your calling program. If you do not
174 C               know initial values of some of the solution
175 C               components, see the explanation of INFO(11).
176 C
177 C  TOUT - Set it to the first point at which a solution
178 C         is desired. You can not take TOUT = T.
179 C         integration either forward in T (TOUT .GT. T) or
180 C         backward in T (TOUT .LT. T) is permitted.
181 C
182 C         The code advances the solution from T to TOUT using
183 C         step sizes which are automatically selected so as to
184 C         achieve the desired accuracy. If you wish, the code will
185 C         return with the solution and its derivative at
186 C         intermediate steps (intermediate-output mode) so that
187 C         you can monitor them, but you still must provide TOUT in
188 C         accord with the basic aim of the code.
189 C
190 C         the first step taken by the code is a critical one
191 C         because it must reflect how fast the solution changes near
192 C         the initial point. The code automatically selects an
193 C         initial step size which is practically always suitable for
194 C         the problem. By using the fact that the code will not step
195 C         past TOUT in the first step, you could, if necessary,
196 C         restrict the length of the initial step size.
197 C
198 C         For some problems it may not be permissable to integrate
199 C         past a point TSTOP because a discontinuity occurs there
200 C         or the solution or its derivative is not defined beyond
201 C         TSTOP. When you have declared a TSTOP point (SEE INFO(4)
202 C         and RWORK(1)), you have told the code not to integrate
203 C         past TSTOP. In this case any TOUT beyond TSTOP is invalid
204 C         input.
205 C
206 C  INFO(*) - Use the INFO array to give the code more details about
207 C            how you want your problem solved. This array should be
208 C            dimensioned of length 15, though DDASRT uses
209 C            only the first eleven entries. You must respond to all of
210 C            the following items which are arranged as questions. The
211 C            simplest use of the code corresponds to answering all
212 C            questions as yes, i.e. setting all entries of INFO to 0.
213 C
214 C       INFO(1) - This parameter enables the code to initialize
215 C              itself. You must set it to indicate the start of every
216 C              new problem.
217 C
218 C          **** Is this the first call for this problem ...
219 C                Yes - Set INFO(1) = 0
220 C                 No - Not applicable here.
221 C                      See below for continuation calls.  ****
222 C
223 C       INFO(2) - How much accuracy you want of your solution
224 C              is specified by the error tolerances RTOL and ATOL.
225 C              The simplest use is to take them both to be scalars.
226 C              To obtain more flexibility, they can both be vectors.
227 C              The code must be told your choice.
228 C
229 C          **** Are both error tolerances RTOL, ATOL scalars ...
230 C                Yes - Set INFO(2) = 0
231 C                      and input scalars for both RTOL and ATOL
232 C                 No - Set INFO(2) = 1
233 C                      and input arrays for both RTOL and ATOL ****
234 C
235 C       INFO(3) - The code integrates from T in the direction
236 C              of TOUT by steps. If you wish, it will return the
237 C              computed solution and derivative at the next
238 C              intermediate step (the intermediate-output mode) or
239 C              TOUT, whichever comes first. This is a good way to
240 C              proceed if you want to see the behavior of the solution.
241 C              If you must have solutions at a great many specific
242 C              TOUT points, this code will compute them efficiently.
243 C
244 C          **** Do you want the solution only at
245 C                TOUT (and not at the next intermediate step) ...
246 C                 Yes - Set INFO(3) = 0
247 C                  No - Set INFO(3) = 1 ****
248 C
249 C       INFO(4) - To handle solutions at a great many specific
250 C              values TOUT efficiently, this code may integrate past
251 C              TOUT and interpolate to obtain the result at TOUT.
252 C              Sometimes it is not possible to integrate beyond some
253 C              point TSTOP because the equation changes there or it is
254 C              not defined past TSTOP. Then you must tell the code
255 C              not to go past.
256 C
257 C           **** Can the integration be carried out without any
258 C                restrictions on the independent variable T ...
259 C                 Yes - Set INFO(4)=0
260 C                  No - Set INFO(4)=1
261 C                       and define the stopping point TSTOP by
262 C                       setting RWORK(1)=TSTOP ****
263 C
264 C       INFO(5) - To solve differential/algebraic problems it is
265 C              necessary to use a matrix of partial derivatives of the
266 C              system of differential equations. If you do not
267 C              provide a subroutine to evaluate it analytically (see
268 C              description of the item JAC in the call list), it will
269 C              be approximated by numerical differencing in this code.
270 C              although it is less trouble for you to have the code
271 C              compute partial derivatives by numerical differencing,
272 C              the solution will be more reliable if you provide the
273 C              derivatives via JAC. Sometimes numerical differencing
274 C              is cheaper than evaluating derivatives in JAC and
275 C              sometimes it is not - this depends on your problem.
276 C
277 C           **** Do you want the code to evaluate the partial
278 C                derivatives automatically by numerical differences ...
279 C                   Yes - Set INFO(5)=0
280 C                    No - Set INFO(5)=1
281 C                  and provide subroutine JAC for evaluating the
282 C                  matrix of partial derivatives ****
283 C
284 C       INFO(6) - DDASRT will perform much better if the matrix of
285 C              partial derivatives, DG/DY + CJ*DG/DYPRIME,
286 C              (here CJ is a scalar determined by DDASRT)
287 C              is banded and the code is told this. In this
288 C              case, the storage needed will be greatly reduced,
289 C              numerical differencing will be performed much cheaper,
290 C              and a number of important algorithms will execute much
291 C              faster. The differential equation is said to have
292 C              half-bandwidths ML (lower) and MU (upper) if equation i
293 C              involves only unknowns Y(J) with
294 C                             I-ML .LE. J .LE. I+MU
295 C              for all I=1,2,...,NEQ. Thus, ML and MU are the widths
296 C              of the lower and upper parts of the band, respectively,
297 C              with the main diagonal being excluded. If you do not
298 C              indicate that the equation has a banded matrix of partial
299 C              derivatives, the code works with a full matrix of NEQ**2
300 C              elements (stored in the conventional way). Computations
301 C              with banded matrices cost less time and storage than with
302 C              full matrices if 2*ML+MU .LT. NEQ. If you tell the
303 C              code that the matrix of partial derivatives has a banded
304 C              structure and you want to provide subroutine JAC to
305 C              compute the partial derivatives, then you must be careful
306 C              to store the elements of the matrix in the special form
307 C              indicated in the description of JAC.
308 C
309 C          **** Do you want to solve the problem using a full
310 C               (dense) matrix (and not a special banded
311 C               structure) ...
312 C                Yes - Set INFO(6)=0
313 C                 No - Set INFO(6)=1
314 C                       and provide the lower (ML) and upper (MU)
315 C                       bandwidths by setting
316 C                       IWORK(1)=ML
317 C                       IWORK(2)=MU ****
318 C
319 C
320 C        INFO(7) -- You can specify a maximum (absolute value of)
321 C              stepsize, so that the code
322 C              will avoid passing over very
323 C              large regions.
324 C
325 C          ****  Do you want the code to decide
326 C                on its own maximum stepsize?
327 C                Yes - Set INFO(7)=0
328 C                 No - Set INFO(7)=1
329 C                      and define HMAX by setting
330 C                      RWORK(2)=HMAX ****
331 C
332 C        INFO(8) -- Differential/algebraic problems
333 C              may occaisionally suffer from
334 C              severe scaling difficulties on the
335 C              first step. If you know a great deal
336 C              about the scaling of your problem, you can
337 C              help to alleviate this problem by
338 C              specifying an initial stepsize H0.
339 C
340 C          ****  Do you want the code to define
341 C                its own initial stepsize?
342 C                Yes - Set INFO(8)=0
343 C                 No - Set INFO(8)=1
344 C                      and define H0 by setting
345 C                      RWORK(3)=H0 ****
346 C
347 C        INFO(9) -- If storage is a severe problem,
348 C              you can save some locations by
349 C              restricting the maximum order MAXORD.
350 C              the default value is 5. for each
351 C              order decrease below 5, the code
352 C              requires NEQ fewer locations, however
353 C              it is likely to be slower. In any
354 C              case, you must have 1 .LE. MAXORD .LE. 5
355 C          ****  Do you want the maximum order to
356 C                default to 5?
357 C                Yes - Set INFO(9)=0
358 C                 No - Set INFO(9)=1
359 C                      and define MAXORD by setting
360 C                      IWORK(3)=MAXORD ****
361 C
362 C        INFO(10) --If you know that the solutions to your equations
363 C               will always be nonnegative, it may help to set this
364 C               parameter. However, it is probably best to
365 C               try the code without using this option first,
366 C               and only to use this option if that doesn't
367 C               work very well.
368 C           ****  Do you want the code to solve the problem without
369 C                 invoking any special nonnegativity constraints?
370 C                  Yes - Set INFO(10)=0
371 C                   No - Set INFO(10)=1
372 C
373 C        INFO(11) --DDASRT normally requires the initial T,
374 C               Y, and YPRIME to be consistent. That is,
375 C               you must have F(T,Y,YPRIME) = 0 at the initial
376 C               time. If you do not know the initial
377 C               derivative precisely, you can let DDASRT try
378 C               to compute it.
379 C          ****   Are the initial T, Y, YPRIME consistent?
380 C                 Yes - Set INFO(11) = 0
381 C                  No - Set INFO(11) = 1,
382 C                       and set YPRIME to an initial approximation
383 C                       to YPRIME.  (If you have no idea what
384 C                       YPRIME should be, set it to zero. Note
385 C                       that the initial Y should be such
386 C                       that there must exist a YPRIME so that
387 C                       F(T,Y,YPRIME) = 0.)
388 C
389 C   RTOL, ATOL -- You must assign relative (RTOL) and absolute (ATOL
390 C               error tolerances to tell the code how accurately you
391 C               want the solution to be computed. They must be defined
392 C               as variables because the code may change them. You
393 C               have two choices --
394 C                     Both RTOL and ATOL are scalars. (INFO(2)=0)
395 C                     Both RTOL and ATOL are vectors. (INFO(2)=1)
396 C               in either case all components must be non-negative.
397 C
398 C               The tolerances are used by the code in a local error
399 C               test at each step which requires roughly that
400 C                     ABS(LOCAL ERROR) .LE. RTOL*ABS(Y)+ATOL
401 C               for each vector component.
402 C               (More specifically, a root-mean-square norm is used to
403 C               measure the size of vectors, and the error test uses the
404 C               magnitude of the solution at the beginning of the step.)
405 C
406 C               The true (global) error is the difference between the
407 C               true solution of the initial value problem and the
408 C               computed approximation. Practically all present day
409 C               codes, including this one, control the local error at
410 C               each step and do not even attempt to control the global
411 C               error directly.
412 C               Usually, but not always, the true accuracy of the
413 C               computed Y is comparable to the error tolerances. This
414 C               code will usually, but not always, deliver a more
415 C               accurate solution if you reduce the tolerances and
416 C               integrate again. By comparing two such solutions you
417 C               can get a fairly reliable idea of the true error in the
418 C               solution at the bigger tolerances.
419 C
420 C               Setting ATOL=0. results in a pure relative error test on
421 C               that component. Setting RTOL=0. results in a pure
422 C               absolute error test on that component. A mixed test
423 C               with non-zero RTOL and ATOL corresponds roughly to a
424 C               relative error test when the solution component is much
425 C               bigger than ATOL and to an absolute error test when the
426 C               solution component is smaller than the threshhold ATOL.
427 C
428 C               The code will not attempt to compute a solution at an
429 C               accuracy unreasonable for the machine being used. It
430 C               will advise you if you ask for too much accuracy and
431 C               inform you as to the maximum accuracy it believes
432 C               possible.
433 C
434 C  RWORK(*) --  Dimension this real work array of length LRW in your
435 C               calling program.
436 C
437 C  LRW -- Set it to the declared length of the RWORK array.
438 C               You must have
439 C                    LRW .GE. 50+(MAXORD+4)*NEQ+NEQ**2+3*NG
440 C               for the full (dense) JACOBIAN case (when INFO(6)=0), or
441 C                    LRW .GE. 50+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQ+3*NG
442 C               for the banded user-defined JACOBIAN case
443 C               (when INFO(5)=1 and INFO(6)=1), or
444 C                     LRW .GE. 50+(MAXORD+4)*NEQ+(2*ML+MU+1)*NEQ
445 C                           +2*(NEQ/(ML+MU+1)+1)+3*NG
446 C               for the banded finite-difference-generated JACOBIAN case
447 C               (when INFO(5)=0 and INFO(6)=1)
448 C
449 C  IWORK(*) --  Dimension this integer work array of length LIW in
450 C               your calling program.
451 C
452 C  LIW -- Set it to the declared length of the IWORK array.
453 C               you must have LIW .GE. 20+NEQ
454 C
455 C  RPAR, IPAR -- These are parameter arrays, of real and integer
456 C               type, respectively. You can use them for communication
457 C               between your program that calls DDASRT and the
458 C               RES subroutine (and the JAC subroutine). They are not
459 C               altered by DDASRT. If you do not need RPAR or IPAR,
460 C               ignore these parameters by treating them as dummy
461 C               arguments. If you do choose to use them, dimension
462 C               them in your calling program and in RES (and in JAC)
463 C               as arrays of appropriate length.
464 C
465 C  JAC -- If you have set INFO(5)=0, you can ignore this parameter
466 C               by treating it as a dummy argument. Otherwise, you must
467 C               provide a subroutine of the form
468 C               JAC(T,Y,YPRIME,PD,CJ,RPAR,IPAR)
469 C               to define the matrix of partial derivatives
470 C               PD=DG/DY+CJ*DG/DYPRIME
471 C               CJ is a scalar which is input to JAC.
472 C               For the given values of T,Y,YPRIME, the
473 C               subroutine must evaluate the non-zero partial
474 C               derivatives for each equation and each solution
475 C               component, and store these values in the
476 C               matrix PD. The elements of PD are set to zero
477 C               before each call to JAC so only non-zero elements
478 C               need to be defined.
479 C
480 C               Subroutine JAC must not alter T,Y,(*),YPRIME(*), or CJ.
481 C               You must declare the name JAC in an
482 C               EXTERNAL STATEMENT in your program that calls
483 C               DDASRT. You must dimension Y, YPRIME and PD
484 C               in JAC.
485 C
486 C               The way you must store the elements into the PD matrix
487 C               depends on the structure of the matrix which you
488 C               indicated by INFO(6).
489 C               *** INFO(6)=0 -- Full (dense) matrix ***
490 C                   Give PD a first dimension of NEQ.
491 C                   When you evaluate the (non-zero) partial derivative
492 C                   of equation I with respect to variable J, you must
493 C                   store it in PD according to
494 C                   PD(I,J) = * DF(I)/DY(J)+CJ*DF(I)/DYPRIME(J)*
495 C               *** INFO(6)=1 -- Banded JACOBIAN with ML lower and MU
496 C                   upper diagonal bands (refer to INFO(6) description
497 C                   of ML and MU) ***
498 C                   Give PD a first dimension of 2*ML+MU+1.
499 C                   when you evaluate the (non-zero) partial derivative
500 C                   of equation I with respect to variable J, you must
501 C                   store it in PD according to
502 C                   IROW = I - J + ML + MU + 1
503 C                   PD(IROW,J) = *DF(I)/DY(J)+CJ*DF(I)/DYPRIME(J)*
504 C               RPAR and IPAR are real and integer parameter arrays
505 C               which you can use for communication between your calling
506 C               program and your JACOBIAN subroutine JAC. They are not
507 C               altered by DDASRT. If you do not need RPAR or IPAR,
508 C               ignore these parameters by treating them as dummy
509 C               arguments. If you do choose to use them, dimension
510 C               them in your calling program and in JAC as arrays of
511 C               appropriate length.
512 C
513 C  G -- This is the name of the subroutine for defining constraint
514 C               functions, whose roots are desired during the
515 C               integration.  It is to have the form
516 C                   SUBROUTINE G(NEQ,T,Y,NG,GOUT,RPAR,IPAR)
517 C                   DIMENSION Y(NEQ),GOUT(NG),
518 C               where NEQ, T, Y and NG are INPUT, and the array GOUT is
519 C               output.  NEQ, T, and Y have the same meaning as in the
520 C               RES routine, and GOUT is an array of length NG.
521 C               For I=1,...,NG, this routine is to load into GOUT(I)
522 C               the value at (T,Y) of the I-th constraint function G(I).
523 C               DDASRT will find roots of the G(I) of odd multiplicity
524 C               (that is, sign changes) as they occur during
525 C               the integration.  G must be declared EXTERNAL in the
526 C               calling program.
527 C
528 C               CAUTION..because of numerical errors in the functions
529 C               G(I) due to roundoff and integration error, DDASRT
530 C               may return false roots, or return the same root at two
531 C               or more nearly equal values of T.  If such false roots
532 C               are suspected, the user should consider smaller error
533 C               tolerances and/or higher precision in the evaluation of
534 C               the G(I).
535 C
536 C               If a root of some G(I) defines the end of the problem,
537 C               the input to DDASRT should nevertheless allow
538 C               integration to a point slightly past that ROOT, so
539 C               that DDASRT can locate the root by interpolation.
540 C
541 C  NG -- The number of constraint functions G(I).  If there are none,
542 C               set NG = 0, and pass a dummy name for G.
543 C
544 C JROOT -- This is an integer array of length NG.  It is used only for
545 C               output.  On a return where one or more roots have been
546 C               found, JROOT(I)=1 If G(I) has a root at T,
547 C               or JROOT(I)=0 if not.
548 C
549 C
550 C
551 C  OPTIONALLY REPLACEABLE NORM ROUTINE:
552 C  DDASRT uses a weighted norm DDANRM to measure the size
553 C  of vectors such as the estimated error in each step.
554 C  A FUNCTION subprogram
555 C    DOUBLE PRECISION FUNCTION DDANRM(NEQ,V,WT,RPAR,IPAR)
556 C    DIMENSION V(NEQ),WT(NEQ)
557 C  is used to define this norm. Here, V is the vector
558 C  whose norm is to be computed, and WT is a vector of
559 C  weights.  A DDANRM routine has been included with DDASRT
560 C  which computes the weighted root-mean-square norm
561 C  given by
562 C    DDANRM=SQRT((1/NEQ)*SUM(V(I)/WT(I))**2)
563 C  this norm is suitable for most problems. In some
564 C  special cases, it may be more convenient and/or
565 C  efficient to define your own norm by writing a function
566 C  subprogram to be called instead of DDANRM. This should
567 C  ,however, be attempted only after careful thought and
568 C  consideration.
569 C
570 C
571 C------OUTPUT-AFTER ANY RETURN FROM DDASRT----
572 C
573 C  The principal aim of the code is to return a computed solution at
574 C  TOUT, although it is also possible to obtain intermediate results
575 C  along the way. To find out whether the code achieved its goal
576 C  or if the integration process was interrupted before the task was
577 C  completed, you must check the IDID parameter.
578 C
579 C
580 C   T -- The solution was successfully advanced to the
581 C               output value of T.
582 C
583 C   Y(*) -- Contains the computed solution approximation at T.
584 C
585 C   YPRIME(*) -- Contains the computed derivative
586 C               approximation at T.
587 C
588 C   IDID -- Reports what the code did.
589 C
590 C                     *** Task completed ***
591 C                Reported by positive values of IDID
592 C
593 C           IDID = 1 -- A step was successfully taken in the
594 C                   intermediate-output mode. The code has not
595 C                   yet reached TOUT.
596 C
597 C           IDID = 2 -- The integration to TSTOP was successfully
598 C                   completed (T=TSTOP) by stepping exactly to TSTOP.
599 C
600 C           IDID = 3 -- The integration to TOUT was successfully
601 C                   completed (T=TOUT) by stepping past TOUT.
602 C                   Y(*) is obtained by interpolation.
603 C                   YPRIME(*) is obtained by interpolation.
604 C
605 C           IDID = 4 -- The integration was successfully completed
606 C                   by finding one or more roots of G at T.
607 C
608 C                    *** Task interrupted ***
609 C                Reported by negative values of IDID
610 C
611 C           IDID = -1 -- A large amount of work has been expended.
612 C                   (About 500 steps)
613 C
614 C           IDID = -2 -- The error tolerances are too stringent.
615 C
616 C           IDID = -3 -- The local error test cannot be satisfied
617 C                   because you specified a zero component in ATOL
618 C                   and the corresponding computed solution
619 C                   component is zero. Thus, a pure relative error
620 C                   test is impossible for this component.
621 C
622 C           IDID = -6 -- DDASRT had repeated error test
623 C                   failures on the last attempted step.
624 C
625 C           IDID = -7 -- The corrector could not converge.
626 C
627 C           IDID = -8 -- The matrix of partial derivatives
628 C                   is singular.
629 C
630 C           IDID = -9 -- The corrector could not converge.
631 C                   there were repeated error test failures
632 C                   in this step.
633 C
634 C           IDID =-10 -- The corrector could not converge
635 C                   because IRES was equal to minus one.
636 C
637 C           IDID =-11 -- IRES equal to -2 was encountered
638 C                   and control is being returned to the
639 C                   calling program.
640 C
641 C           IDID =-12 -- DDASRT failed to compute the initial
642 C                   YPRIME.
643 C
644 C
645 C
646 C           IDID = -13,..,-32 -- Not applicable for this code
647 C
648 C                    *** Task terminated ***
649 C                Reported by the value of IDID=-33
650 C
651 C           IDID = -33 -- The code has encountered trouble from which
652 C                   it cannot recover. A message is printed
653 C                   explaining the trouble and control is returned
654 C                   to the calling program. For example, this occurs
655 C                   when invalid input is detected.
656 C
657 C   RTOL, ATOL -- These quantities remain unchanged except when
658 C               IDID = -2. In this case, the error tolerances have been
659 C               increased by the code to values which are estimated to
660 C               be appropriate for continuing the integration. However,
661 C               the reported solution at T was obtained using the input
662 C               values of RTOL and ATOL.
663 C
664 C   RWORK, IWORK -- Contain information which is usually of no
665 C               interest to the user but necessary for subsequent calls.
666 C               However, you may find use for
667 C
668 C               RWORK(3)--Which contains the step size H to be
669 C                       attempted on the next step.
670 C
671 C               RWORK(4)--Which contains the current value of the
672 C                       independent variable, i.e., the farthest point
673 C                       integration has reached. This will be different
674 C                       from T only when interpolation has been
675 C                       performed (IDID=3).
676 C
677 C               RWORK(7)--Which contains the stepsize used
678 C                       on the last successful step.
679 C
680 C               IWORK(7)--Which contains the order of the method to
681 C                       be attempted on the next step.
682 C
683 C               IWORK(8)--Which contains the order of the method used
684 C                       on the last step.
685 C
686 C               IWORK(11)--Which contains the number of steps taken so
687 C                        far.
688 C
689 C               IWORK(12)--Which contains the number of calls to RES
690 C                        so far.
691 C
692 C               IWORK(13)--Which contains the number of evaluations of
693 C                        the matrix of partial derivatives needed so
694 C                        far.
695 C
696 C               IWORK(14)--Which contains the total number
697 C                        of error test failures so far.
698 C
699 C               IWORK(15)--Which contains the total number
700 C                        of convergence test failures so far.
701 C                        (includes singular iteration matrix
702 C                        failures.)
703 C
704 C               IWORK(16)--Which contains the total number of calls
705 C                        to the constraint function g so far
706 C
707 C
708 C
709 C   INPUT -- What to do to continue the integration
710 C            (calls after the first)                **
711 C
712 C     This code is organized so that subsequent calls to continue the
713 C     integration involve little (if any) additional effort on your
714 C     part. You must monitor the IDID parameter in order to determine
715 C     what to do next.
716 C
717 C     Recalling that the principal task of the code is to integrate
718 C     from T to TOUT (the interval mode), usually all you will need
719 C     to do is specify a new TOUT upon reaching the current TOUT.
720 C
721 C     Do not alter any quantity not specifically permitted below,
722 C     in particular do not alter NEQ,T,Y(*),YPRIME(*),RWORK(*),IWORK(*)
723 C     or the differential equation in subroutine RES. Any such
724 C     alteration constitutes a new problem and must be treated as such,
725 C     i.e., you must start afresh.
726 C
727 C     You cannot change from vector to scalar error control or vice
728 C     versa (INFO(2)), but you can change the size of the entries of
729 C     RTOL, ATOL. Increasing a tolerance makes the equation easier
730 C     to integrate. Decreasing a tolerance will make the equation
731 C     harder to integrate and should generally be avoided.
732 C
733 C     You can switch from the intermediate-output mode to the
734 C     interval mode (INFO(3)) or vice versa at any time.
735 C
736 C     If it has been necessary to prevent the integration from going
737 C     past a point TSTOP (INFO(4), RWORK(1)), keep in mind that the
738 C     code will not integrate to any TOUT beyond the currently
739 C     specified TSTOP. Once TSTOP has been reached you must change
740 C     the value of TSTOP or set INFO(4)=0. You may change INFO(4)
741 C     or TSTOP at any time but you must supply the value of TSTOP in
742 C     RWORK(1) whenever you set INFO(4)=1.
743 C
744 C     Do not change INFO(5), INFO(6), IWORK(1), or IWORK(2)
745 C     unless you are going to restart the code.
746 C
747 C                    *** Following a completed task ***
748 C     If
749 C     IDID = 1, call the code again to continue the integration
750 C                  another step in the direction of TOUT.
751 C
752 C     IDID = 2 or 3, define a new TOUT and call the code again.
753 C                  TOUT must be different from T. You cannot change
754 C                  the direction of integration without restarting.
755 C
756 C     IDID = 4, call the code again to continue the integration
757 C                  another step in the direction of TOUT.  You may
758 C                  change the functions in G after a return with IDID=4,
759 C                  but the number of constraint functions NG must remain
760 C                  the same.  If you wish to change
761 C                  the functions in RES or in G, then you
762 C                  must restart the code.
763 C
764 C                    *** Following an interrupted task ***
765 C                  To show the code that you realize the task was
766 C                  interrupted and that you want to continue, you
767 C                  must take appropriate action and set INFO(1) = 1
768 C     If
769 C     IDID = -1, The code has taken about 500 steps.
770 C                  If you want to continue, set INFO(1) = 1 and
771 C                  call the code again. An additional 500 steps
772 C                  will be allowed.
773 C
774 C     IDID = -2, The error tolerances RTOL, ATOL have been
775 C                  increased to values the code estimates appropriate
776 C                  for continuing. You may want to change them
777 C                  yourself. If you are sure you want to continue
778 C                  with relaxed error tolerances, set INFO(1)=1 and
779 C                  call the code again.
780 C
781 C     IDID = -3, A solution component is zero and you set the
782 C                  corresponding component of ATOL to zero. If you
783 C                  are sure you want to continue, you must first
784 C                  alter the error criterion to use positive values
785 C                  for those components of ATOL corresponding to zero
786 C                  solution components, then set INFO(1)=1 and call
787 C                  the code again.
788 C
789 C     IDID = -4,-5  --- Cannot occur with this code.
790 C
791 C     IDID = -6, Repeated error test failures occurred on the
792 C                  last attempted step in DDASRT. A singularity in the
793 C                  solution may be present. If you are absolutely
794 C                  certain you want to continue, you should restart
795 C                  the integration. (Provide initial values of Y and
796 C                  YPRIME which are consistent)
797 C
798 C     IDID = -7, Repeated convergence test failures occurred
799 C                  on the last attempted step in DDASRT. An inaccurate
800 C                  or ill-conditioned JACOBIAN may be the problem. If
801 C                  you are absolutely certain you want to continue, you
802 C                  should restart the integration.
803 C
804 C     IDID = -8, The matrix of partial derivatives is singular.
805 C                  Some of your equations may be redundant.
806 C                  DDASRT cannot solve the problem as stated.
807 C                  It is possible that the redundant equations
808 C                  could be removed, and then DDASRT could
809 C                  solve the problem. It is also possible
810 C                  that a solution to your problem either
811 C                  does not exist or is not unique.
812 C
813 C     IDID = -9, DDASRT had multiple convergence test
814 C                  failures, preceeded by multiple error
815 C                  test failures, on the last attempted step.
816 C                  It is possible that your problem
817 C                  is ill-posed, and cannot be solved
818 C                  using this code. Or, there may be a
819 C                  discontinuity or a singularity in the
820 C                  solution. If you are absolutely certain
821 C                  you want to continue, you should restart
822 C                  the integration.
823 C
824 C    IDID =-10, DDASRT had multiple convergence test failures
825 C                  because IRES was equal to minus one.
826 C                  If you are absolutely certain you want
827 C                  to continue, you should restart the
828 C                  integration.
829 C
830 C    IDID =-11, IRES=-2 was encountered, and control is being
831 C                  returned to the calling program.
832 C
833 C    IDID =-12, DDASRT failed to compute the initial YPRIME.
834 C               This could happen because the initial
835 C               approximation to YPRIME was not very good, or
836 C               if a YPRIME consistent with the initial Y
837 C               does not exist. The problem could also be caused
838 C               by an inaccurate or singular iteration matrix.
839 C
840 C
841 C
842 C     IDID = -13,..,-32 --- Cannot occur with this code.
843 C
844 C                       *** Following a terminated task ***
845 C     If IDID= -33, you cannot continue the solution of this
846 C                  problem. An attempt to do so will result in your
847 C                  run being terminated.
848 C
849 C  ---------------------------------------------------------------------
850 C
851 C***REFERENCE
852 C      K. E. Brenan, S. L. Campbell, and L. R. Petzold, Numerical 
853 C      Solution of Initial-Value Problems in Differential-Algebraic
854 C      Equations, Elsevier, New York, 1989.
855 C
856 C***ROUTINES CALLED  DDASTP,DDAINI,DDANRM,DDAWTS,DDATRP,DRCHEK,DROOTS,
857 C                    XERRWV,D1MACH
858 C***END PROLOGUE  DDASRT
859 C
860 C**End
861 C
862       
863       
864       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
865
866       include 'stack.h'
867       
868       LOGICAL DONE
869       EXTERNAL RES, JAC, G
870       DIMENSION Y(*),YPRIME(*)
871       DIMENSION INFO(15)
872       DIMENSION RWORK(*),IWORK(*)
873       DIMENSION RTOL(*),ATOL(*)
874       DIMENSION RPAR(*),IPAR(*)
875       CHARACTER MSG*80
876 C
877 C     SET POINTERS INTO IWORK
878       PARAMETER (LML=1, LMU=2, LMXORD=3, LMTYPE=4, LNST=11,
879      *  LNRE=12, LNJE=13, LETF=14, LCTF=15, LNGE=16, LNPD=17,
880      *  LIRFND=18, LIPVT=21, LJCALC=5, LPHASE=6, LK=7, LKOLD=8,
881      *  LNS=9, LNSTL=10, LIWM=1)
882 C
883 C     SET RELATIVE OFFSET INTO RWORK
884       PARAMETER (NPD=1)
885 C
886 C     SET POINTERS INTO RWORK
887       PARAMETER (LTSTOP=1, LHMAX=2, LH=3, LTN=4,
888      *  LCJ=5, LCJOLD=6, LHOLD=7, LS=8, LROUND=9,
889      *  LALPHA=11, LBETA=17, LGAMMA=23,
890      *  LPSI=29, LSIGMA=35, LT0=41, LTLAST=42, LALPHR=43, LX2=44,
891      *  LDELTA=51)
892 C
893 C***FIRST EXECUTABLE STATEMENT  DDASRT
894      
895       IF(INFO(1).NE.0)GO TO 100
896 C
897 C-----------------------------------------------------------------------
898 C     THIS BLOCK IS EXECUTED FOR THE INITIAL CALL ONLY.
899 C     IT CONTAINS CHECKING OF INPUTS AND INITIALIZATIONS.
900 C-----------------------------------------------------------------------
901 C
902 C     FIRST CHECK INFO ARRAY TO MAKE SURE ALL ELEMENTS OF INFO
903 C     ARE EITHER ZERO OR ONE.
904       DO 10 I=2,11
905          IF(INFO(I).NE.0.AND.INFO(I).NE.1)GO TO 701
906 10       CONTINUE
907 C
908       IF(NEQ.LE.0)GO TO 702
909 C
910 C     CHECK AND COMPUTE MAXIMUM ORDER
911       MXORD=5
912       IF(INFO(9).EQ.0)GO TO 20
913          MXORD=IWORK(LMXORD)
914          IF(MXORD.LT.1.OR.MXORD.GT.5)GO TO 703
915 20       IWORK(LMXORD)=MXORD
916 C
917 C     COMPUTE MTYPE,LENPD,LENRW.CHECK ML AND MU.
918       IF(INFO(6).NE.0)GO TO 40
919          LENPD=NEQ**2
920          LENRW=50+(IWORK(LMXORD)+4)*NEQ+LENPD
921          IF(INFO(5).NE.0)GO TO 30
922             IWORK(LMTYPE)=2
923             GO TO 60
924 30          IWORK(LMTYPE)=1
925             GO TO 60
926 40    IF(IWORK(LML).LT.0.OR.IWORK(LML).GE.NEQ)GO TO 717
927       IF(IWORK(LMU).LT.0.OR.IWORK(LMU).GE.NEQ)GO TO 718
928       LENPD=(2*IWORK(LML)+IWORK(LMU)+1)*NEQ
929       IF(INFO(5).NE.0)GO TO 50
930          IWORK(LMTYPE)=5
931          MBAND=IWORK(LML)+IWORK(LMU)+1
932          MSAVE=(NEQ/MBAND)+1
933          LENRW=50+(IWORK(LMXORD)+4)*NEQ+LENPD+2*MSAVE
934          GO TO 60
935 50       IWORK(LMTYPE)=4
936          LENRW=50+(IWORK(LMXORD)+4)*NEQ+LENPD
937 C
938 C     CHECK LENGTHS OF RWORK AND IWORK
939 60    LENIW=20+NEQ
940       IWORK(LNPD)=LENPD
941       IF(LRW.LT.LENRW)GO TO 704
942       IF(LIW.LT.LENIW)GO TO 705
943 C
944 C     CHECK TO SEE THAT TOUT IS DIFFERENT FROM T
945 C     Also check to see that NG is larger than 0.
946       IF(TOUT .EQ. T)GO TO 719
947       IF(NG .LT. 0) GO TO 730
948 C
949 C     CHECK HMAX
950       IF(INFO(7).EQ.0)GO TO 70
951          HMAX=RWORK(LHMAX)
952          IF(HMAX.LE.0.0D0)GO TO 710
953 70    CONTINUE
954 C
955 C     INITIALIZE COUNTERS
956       IWORK(LNST)=0
957       IWORK(LNRE)=0
958       IWORK(LNJE)=0
959       IWORK(LNGE)=0
960 C
961       IWORK(LNSTL)=0
962       IDID=1
963       GO TO 200
964 C
965 C-----------------------------------------------------------------------
966 C     THIS BLOCK IS FOR CONTINUATION CALLS
967 C     ONLY. HERE WE CHECK INFO(1),AND IF THE
968 C     LAST STEP WAS INTERRUPTED WE CHECK WHETHER
969 C     APPROPRIATE ACTION WAS TAKEN.
970 C-----------------------------------------------------------------------
971 C
972 100   CONTINUE
973       IF(INFO(1).EQ.1)GO TO 110
974       IF(INFO(1).NE.-1)GO TO 701
975 C     IF WE ARE HERE, THE LAST STEP WAS INTERRUPTED
976 C     BY AN ERROR CONDITION FROM DDASTP,AND
977 C     APPROPRIATE ACTION WAS NOT TAKEN. THIS
978 C     IS A FATAL ERROR.
979       MSG = 'DASSL--  THE LAST STEP TERMINATED WITH A NEGATIVE'
980       CALL XERRWV(MSG,49,201,0,0,0,0,0,0.0D0,0.0D0)
981       MSG = 'DASSL--  VALUE (=I1) OF IDID AND NO APPROPRIATE'
982       CALL XERRWV(MSG,47,202,0,1,IDID,0,0,0.0D0,0.0D0)
983       MSG = 'DASSL--  ACTION WAS TAKEN. RUN TERMINATED'
984       CALL XERRWV(MSG,41,203,1,0,0,0,0,0.0D0,0.0D0)
985       RETURN
986 110   CONTINUE
987       IWORK(LNSTL)=IWORK(LNST)
988 C
989 C-----------------------------------------------------------------------
990 C     THIS BLOCK IS EXECUTED ON ALL CALLS.
991 C     THE ERROR TOLERANCE PARAMETERS ARE
992 C     CHECKED, AND THE WORK ARRAY POINTERS
993 C     ARE SET.
994 C-----------------------------------------------------------------------
995 C
996 200   CONTINUE
997 C     CHECK RTOL,ATOL
998       NZFLG=0
999       RTOLI=RTOL(1)
1000       ATOLI=ATOL(1)
1001       DO 210 I=1,NEQ
1002          IF(INFO(2).EQ.1)RTOLI=RTOL(I)
1003          IF(INFO(2).EQ.1)ATOLI=ATOL(I)
1004          IF(RTOLI.GT.0.0D0.OR.ATOLI.GT.0.0D0)NZFLG=1
1005          IF(RTOLI.LT.0.0D0)GO TO 706
1006          IF(ATOLI.LT.0.0D0)GO TO 707
1007 210      CONTINUE
1008       IF(NZFLG.EQ.0)GO TO 708
1009 C
1010 C     SET UP RWORK STORAGE.IWORK STORAGE IS FIXED
1011 C     IN DATA STATEMENT.
1012       LG0=LDELTA+NEQ
1013       LG1=LG0+NG
1014       LGX=LG1+NG
1015       LE=LGX+NG
1016       LWT=LE+NEQ
1017       LPHI=LWT+NEQ
1018       LPD=LPHI+(IWORK(LMXORD)+1)*NEQ
1019       LWM=LPD
1020       NTEMP=NPD+IWORK(LNPD)
1021       IF(INFO(1).EQ.1)GO TO 400
1022 C
1023 C-----------------------------------------------------------------------
1024 C     THIS BLOCK IS EXECUTED ON THE INITIAL CALL
1025 C     ONLY. SET THE INITIAL STEP SIZE, AND
1026 C     THE ERROR WEIGHT VECTOR, AND PHI.
1027 C     COMPUTE INITIAL YPRIME, IF NECESSARY.
1028 C-----------------------------------------------------------------------
1029 C
1030 300   CONTINUE
1031       TN=T
1032       IDID=1
1033 C
1034 C     SET ERROR WEIGHT VECTOR WT
1035       CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR)
1036 c      if(ierror.gt.0) return
1037       DO 305 I = 1,NEQ
1038          IF(RWORK(LWT+I-1).LE.0.0D0) GO TO 713
1039 305      CONTINUE
1040 C
1041 C     COMPUTE UNIT ROUNDOFF AND HMIN
1042       UROUND = DLAMCH('P')
1043       RWORK(LROUND) = UROUND
1044       HMIN = 4.0D0*UROUND*DMAX1(DABS(T),DABS(TOUT))
1045 C
1046 C     CHECK INITIAL INTERVAL TO SEE THAT IT IS LONG ENOUGH
1047       TDIST = DABS(TOUT - T)
1048       IF(TDIST .LT. HMIN) GO TO 714
1049 C
1050 C     CHECK H0, IF THIS WAS INPUT
1051       IF (INFO(8) .EQ. 0) GO TO 310
1052          HO = RWORK(LH)
1053          IF ((TOUT - T)*HO .LT. 0.0D0) GO TO 711
1054          IF (HO .EQ. 0.0D0) GO TO 712
1055          GO TO 320
1056 310    CONTINUE
1057 C
1058 C     COMPUTE INITIAL STEPSIZE, TO BE USED BY EITHER
1059 C     DDASTP OR DDAINI, DEPENDING ON INFO(11)
1060       HO = 0.001D0*TDIST
1061       YPNORM = DDANRM(NEQ,YPRIME,RWORK(LWT),RPAR,IPAR)
1062       IF (YPNORM .GT. 0.5D0/HO) HO = 0.5D0/YPNORM
1063       HO = DSIGN(HO,TOUT-T)
1064 C     ADJUST HO IF NECESSARY TO MEET HMAX BOUND
1065 320   IF (INFO(7) .EQ. 0) GO TO 330
1066          RH = DABS(HO)/RWORK(LHMAX)
1067          IF (RH .GT. 1.0D0) HO = HO/RH
1068 C     COMPUTE TSTOP, IF APPLICABLE
1069 330   IF (INFO(4) .EQ. 0) GO TO 340
1070          TSTOP = RWORK(LTSTOP)
1071          IF ((TSTOP - T)*HO .LT. 0.0D0) GO TO 715
1072          IF ((T + HO - TSTOP)*HO .GT. 0.0D0) HO = TSTOP - T
1073          IF ((TSTOP - TOUT)*HO .LT. 0.0D0) GO TO 709
1074 C
1075 C     COMPUTE INITIAL DERIVATIVE, UPDATING TN AND Y, IF APPLICABLE
1076 340   IF (INFO(11) .EQ. 0) GO TO 350
1077       CALL DDAINI(TN,Y,YPRIME,NEQ,
1078      *  RES,JAC,HO,RWORK(LWT),IDID,RPAR,IPAR,
1079      *  RWORK(LPHI),RWORK(LDELTA),RWORK(LE),
1080      *  RWORK(LWM),IWORK(LIWM),HMIN,RWORK(LROUND),
1081      *  INFO(10),NTEMP)
1082       if(ierror.gt.0) return
1083       IF (IDID .LT. 0) GO TO 390
1084 C
1085 C     LOAD H WITH H0.  STORE H IN RWORK(LH)
1086 350   H = HO
1087       RWORK(LH) = H
1088 C
1089 C     LOAD Y AND H*YPRIME INTO PHI(*,1) AND PHI(*,2)
1090 360   ITEMP = LPHI + NEQ
1091       DO 370 I = 1,NEQ
1092          RWORK(LPHI + I - 1) = Y(I)
1093 370      RWORK(ITEMP + I - 1) = H*YPRIME(I)
1094 C
1095 C     INITIALIZE T0 IN RWORK AND CHECK FOR A ZERO OF G NEAR THE
1096 C     INITIAL T.
1097 C
1098       RWORK(LT0) = T
1099       IWORK(LIRFND) = 0
1100       RWORK(LPSI)=H
1101       RWORK(LPSI+1)=2.0D0*H
1102       IWORK(LKOLD)=1
1103       IF(NG .EQ. 0) GO TO 390
1104       CALL DRCHEK(1,G,NG,NEQ,T,TOUT,Y,RWORK(LE),RWORK(LPHI),
1105      *  RWORK(LPSI),IWORK(LKOLD),RWORK(LG0),RWORK(LG1),
1106      *  RWORK(LGX),JROOT,IRT,RWORK(LROUND),INFO(3),
1107      *  RWORK,IWORK,RPAR,IPAR)
1108       if(ierror.gt.0) return
1109       IF(IRT .NE. 0) GO TO 732
1110 C
1111 C     Check for a root in the interval (T0,TN], unless DDASRT
1112 C     did not have to initialize YPRIME.
1113 C
1114       IF(NG .EQ. 0 .OR. INFO(11) .EQ. 0) GO TO 390
1115       CALL DRCHEK(3,G,NG,NEQ,TN,TOUT,Y,RWORK(LE),RWORK(LPHI),
1116      *  RWORK(LPSI),IWORK(LKOLD),RWORK(LG0),RWORK(LG1),
1117      *  RWORK(LGX),JROOT,IRT,RWORK(LROUND),INFO(3),
1118      *  RWORK,IWORK,RPAR,IPAR)
1119       if (ierror.gt.0) return
1120       IF(IRT .NE. 1) GO TO 390
1121       IWORK(LIRFND) = 1
1122       IDID = 4
1123       T = RWORK(LT0)
1124 c*SS* 1997 next line added to return current value of yprime
1125       call dcopy(neq,RWORK(LE),1,YPRIME,1)
1126       GO TO 580
1127 C
1128 390   GO TO 500
1129 C
1130 C-------------------------------------------------------
1131 C     THIS BLOCK IS FOR CONTINUATION CALLS ONLY. ITS
1132 C     PURPOSE IS TO CHECK STOP CONDITIONS BEFORE
1133 C     TAKING A STEP.
1134 C     ADJUST H IF NECESSARY TO MEET HMAX BOUND
1135 C-------------------------------------------------------
1136 C
1137 400   CONTINUE
1138       UROUND=RWORK(LROUND)
1139       DONE = .FALSE.
1140       TN=RWORK(LTN)
1141       H=RWORK(LH)
1142       IF(NG .EQ. 0) GO TO 405
1143 C
1144 C     Check for a zero of G near TN.
1145 C
1146       CALL DRCHEK(2,G,NG,NEQ,TN,TOUT,Y,RWORK(LE),RWORK(LPHI),
1147      *  RWORK(LPSI),IWORK(LKOLD),RWORK(LG0),RWORK(LG1),
1148      *  RWORK(LGX),JROOT,IRT,RWORK(LROUND),INFO(3),
1149      *  RWORK,IWORK,RPAR,IPAR)
1150       if(ierror.gt.0) return
1151       IF(IRT .NE. 1) GO TO 405
1152       IWORK(LIRFND) = 1
1153       IDID = 4
1154       T = RWORK(LT0)
1155 c*SS* 1997 next line added to return current value of yprime
1156       call dcopy(neq,RWORK(LE),1,YPRIME,1)
1157       DONE = .TRUE.
1158       GO TO 490
1159 C
1160 405   CONTINUE
1161       IF(INFO(7) .EQ. 0) GO TO 410
1162          RH = DABS(H)/RWORK(LHMAX)
1163          IF(RH .GT. 1.0D0) H = H/RH
1164 410   CONTINUE
1165       IF(T .EQ. TOUT) GO TO 719
1166       IF((T - TOUT)*H .GT. 0.0D0) GO TO 711
1167       IF(INFO(4) .EQ. 1) GO TO 430
1168       IF(INFO(3) .EQ. 1) GO TO 420
1169       IF((TN-TOUT)*H.LT.0.0D0)GO TO 490
1170       CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
1171      *  RWORK(LPHI),RWORK(LPSI))
1172       if(ierror.gt.0) return
1173       T=TOUT
1174       IDID = 3
1175       DONE = .TRUE.
1176       GO TO 490
1177 420   IF((TN-T)*H .LE. 0.0D0) GO TO 490
1178       IF((TN - TOUT)*H .GT. 0.0D0) GO TO 425
1179       CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD),
1180      *  RWORK(LPHI),RWORK(LPSI))
1181       if(ierror.gt.0) return
1182       T = TN
1183       IDID = 1
1184       DONE = .TRUE.
1185       GO TO 490
1186 425   CONTINUE
1187       CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
1188      *  RWORK(LPHI),RWORK(LPSI))
1189       if(ierror.gt.0) return
1190       T = TOUT
1191       IDID = 3
1192       DONE = .TRUE.
1193       GO TO 490
1194 430   IF(INFO(3) .EQ. 1) GO TO 440
1195       TSTOP=RWORK(LTSTOP)
1196       IF((TN-TSTOP)*H.GT.0.0D0) GO TO 715
1197       IF((TSTOP-TOUT)*H.LT.0.0D0)GO TO 709
1198       IF((TN-TOUT)*H.LT.0.0D0)GO TO 450
1199       CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
1200      *   RWORK(LPHI),RWORK(LPSI))
1201       if(ierror.gt.0) return
1202       T=TOUT
1203       IDID = 3
1204       DONE = .TRUE.
1205       GO TO 490
1206 440   TSTOP = RWORK(LTSTOP)
1207       IF((TN-TSTOP)*H .GT. 0.0D0) GO TO 715
1208       IF((TSTOP-TOUT)*H .LT. 0.0D0) GO TO 709
1209       IF((TN-T)*H .LE. 0.0D0) GO TO 450
1210       IF((TN - TOUT)*H .GT. 0.0D0) GO TO 445
1211       CALL DDATRP(TN,TN,Y,YPRIME,NEQ,IWORK(LKOLD),
1212      *  RWORK(LPHI),RWORK(LPSI))
1213       if(ierror.gt.0) return
1214       T = TN
1215       IDID = 1
1216       DONE = .TRUE.
1217       GO TO 490
1218 445   CONTINUE
1219       CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,IWORK(LKOLD),
1220      *  RWORK(LPHI),RWORK(LPSI))
1221       if(ierror.gt.0) return
1222       T = TOUT
1223       IDID = 3
1224       DONE = .TRUE.
1225       GO TO 490
1226 450   CONTINUE
1227 C     CHECK WHETHER WE ARE WITH IN ROUNDOFF OF TSTOP
1228       IF(DABS(TN-TSTOP).GT.100.0D0*UROUND*
1229      *   (DABS(TN)+DABS(H)))GO TO 460
1230       CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,IWORK(LKOLD),
1231      *  RWORK(LPHI),RWORK(LPSI))
1232       if(ierror.gt.0) return
1233       IDID=2
1234       T=TSTOP
1235       DONE = .TRUE.
1236       GO TO 490
1237 460   TNEXT=TN+H
1238       IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 490
1239       H=TSTOP-TN
1240       RWORK(LH)=H
1241 C
1242 490   IF (DONE) GO TO 590
1243 C
1244 C-------------------------------------------------------
1245 C     THE NEXT BLOCK CONTAINS THE CALL TO THE
1246 C     ONE-STEP INTEGRATOR DDASTP.
1247 C     THIS IS A LOOPING POINT FOR THE INTEGRATION STEPS.
1248 C     CHECK FOR TOO MANY STEPS.
1249 C     UPDATE WT.
1250 C     CHECK FOR TOO MUCH ACCURACY REQUESTED.
1251 C     COMPUTE MINIMUM STEPSIZE.
1252 C-------------------------------------------------------
1253 C
1254 500   CONTINUE
1255 C     CHECK FOR FAILURE TO COMPUTE INITIAL YPRIME
1256       IF (IDID .EQ. -12) GO TO 527
1257 C
1258 C     CHECK FOR TOO MANY STEPS
1259       IF((IWORK(LNST)-IWORK(LNSTL)).LT.500)
1260      *   GO TO 510
1261            IDID=-1
1262            GO TO 527
1263 C
1264 C     UPDATE WT
1265 510   CALL DDAWTS(NEQ,INFO(2),RTOL,ATOL,RWORK(LPHI),
1266      *  RWORK(LWT),RPAR,IPAR)
1267       if(ierror.gt.0) return
1268       DO 520 I=1,NEQ
1269          IF(RWORK(I+LWT-1).GT.0.0D0)GO TO 520
1270            IDID=-3
1271            GO TO 527
1272 520   CONTINUE
1273 C
1274 C     TEST FOR TOO MUCH ACCURACY REQUESTED.
1275       R=DDANRM(NEQ,RWORK(LPHI),RWORK(LWT),RPAR,IPAR)*
1276      *   100.0D0*UROUND
1277       IF(R.LE.1.0D0)GO TO 525
1278 C     MULTIPLY RTOL AND ATOL BY R AND RETURN
1279       IF(INFO(2).EQ.1)GO TO 523
1280            RTOL(1)=R*RTOL(1)
1281            ATOL(1)=R*ATOL(1)
1282            IDID=-2
1283            GO TO 527
1284 523   DO 524 I=1,NEQ
1285            RTOL(I)=R*RTOL(I)
1286 524        ATOL(I)=R*ATOL(I)
1287       IDID=-2
1288       GO TO 527
1289 525   CONTINUE
1290 C
1291 C     COMPUTE MINIMUM STEPSIZE
1292       HMIN=4.0D0*UROUND*DMAX1(DABS(TN),DABS(TOUT))
1293 C
1294 C     TEST H VS. HMAX
1295       IF (INFO(7) .EQ. 0) GO TO 526
1296          RH = ABS(H)/RWORK(LHMAX)
1297          IF (RH .GT. 1.0D0) H = H/RH
1298 526   CONTINUE     
1299 C
1300       CALL DDASTP(TN,Y,YPRIME,NEQ,
1301      *   RES,JAC,H,RWORK(LWT),INFO(1),IDID,RPAR,IPAR,
1302      *   RWORK(LPHI),RWORK(LDELTA),RWORK(LE),
1303      *   RWORK(LWM),IWORK(LIWM),
1304      *   RWORK(LALPHA),RWORK(LBETA),RWORK(LGAMMA),
1305      *   RWORK(LPSI),RWORK(LSIGMA),
1306      *   RWORK(LCJ),RWORK(LCJOLD),RWORK(LHOLD),
1307      *   RWORK(LS),HMIN,RWORK(LROUND),
1308      *   IWORK(LPHASE),IWORK(LJCALC),IWORK(LK),
1309      *   IWORK(LKOLD),IWORK(LNS),INFO(10),NTEMP)
1310       if(ierror.gt.0) return
1311 527   IF(IDID.LT.0)GO TO 600
1312 C
1313 C--------------------------------------------------------
1314 C     THIS BLOCK HANDLES THE CASE OF A SUCCESSFUL RETURN
1315 C     FROM DDASTP (IDID=1).  TEST FOR STOP CONDITIONS.
1316 C--------------------------------------------------------
1317 C
1318       IF(NG .EQ. 0) GO TO 529
1319 C
1320 C     Check for a zero of G near TN.
1321 C
1322       CALL DRCHEK(3,G,NG,NEQ,TN,TOUT,Y,RWORK(LE),RWORK(LPHI),
1323      *  RWORK(LPSI),IWORK(LKOLD),RWORK(LG0),RWORK(LG1),
1324      *  RWORK(LGX),JROOT,IRT,RWORK(LROUND),INFO(3),
1325      *  RWORK,IWORK,RPAR,IPAR)
1326       if(ierror.gt.0) return
1327       IF(IRT .NE. 1) GO TO 529
1328       IWORK(LIRFND) = 1
1329       IDID = 4
1330       T = RWORK(LT0)
1331 c*SS* 1997 next line added to return current value of yprime
1332       call dcopy(neq,RWORK(LE),1,YPRIME,1)
1333       GO TO 580
1334 C
1335 529   CONTINUE
1336       IF(INFO(4).NE.0)GO TO 540
1337            IF(INFO(3).NE.0)GO TO 530
1338              IF((TN-TOUT)*H.LT.0.0D0)GO TO 500
1339              CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
1340      *         IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1341       if(ierror.gt.0) return
1342              IDID=3
1343              T=TOUT
1344              GO TO 580
1345 530          IF((TN-TOUT)*H.GE.0.0D0)GO TO 535
1346              T=TN
1347              IDID=1
1348              GO TO 580
1349 535          CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
1350      *         IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1351       if(ierror.gt.0) return
1352              IDID=3
1353              T=TOUT
1354              GO TO 580
1355 540   IF(INFO(3).NE.0)GO TO 550
1356       IF((TN-TOUT)*H.LT.0.0D0)GO TO 542
1357          CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
1358      *     IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1359       if(ierror.gt.0) return
1360          T=TOUT
1361          IDID=3
1362          GO TO 580
1363 542   IF(DABS(TN-TSTOP).LE.100.0D0*UROUND*
1364      *   (DABS(TN)+DABS(H)))GO TO 545
1365       TNEXT=TN+H
1366       IF((TNEXT-TSTOP)*H.LE.0.0D0)GO TO 500
1367       H=TSTOP-TN
1368       GO TO 500
1369 545   CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,
1370      *  IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1371       if(ierror.gt.0) return
1372       IDID=2
1373       T=TSTOP
1374       GO TO 580
1375 550   IF((TN-TOUT)*H.GE.0.0D0)GO TO 555
1376       IF(DABS(TN-TSTOP).LE.100.0D0*UROUND*(DABS(TN)+DABS(H)))GO TO 552
1377       T=TN
1378       IDID=1
1379       GO TO 580
1380 552   CALL DDATRP(TN,TSTOP,Y,YPRIME,NEQ,
1381      *  IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1382       if(ierror.gt.0) return
1383       IDID=2
1384       T=TSTOP
1385       GO TO 580
1386 555   CALL DDATRP(TN,TOUT,Y,YPRIME,NEQ,
1387      *   IWORK(LKOLD),RWORK(LPHI),RWORK(LPSI))
1388       if(ierror.gt.0) return
1389       T=TOUT
1390       IDID=3
1391 580   CONTINUE
1392 C
1393 C--------------------------------------------------------
1394 C     ALL SUCCESSFUL RETURNS FROM DDASRT ARE MADE FROM
1395 C     THIS BLOCK.
1396 C--------------------------------------------------------
1397 C
1398 590   CONTINUE
1399       RWORK(LTN)=TN
1400       RWORK(LH)=H
1401       RWORK(LTLAST) = T
1402       RETURN
1403 C
1404 C-----------------------------------------------------------------------
1405 C     THIS BLOCK HANDLES ALL UNSUCCESSFUL
1406 C     RETURNS OTHER THAN FOR ILLEGAL INPUT.
1407 C-----------------------------------------------------------------------
1408 C
1409 600   CONTINUE
1410       ITEMP=-IDID
1411       GO TO (610,620,630,690,690,640,650,660,670,675,
1412      *  680,685), ITEMP
1413 C
1414 C     THE MAXIMUM NUMBER OF STEPS WAS TAKEN BEFORE
1415 C     REACHING TOUT
1416 610   MSG = 'DASSL--  AT CURRENT T (=R1)  500 STEPS'
1417       CALL XERRWV(MSG,38,610,0,0,0,0,1,TN,0.0D0)
1418       MSG = 'DASSL--  TAKEN ON THIS CALL BEFORE REACHING TOUT'
1419       CALL XERRWV(MSG,48,611,0,0,0,0,0,0.0D0,0.0D0)
1420       GO TO 690
1421 C
1422 C     TOO MUCH ACCURACY FOR MACHINE PRECISION
1423 620   MSG = 'DASSL--  AT T (=R1) TOO MUCH ACCURACY REQUESTED'
1424       CALL XERRWV(MSG,47,620,0,0,0,0,1,TN,0.0D0)
1425       MSG = 'DASSL--  FOR PRECISION OF MACHINE. RTOL AND ATOL'
1426       CALL XERRWV(MSG,48,621,0,0,0,0,0,0.0D0,0.0D0)
1427       MSG = 'DASSL--  WERE INCREASED TO APPROPRIATE VALUES'
1428       CALL XERRWV(MSG,45,622,0,0,0,0,0,0.0D0,0.0D0)
1429 C
1430       GO TO 690
1431 C     WT(I) .LE. 0.0D0 FOR SOME I (NOT AT START OF PROBLEM)
1432 630   MSG = 'DASSL--  AT T (=R1) SOME ELEMENT OF WT'
1433       CALL XERRWV(MSG,38,630,0,0,0,0,1,TN,0.0D0)
1434       MSG = 'DASSL--  HAS BECOME .LE. 0.0'
1435       CALL XERRWV(MSG,28,631,0,0,0,0,0,0.0D0,0.0D0)
1436       GO TO 690
1437 C
1438 C     ERROR TEST FAILED REPEATEDLY OR WITH H=HMIN
1439 640   MSG = 'DASSL--  AT T (=R1) AND STEPSIZE H (=R2) THE'
1440       CALL XERRWV(MSG,44,640,0,0,0,0,2,TN,H)
1441       MSG='DASSL--  ERROR TEST FAILED REPEATEDLY OR WITH ABS(H)=HMIN'
1442       CALL XERRWV(MSG,57,641,0,0,0,0,0,0.0D0,0.0D0)
1443       GO TO 690
1444 C
1445 C     CORRECTOR CONVERGENCE FAILED REPEATEDLY OR WITH H=HMIN
1446 650   MSG = 'DASSL--  AT T (=R1) AND STEPSIZE H (=R2) THE'
1447       CALL XERRWV(MSG,44,650,0,0,0,0,2,TN,H)
1448       MSG = 'DASSL--  CORRECTOR FAILED TO CONVERGE REPEATEDLY'
1449       CALL XERRWV(MSG,48,651,0,0,0,0,0,0.0D0,0.0D0)
1450       MSG = 'DASSL--  OR WITH ABS(H)=HMIN'
1451       CALL XERRWV(MSG,28,652,0,0,0,0,0,0.0D0,0.0D0)
1452       GO TO 690
1453 C
1454 C     THE ITERATION MATRIX IS SINGULAR
1455 660   MSG = 'DASSL--  AT T (=R1) AND STEPSIZE H (=R2) THE'
1456       CALL XERRWV(MSG,44,660,0,0,0,0,2,TN,H)
1457       MSG = 'DASSL--  ITERATION MATRIX IS SINGULAR'
1458       CALL XERRWV(MSG,37,661,0,0,0,0,0,0.0D0,0.0D0)
1459       GO TO 690
1460 C
1461 C     CORRECTOR FAILURE PRECEEDED BY ERROR TEST FAILURES.
1462 670   MSG = 'DASSL--  AT T (=R1) AND STEPSIZE H (=R2) THE'
1463       CALL XERRWV(MSG,44,670,0,0,0,0,2,TN,H)
1464       MSG = 'DASSL--  CORRECTOR COULD NOT CONVERGE.  ALSO, THE'
1465       CALL XERRWV(MSG,49,671,0,0,0,0,0,0.0D0,0.0D0)
1466       MSG = 'DASSL--  ERROR TEST FAILED REPEATEDLY.'
1467       CALL XERRWV(MSG,38,672,0,0,0,0,0,0.0D0,0.0D0)
1468       GO TO 690
1469 C
1470 C     CORRECTOR FAILURE BECAUSE IRES = -1
1471 675   MSG = 'DASSL--  AT T (=R1) AND STEPSIZE H (=R2) THE'
1472       CALL XERRWV(MSG,44,675,0,0,0,0,2,TN,H)
1473       MSG = 'DASSL--  CORRECTOR COULD NOT CONVERGE BECAUSE'
1474       CALL XERRWV(MSG,45,676,0,0,0,0,0,0.0D0,0.0D0)
1475       MSG = 'DASSL--  IRES WAS EQUAL TO MINUS ONE'
1476       CALL XERRWV(MSG,36,677,0,0,0,0,0,0.0D0,0.0D0)
1477       GO TO 690
1478 C
1479 C     FAILURE BECAUSE IRES = -2
1480 680   MSG = 'DASSL--  AT T (=R1) AND STEPSIZE H (=R2)'
1481       CALL XERRWV(MSG,40,680,0,0,0,0,2,TN,H)
1482       MSG = 'DASSL--  IRES WAS EQUAL TO MINUS TWO'
1483       CALL XERRWV(MSG,36,681,0,0,0,0,0,0.0D0,0.0D0)
1484       GO TO 690
1485 C
1486 C     FAILED TO COMPUTE INITIAL YPRIME
1487 685   MSG = 'DASSL--  AT T (=R1) AND STEPSIZE H (=R2) THE'
1488       CALL XERRWV(MSG,44,685,0,0,0,0,2,TN,HO)
1489       MSG = 'DASSL--  INITIAL YPRIME COULD NOT BE COMPUTED'
1490       CALL XERRWV(MSG,45,686,0,0,0,0,0,0.0D0,0.0D0)
1491       GO TO 690
1492 690   CONTINUE
1493       INFO(1)=-1
1494       T=TN
1495       RWORK(LTN)=TN
1496       RWORK(LH)=H
1497       RETURN
1498 C-----------------------------------------------------------------------
1499 C     THIS BLOCK HANDLES ALL ERROR RETURNS DUE
1500 C     TO ILLEGAL INPUT, AS DETECTED BEFORE CALLING
1501 C     DDASTP. FIRST THE ERROR MESSAGE ROUTINE IS
1502 C     CALLED. IF THIS HAPPENS TWICE IN
1503 C     SUCCESSION, EXECUTION IS TERMINATED
1504 C
1505 C-----------------------------------------------------------------------
1506 701   MSG = 'DASSL--  SOME ELEMENT OF INFO VECTOR IS NOT ZERO OR ONE'
1507       CALL XERRWV(MSG,55,1,0,0,0,0,0,0.0D0,0.0D0)
1508       GO TO 750
1509 702   MSG = 'DASSL--  NEQ (=I1) .LE. 0'
1510       CALL XERRWV(MSG,25,2,0,1,NEQ,0,0,0.0D0,0.0D0)
1511       GO TO 750
1512 703   MSG = 'DASSL--  MAXORD (=I1) NOT IN RANGE'
1513       CALL XERRWV(MSG,34,3,0,1,MXORD,0,0,0.0D0,0.0D0)
1514       GO TO 750
1515 704   MSG='DASSL--  RWORK LENGTH NEEDED, LENRW (=I1), EXCEEDS LRW (=I2)'
1516       CALL XERRWV(MSG,60,4,0,2,LENRW,LRW,0,0.0D0,0.0D0)
1517       GO TO 750
1518 705   MSG='DASSL--  IWORK LENGTH NEEDED, LENIW (=I1), EXCEEDS LIW (=I2)'
1519       CALL XERRWV(MSG,60,5,0,2,LENIW,LIW,0,0.0D0,0.0D0)
1520       GO TO 750
1521 706   MSG = 'DASSL--  SOME ELEMENT OF RTOL IS .LT. 0'
1522       CALL XERRWV(MSG,39,6,0,0,0,0,0,0.0D0,0.0D0)
1523       GO TO 750
1524 707   MSG = 'DASSL--  SOME ELEMENT OF ATOL IS .LT. 0'
1525       CALL XERRWV(MSG,39,7,0,0,0,0,0,0.0D0,0.0D0)
1526       GO TO 750
1527 708   MSG = 'DASSL--  ALL ELEMENTS OF RTOL AND ATOL ARE ZERO'
1528       CALL XERRWV(MSG,47,8,0,0,0,0,0,0.0D0,0.0D0)
1529       GO TO 750
1530 709   MSG='DASSL--  INFO(4) = 1 AND TSTOP (=R1) BEHIND TOUT (=R2)'
1531       CALL XERRWV(MSG,54,9,0,0,0,0,2,TSTOP,TOUT)
1532       GO TO 750
1533 710   MSG = 'DASSL--  HMAX (=R1) .LT. 0.0'
1534       CALL XERRWV(MSG,28,10,0,0,0,0,1,HMAX,0.0D0)
1535       GO TO 750
1536 711   MSG = 'DASSL--  TOUT (=R1) BEHIND T (=R2)'
1537       CALL XERRWV(MSG,34,11,0,0,0,0,2,TOUT,T)
1538       GO TO 750
1539 712   MSG = 'DASSL--  INFO(8)=1 AND H0=0.0'
1540       CALL XERRWV(MSG,29,12,0,0,0,0,0,0.0D0,0.0D0)
1541       GO TO 750
1542 713   MSG = 'DASSL--  SOME ELEMENT OF WT IS .LE. 0.0'
1543       CALL XERRWV(MSG,39,13,0,0,0,0,0,0.0D0,0.0D0)
1544       GO TO 750
1545 714   MSG='DASSL-- TOUT (=R1) TOO CLOSE TO T (=R2) TO START INTEGRATION'
1546       CALL XERRWV(MSG,60,14,0,0,0,0,2,TOUT,T)
1547       GO TO 750
1548 715   MSG = 'DASSL--  INFO(4)=1 AND TSTOP (=R1) BEHIND T (=R2)'
1549       CALL XERRWV(MSG,49,15,0,0,0,0,2,TSTOP,T)
1550       GO TO 750
1551 717   MSG = 'DASSL--  ML (=I1) ILLEGAL. EITHER .LT. 0 OR .GT. NEQ'
1552       CALL XERRWV(MSG,52,17,0,1,IWORK(LML),0,0,0.0D0,0.0D0)
1553       GO TO 750
1554 718   MSG = 'DASSL--  MU (=I1) ILLEGAL. EITHER .LT. 0 OR .GT. NEQ'
1555       CALL XERRWV(MSG,52,18,0,1,IWORK(LMU),0,0,0.0D0,0.0D0)
1556       GO TO 750
1557 719   MSG = 'DASSL--  TOUT (=R1) IS EQUAL TO T (=R2)'
1558       CALL XERRWV(MSG,39,19,0,0,0,0,2,TOUT,T)
1559       GO TO 750
1560 730   MSG = 'DASSL--  NG (=I1) .LT. 0'
1561       CALL XERRWV(MSG,24,30,1,1,NG,0,0,0.0D0,0.0D0)
1562       GO TO 750
1563 732   MSG = 'DASSL--  ONE OR MORE COMPONENTS OF G HAS A ROOT'
1564       CALL XERRWV(MSG,47,32,1,0,0,0,0,0.0D0,0.0D0)
1565       MSG = '         TOO NEAR TO THE INITIAL POINT'
1566       CALL XERRWV(MSG,38,32,1,0,0,0,0,0.0D0,0.0D0)
1567 750   IF(INFO(1).EQ.-1) GO TO 760
1568       INFO(1)=-1
1569       IDID=-33
1570       RETURN
1571 760   MSG = 'DASSL--  REPEATED OCCURRENCES OF ILLEGAL INPUT'
1572       CALL XERRWV(MSG,46,801,0,0,0,0,0,0.0D0,0.0D0)
1573 770   MSG = 'DASSL--  RUN TERMINATED. APPARENT INFINITE LOOP'
1574       CALL XERRWV(MSG,47,802,1,0,0,0,0,0.0D0,0.0D0)
1575       RETURN
1576 C-----------END OF SUBROUTINE DDASRT------------------------------------
1577       END
1578       SUBROUTINE DRCHEK (JOB, G, NG, NEQ, TN, TOUT, Y, YP, PHI, PSI,
1579      *  KOLD, G0, G1, GX, JROOT, IRT, UROUND, INFO3, RWORK, IWORK,
1580      *  RPAR, IPAR)
1581 C
1582 C***BEGIN PROLOGUE  DRCHEK
1583 C***REFER TO DDASRT
1584 C***ROUTINES CALLED  DDATRP, DROOTS, DCOPY
1585 C***DATE WRITTEN   821001   (YYMMDD)
1586 C***REVISION DATE  900926   (YYMMDD)
1587 C***END PROLOGUE  DRCHEK
1588 C
1589       
1590       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1591
1592       include 'stack.h'
1593
1594       PARAMETER (LNGE=16, LIRFND=18, LLAST=19, LIMAX=20,
1595      *           LT0=41, LTLAST=42, LALPHR=43, LX2=44)
1596       EXTERNAL G
1597       INTEGER JOB, NG, NEQ, KOLD, JROOT, IRT, INFO3, IWORK, IPAR
1598       DOUBLE PRECISION TN, TOUT, Y, YP, PHI, PSI, G0, G1, GX, UROUND,
1599      *  RWORK, RPAR
1600       DIMENSION  Y(*), YP(*), PHI(NEQ,*), PSI(*),
1601      1  G0(*), G1(*), GX(*), JROOT(*), RWORK(*), IWORK(*)
1602       INTEGER I, JFLAG
1603       DOUBLE PRECISION H
1604       DOUBLE PRECISION HMING, T1, TEMP1, TEMP2, X
1605       LOGICAL ZROOT
1606 C-----------------------------------------------------------------------
1607 C THIS ROUTINE CHECKS FOR THE PRESENCE OF A ROOT IN THE
1608 C VICINITY OF THE CURRENT T, IN A MANNER DEPENDING ON THE
1609 C INPUT FLAG JOB.  IT CALLS SUBROUTINE DROOTS TO LOCATE THE ROOT
1610 C AS PRECISELY AS POSSIBLE.
1611 C
1612 C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, DRCHEK
1613 C USES THE FOLLOWING FOR COMMUNICATION..
1614 C JOB    = INTEGER FLAG INDICATING TYPE OF CALL..
1615 C          JOB = 1 MEANS THE PROBLEM IS BEING INITIALIZED, AND DRCHEK
1616 C                  IS TO LOOK FOR A ROOT AT OR VERY NEAR THE INITIAL T.
1617 C          JOB = 2 MEANS A CONTINUATION CALL TO THE SOLVER WAS JUST
1618 C                  MADE, AND DRCHEK IS TO CHECK FOR A ROOT IN THE
1619 C                  RELEVANT PART OF THE STEP LAST TAKEN.
1620 C          JOB = 3 MEANS A SUCCESSFUL STEP WAS JUST TAKEN, AND DRCHEK
1621 C                  IS TO LOOK FOR A ROOT IN THE INTERVAL OF THE STEP.
1622 C G0     = ARRAY OF LENGTH NG, CONTAINING THE VALUE OF G AT T = T0.
1623 C          G0 IS INPUT FOR JOB .GE. 2 AND ON OUTPUT IN ALL CASES.
1624 C G1,GX  = ARRAYS OF LENGTH NG FOR WORK SPACE.
1625 C IRT    = COMPLETION FLAG..
1626 C          IRT = 0  MEANS NO ROOT WAS FOUND.
1627 C          IRT = -1 MEANS JOB = 1 AND A ROOT WAS FOUND TOO NEAR TO T.
1628 C          IRT = 1  MEANS A LEGITIMATE ROOT WAS FOUND (JOB = 2 OR 3).
1629 C                   ON RETURN, T0 IS THE ROOT LOCATION, AND Y IS THE
1630 C                   CORRESPONDING SOLUTION VECTOR.
1631 C T0     = VALUE OF T AT ONE ENDPOINT OF INTERVAL OF INTEREST.  ONLY
1632 C          ROOTS BEYOND T0 IN THE DIRECTION OF INTEGRATION ARE SOUGHT.
1633 C          T0 IS INPUT IF JOB .GE. 2, AND OUTPUT IN ALL CASES.
1634 C          T0 IS UPDATED BY DRCHEK, WHETHER A ROOT IS FOUND OR NOT.
1635 C          STORED IN THE GLOBAL ARRAY RWORK.
1636 C TLAST  = LAST VALUE OF T RETURNED BY THE SOLVER (INPUT ONLY).
1637 C          STORED IN THE GLOBAL ARRAY RWORK.
1638 C TOUT   = FINAL OUTPUT TIME FOR THE SOLVER.
1639 C IRFND  = INPUT FLAG SHOWING WHETHER THE LAST STEP TAKEN HAD A ROOT.
1640 C          IRFND = 1 IF IT DID, = 0 IF NOT.
1641 C          STORED IN THE GLOBAL ARRAY IWORK.
1642 C INFO3  = COPY OF INFO(3) (INPUT ONLY).
1643 C-----------------------------------------------------------------------
1644 C     
1645       H = PSI(1)
1646       IRT = 0
1647       DO 10 I = 1,NG
1648  10     JROOT(I) = 0
1649       HMING = (DABS(TN) + DABS(H))*UROUND*100.0D0
1650 C
1651       GO TO (100, 200, 300), JOB
1652 C
1653 C EVALUATE G AT INITIAL T (STORED IN RWORK(LT0)), AND CHECK FOR
1654 C ZERO VALUES.----------------------------------------------------------
1655  100  CONTINUE
1656       CALL DDATRP(TN,RWORK(LT0),Y,YP,NEQ,KOLD,PHI,PSI)
1657       if(ierror.gt.0) return
1658       CALL G (NEQ, RWORK(LT0), Y, NG, G0, RPAR, IPAR)
1659       if(ierror.gt.0) return
1660       IWORK(LNGE) = 1
1661       ZROOT = .FALSE.
1662       DO 110 I = 1,NG
1663  110    IF (DABS(G0(I)) .LE. 0.0D0) ZROOT = .TRUE.
1664       IF (.NOT. ZROOT) GO TO 190
1665 C G HAS A ZERO AT T.  LOOK AT G AT T + (SMALL INCREMENT). --------------
1666       TEMP1 = DSIGN(HMING,H)
1667       RWORK(LT0) = RWORK(LT0) + TEMP1
1668       TEMP2 = TEMP1/H
1669       DO 120 I = 1,NEQ
1670  120    Y(I) = Y(I) + TEMP2*PHI(I,2)
1671       CALL G (NEQ, RWORK(LT0), Y, NG, G0, RPAR, IPAR)
1672       if(ierror.gt.0) return
1673       IWORK(LNGE) = IWORK(LNGE) + 1
1674       ZROOT = .FALSE.
1675       DO 130 I = 1,NG
1676  130    IF (DABS(G0(I)) .LE. 0.0D0) ZROOT = .TRUE.
1677       IF (.NOT. ZROOT) GO TO 190
1678 C G HAS A ZERO AT T AND ALSO CLOSE TO T.  TAKE ERROR RETURN. -----------
1679       IRT = -1
1680       RETURN
1681 C
1682  190  CONTINUE
1683       RETURN
1684 C
1685 C
1686  200  CONTINUE
1687       IF (IWORK(LIRFND) .EQ. 0) GO TO 260
1688 C IF A ROOT WAS FOUND ON THE PREVIOUS STEP, EVALUATE G0 = G(T0). -------
1689       CALL DDATRP (TN, RWORK(LT0), Y, YP, NEQ, KOLD, PHI, PSI)
1690       if(ierror.gt.0) return
1691       CALL G (NEQ, RWORK(LT0), Y, NG, G0, RPAR, IPAR)
1692       if(ierror.gt.0) return
1693       IWORK(LNGE) = IWORK(LNGE) + 1
1694       ZROOT = .FALSE.
1695       DO 210 I = 1,NG
1696  210    IF (DABS(G0(I)) .LE. 0.0D0) ZROOT = .TRUE.
1697       IF (.NOT. ZROOT) GO TO 260
1698 C G HAS A ZERO AT T0.  LOOK AT G AT T + (SMALL INCREMENT). -------------
1699       TEMP1 = DSIGN(HMING,H)
1700       RWORK(LT0) = RWORK(LT0) + TEMP1
1701       IF ((RWORK(LT0) - TN)*H .LT. 0.0D0) GO TO 230
1702       TEMP2 = TEMP1/H
1703       DO 220 I = 1,NEQ
1704  220    Y(I) = Y(I) + TEMP2*PHI(I,2)
1705       GO TO 240
1706  230  CALL DDATRP (TN, RWORK(LT0), Y, YP, NEQ, KOLD, PHI, PSI)
1707       if(ierror.gt.0) return
1708  240  CALL G (NEQ, RWORK(LT0), Y, NG, G0, RPAR, IPAR)
1709       if(ierror.gt.0) return
1710       IWORK(LNGE) = IWORK(LNGE) + 1
1711       ZROOT = .FALSE.
1712       DO 250 I = 1,NG
1713         IF (DABS(G0(I)) .GT. 0.0D0) GO TO 250
1714         JROOT(I) = 1
1715         ZROOT = .TRUE.
1716  250    CONTINUE
1717       IF (.NOT. ZROOT) GO TO 260
1718 C G HAS A ZERO AT T0 AND ALSO CLOSE TO T0.  RETURN ROOT. ---------------
1719       IRT = 1
1720       RETURN
1721 C     HERE, G0 DOES NOT HAVE A ROOT
1722 C G0 HAS NO ZERO COMPONENTS.  PROCEED TO CHECK RELEVANT INTERVAL. ------
1723  260  IF (TN .EQ. RWORK(LTLAST)) GO TO 390
1724 C
1725  300  CONTINUE
1726 C SET T1 TO TN OR TOUT, WHICHEVER COMES FIRST, AND GET G AT T1. --------
1727       IF (INFO3 .EQ. 1) GO TO 310
1728       IF ((TOUT - TN)*H .GE. 0.0D0) GO TO 310
1729       T1 = TOUT
1730       IF ((T1 - RWORK(LT0))*H .LE. 0.0D0) GO TO 390
1731       CALL DDATRP (TN, T1, Y, YP, NEQ, KOLD, PHI, PSI)
1732       if(ierror.gt.0) return
1733       GO TO 330
1734  310  T1 = TN
1735       DO 320 I = 1,NEQ
1736  320    Y(I) = PHI(I,1)
1737  330  CALL G (NEQ, T1, Y, NG, G1, RPAR, IPAR)
1738       if(ierror.gt.0) return
1739       IWORK(LNGE) = IWORK(LNGE) + 1
1740 C CALL DROOTS TO SEARCH FOR ROOT IN INTERVAL FROM T0 TO T1. ------------
1741       JFLAG = 0
1742  350  CONTINUE
1743       CALL DROOTS (NG, HMING, JFLAG, RWORK(LT0), T1, G0, G1, GX, X,
1744      *             JROOT, IWORK(LIMAX), IWORK(LLAST), RWORK(LALPHR),
1745      *             RWORK(LX2))
1746       if(ierror.gt.0) return
1747       IF (JFLAG .GT. 1) GO TO 360
1748       CALL DDATRP (TN, X, Y, YP, NEQ, KOLD, PHI, PSI)
1749       if(ierror.gt.0) return
1750       CALL G (NEQ, X, Y, NG, GX, RPAR, IPAR)
1751       if(ierror.gt.0) return
1752       IWORK(LNGE) = IWORK(LNGE) + 1
1753       GO TO 350
1754  360  RWORK(LT0) = X
1755       CALL DCOPY (NG, GX, 1, G0, 1)
1756       IF (JFLAG .EQ. 4) GO TO 390
1757 C FOUND A ROOT.  INTERPOLATE TO X AND RETURN. --------------------------
1758       CALL DDATRP (TN, X, Y, YP, NEQ, KOLD, PHI, PSI)
1759       if(ierror.gt.0) return
1760       IRT = 1
1761       RETURN
1762 C
1763  390  CONTINUE
1764       RETURN
1765 C---------------------- END OF SUBROUTINE DRCHEK -----------------------
1766       END
1767       SUBROUTINE DROOTS (NG, HMIN, JFLAG, X0, X1, G0, G1, GX, X, JROOT,
1768      *                   IMAX, LAST, ALPHA, X2)
1769 c      subroutine roots (ng, hmin, jflag, x0, x1, g0, g1, gx, x, jroot)
1770 C
1771 C***BEGIN PROLOGUE  DROOTS
1772 C***REFER TO DDASRT
1773 C***ROUTINES CALLED  DCOPY
1774 C***DATE WRITTEN   821001   (YYMMDD)
1775 C***REVISION DATE  900926   (YYMMDD)
1776 C***END PROLOGUE  DROOTS
1777 C
1778       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1779
1780       include 'stack.h'
1781       
1782       INTEGER NG, JFLAG, JROOT, IMAX, LAST
1783       DOUBLE PRECISION HMIN, X0, X1, G0, G1, GX, X, ALPHA, X2
1784       DIMENSION G0(NG), G1(NG), GX(NG), JROOT(NG)
1785 C-----------------------------------------------------------------------
1786 C THIS SUBROUTINE FINDS THE LEFTMOST ROOT OF A SET OF ARBITRARY
1787 C FUNCTIONS GI(X) (I = 1,...,NG) IN AN INTERVAL (X0,X1).  ONLY ROOTS
1788 C OF ODD MULTIPLICITY (I.E. CHANGES OF SIGN OF THE GI) ARE FOUND.
1789 C HERE THE SIGN OF X1 - X0 IS ARBITRARY, BUT IS CONSTANT FOR A GIVEN
1790 C PROBLEM, AND -LEFTMOST- MEANS NEAREST TO X0.
1791 C THE VALUES OF THE VECTOR-VALUED FUNCTION G(X) = (GI, I=1...NG)
1792 C ARE COMMUNICATED THROUGH THE CALL SEQUENCE OF DROOTS.
1793 C THE METHOD USED IS THE ILLINOIS ALGORITHM.
1794 C
1795 C REFERENCE..
1796 C KATHIE L. HIEBERT AND LAWRENCE F. SHAMPINE, IMPLICITLY DEFINED
1797 C OUTPUT POINTS FOR SOLUTIONS OF ODE-S, SANDIA REPORT SAND80-0180,
1798 C FEBRUARY, 1980.
1799 C
1800 C DESCRIPTION OF PARAMETERS.
1801 C
1802 C NG     = NUMBER OF FUNCTIONS GI, OR THE NUMBER OF COMPONENTS OF
1803 C          THE VECTOR VALUED FUNCTION G(X).  INPUT ONLY.
1804 C
1805 C HMIN   = RESOLUTION PARAMETER IN X.  INPUT ONLY.  WHEN A ROOT IS
1806 C          FOUND, IT IS LOCATED ONLY TO WITHIN AN ERROR OF HMIN IN X.
1807 C          TYPICALLY, HMIN SHOULD BE SET TO SOMETHING ON THE ORDER OF
1808 C               100 * UROUND * MAX(ABS(X0),ABS(X1)),
1809 C          WHERE UROUND IS THE UNIT ROUNDOFF OF THE MACHINE.
1810 C
1811 C JFLAG  = INTEGER FLAG FOR INPUT AND OUTPUT COMMUNICATION.
1812 C
1813 C          ON INPUT, SET JFLAG = 0 ON THE FIRST CALL FOR THE PROBLEM,
1814 C          AND LEAVE IT UNCHANGED UNTIL THE PROBLEM IS COMPLETED.
1815 C          (THE PROBLEM IS COMPLETED WHEN JFLAG .GE. 2 ON RETURN.)
1816 C
1817 C          ON OUTPUT, JFLAG HAS THE FOLLOWING VALUES AND MEANINGS..
1818 C          JFLAG = 1 MEANS DROOTS NEEDS A VALUE OF G(X).  SET GX = G(X)
1819 C                    AND CALL DROOTS AGAIN.
1820 C          JFLAG = 2 MEANS A ROOT HAS BEEN FOUND.  THE ROOT IS
1821 C                    AT X, AND GX CONTAINS G(X).  (ACTUALLY, X IS THE
1822 C                    RIGHTMOST APPROXIMATION TO THE ROOT ON AN INTERVAL
1823 C                    (X0,X1) OF SIZE HMIN OR LESS.)
1824 C          JFLAG = 3 MEANS X = X1 IS A ROOT, WITH ONE OR MORE OF THE GI
1825 C                    BEING ZERO AT X1 AND NO SIGN CHANGES IN (X0,X1).
1826 C                    GX CONTAINS G(X) ON OUTPUT.
1827 C          JFLAG = 4 MEANS NO ROOTS (OF ODD MULTIPLICITY) WERE
1828 C                    FOUND IN (X0,X1) (NO SIGN CHANGES).
1829 C
1830 C X0,X1  = ENDPOINTS OF THE INTERVAL WHERE ROOTS ARE SOUGHT.
1831 C          X1 AND X0 ARE INPUT WHEN JFLAG = 0 (FIRST CALL), AND
1832 C          MUST BE LEFT UNCHANGED BETWEEN CALLS UNTIL THE PROBLEM IS
1833 C          COMPLETED.  X0 AND X1 MUST BE DISTINCT, BUT X1 - X0 MAY BE
1834 C          OF EITHER SIGN.  HOWEVER, THE NOTION OF -LEFT- AND -RIGHT-
1835 C          WILL BE USED TO MEAN NEARER TO X0 OR X1, RESPECTIVELY.
1836 C          WHEN JFLAG .GE. 2 ON RETURN, X0 AND X1 ARE OUTPUT, AND
1837 C          ARE THE ENDPOINTS OF THE RELEVANT INTERVAL.
1838 C
1839 C G0,G1  = ARRAYS OF LENGTH NG CONTAINING THE VECTORS G(X0) AND G(X1),
1840 C          RESPECTIVELY.  WHEN JFLAG = 0, G0 AND G1 ARE INPUT AND
1841 C          NONE OF THE G0(I) SHOULD BE BE ZERO.
1842 C          WHEN JFLAG .GE. 2 ON RETURN, G0 AND G1 ARE OUTPUT.
1843 C
1844 C GX     = ARRAY OF LENGTH NG CONTAINING G(X).  GX IS INPUT
1845 C          WHEN JFLAG = 1, AND OUTPUT WHEN JFLAG .GE. 2.
1846 C
1847 C X      = INDEPENDENT VARIABLE VALUE.  OUTPUT ONLY.
1848 C          WHEN JFLAG = 1 ON OUTPUT, X IS THE POINT AT WHICH G(X)
1849 C          IS TO BE EVALUATED AND LOADED INTO GX.
1850 C          WHEN JFLAG = 2 OR 3, X IS THE ROOT.
1851 C          WHEN JFLAG = 4, X IS THE RIGHT ENDPOINT OF THE INTERVAL, X1.
1852 C
1853 C JROOT  = INTEGER ARRAY OF LENGTH NG.  OUTPUT ONLY.
1854 C          WHEN JFLAG = 2 OR 3, JROOT INDICATES WHICH COMPONENTS
1855 C          OF G(X) HAVE A ROOT AT X.  JROOT(I) IS 1 IF THE I-TH
1856 C          COMPONENT HAS A ROOT, AND JROOT(I) = 0 OTHERWISE.
1857 C
1858 C IMAX, LAST, ALPHA, X2 =
1859 C          BOOKKEEPING VARIABLES WHICH MUST BE SAVED FROM CALL
1860 C          TO CALL.  THEY ARE SAVED INSIDE THE CALLING ROUTINE,
1861 C          BUT THEY ARE USED ONLY WITHIN THIS ROUTINE.
1862 C-----------------------------------------------------------------------
1863       INTEGER I, IMXOLD, NXLAST
1864       DOUBLE PRECISION T2, TMAX, ZERO
1865       LOGICAL ZROOT, SGNCHG, XROOT
1866       DATA ZERO/0.0D0/
1867 C
1868       IF (JFLAG .EQ. 1) GO TO 200
1869 C JFLAG .NE. 1.  CHECK FOR CHANGE IN SIGN OF G OR ZERO AT X1. ----------
1870       IMAX = 0
1871       TMAX = ZERO
1872       ZROOT = .FALSE.
1873       DO 120 I = 1,NG
1874         IF (DABS(G1(I)) .GT. ZERO) GO TO 110
1875         ZROOT = .TRUE.
1876         GO TO 120
1877 C AT THIS POINT, G0(I) HAS BEEN CHECKED AND CANNOT BE ZERO. ------------
1878  110    IF (DSIGN(1.0D0,G0(I)) .EQ. DSIGN(1.0D0,G1(I))) GO TO 120
1879           T2 = DABS(G1(I)/(G1(I)-G0(I)))
1880           IF (T2 .LE. TMAX) GO TO 120
1881             TMAX = T2
1882             IMAX = I
1883  120    CONTINUE
1884       IF (IMAX .GT. 0) GO TO 130
1885       SGNCHG = .FALSE.
1886       GO TO 140
1887  130  SGNCHG = .TRUE.
1888  140  IF (.NOT. SGNCHG) GO TO 400
1889 C THERE IS A SIGN CHANGE.  FIND THE FIRST ROOT IN THE INTERVAL. --------
1890       XROOT = .FALSE.
1891       NXLAST = 0
1892       LAST = 1
1893 C
1894 C REPEAT UNTIL THE FIRST ROOT IN THE INTERVAL IS FOUND.  LOOP POINT. ---
1895  150  CONTINUE
1896       IF (XROOT) GO TO 300
1897       IF (NXLAST .EQ. LAST) GO TO 160
1898       ALPHA = 1.0D0
1899       GO TO 180
1900  160  IF (LAST .EQ. 0) GO TO 170
1901       ALPHA = 0.5D0*ALPHA
1902       GO TO 180
1903  170  ALPHA = 2.0D0*ALPHA
1904  180  X2 = X1 - (X1-X0)*G1(IMAX)/(G1(IMAX) - ALPHA*G0(IMAX))
1905       IF ((DABS(X2-X0) .LT. HMIN) .AND.
1906      1   (DABS(X1-X0) .GT. 10.0D0*HMIN)) X2 = X0 + 0.1D0*(X1-X0)
1907       JFLAG = 1
1908       X = X2
1909 C RETURN TO THE CALLING ROUTINE TO GET A VALUE OF GX = G(X). -----------
1910       RETURN
1911 C CHECK TO SEE IN WHICH INTERVAL G CHANGES SIGN. -----------------------
1912  200  IMXOLD = IMAX
1913       IMAX = 0
1914       TMAX = ZERO
1915       ZROOT = .FALSE.
1916       DO 220 I = 1,NG
1917         IF (DABS(GX(I)) .GT. ZERO) GO TO 210
1918         ZROOT = .TRUE.
1919         GO TO 220
1920 C NEITHER G0(I) NOR GX(I) CAN BE ZERO AT THIS POINT. -------------------
1921  210    IF (DSIGN(1.0D0,G0(I)) .EQ. DSIGN(1.0D0,GX(I))) GO TO 220
1922           T2 = DABS(GX(I)/(GX(I) - G0(I)))
1923           IF (T2 .LE. TMAX) GO TO 220
1924             TMAX = T2
1925             IMAX = I
1926  220    CONTINUE
1927       IF (IMAX .GT. 0) GO TO 230
1928       SGNCHG = .FALSE.
1929       IMAX = IMXOLD
1930       GO TO 240
1931  230  SGNCHG = .TRUE.
1932  240  NXLAST = LAST
1933       IF (.NOT. SGNCHG) GO TO 250
1934 C SIGN CHANGE BETWEEN X0 AND X2, SO REPLACE X1 WITH X2. ----------------
1935       X1 = X2
1936       CALL DCOPY (NG, GX, 1, G1, 1)
1937       LAST = 1
1938       XROOT = .FALSE.
1939       GO TO 270
1940  250  IF (.NOT. ZROOT) GO TO 260
1941 C ZERO VALUE AT X2 AND NO SIGN CHANGE IN (X0,X2), SO X2 IS A ROOT. -----
1942       X1 = X2
1943       CALL DCOPY (NG, GX, 1, G1, 1)
1944       XROOT = .TRUE.
1945       GO TO 270
1946 C NO SIGN CHANGE BETWEEN X0 AND X2.  REPLACE X0 WITH X2. ---------------
1947  260  CONTINUE
1948       CALL DCOPY (NG, GX, 1, G0, 1)
1949       X0 = X2
1950       LAST = 0
1951       XROOT = .FALSE.
1952  270  IF (DABS(X1-X0) .LE. HMIN) XROOT = .TRUE.
1953       GO TO 150
1954 C
1955 C RETURN WITH X1 AS THE ROOT.  SET JROOT.  SET X = X1 AND GX = G1. -----
1956  300  JFLAG = 2
1957       X = X1
1958       CALL DCOPY (NG, G1, 1, GX, 1)
1959       DO 320 I = 1,NG
1960         JROOT(I) = 0
1961         IF (DABS(G1(I)) .GT. ZERO) GO TO 310
1962           JROOT(I) = 1
1963           GO TO 320
1964  310    IF (DSIGN(1.0D0,G0(I)) .NE. DSIGN(1.0D0,G1(I))) JROOT(I) = 1
1965  320    CONTINUE
1966       RETURN
1967 C
1968 C NO SIGN CHANGE IN THE INTERVAL.  CHECK FOR ZERO AT RIGHT ENDPOINT. ---
1969  400  IF (.NOT. ZROOT) GO TO 420
1970 C
1971 C ZERO VALUE AT X1 AND NO SIGN CHANGE IN (X0,X1).  RETURN JFLAG = 3. ---
1972       X = X1
1973       CALL DCOPY (NG, G1, 1, GX, 1)
1974       DO 410 I = 1,NG
1975         JROOT(I) = 0
1976         IF (DABS(G1(I)) .LE. ZERO) JROOT (I) = 1
1977  410  CONTINUE
1978       JFLAG = 3
1979       RETURN
1980 C
1981 C NO SIGN CHANGES IN THIS INTERVAL.  SET X = X1, RETURN JFLAG = 4. -----
1982  420  CALL DCOPY (NG, G1, 1, GX, 1)
1983       X = X1
1984       JFLAG = 4
1985       RETURN
1986 C---------------------- END OF SUBROUTINE DROOTS -----------------------
1987       END