/[ascend]/trunk/lsod/lsode.f
ViewVC logotype

Annotation of /trunk/lsod/lsode.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1 - (hide annotations) (download)
Fri Oct 29 20:54:12 2004 UTC (17 years, 1 month ago) by aw0a
File size: 118049 byte(s)
Setting up web subdirectory in repository
1 aw0a 1 subroutine lsode (f, neq, y, t, tout, itol, rtol, atol, itask,
2     1 istate, iopt, rwork, lrw, iwork, liw, jac, mf)
3     external f, jac, xascwv
4     integer neq, itol, itask, istate, iopt, lrw, iwork, liw, mf
5     double precision y, t, tout, rtol, atol, rwork
6     dimension neq(1), y(1), rtol(1), atol(1), rwork(lrw), iwork(liw)
7     c-----------------------------------------------------------------------
8     c this is the march 30, 1987 version of
9     c lsode.. livermore solver for ordinary differential equations.
10     c this version is in double precision.
11     c
12     c lsode solves the initial value problem for stiff or nonstiff
13     c systems of first order ode-s,
14     c dy/dt = f(t,y) , or, in component form,
15     c dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(neq)) (i = 1,...,neq).
16     c lsode is a package based on the gear and gearb packages, and on the
17     c october 23, 1978 version of the tentative odepack user interface
18     c standard, with minor modifications.
19     c-----------------------------------------------------------------------
20     c reference..
21     c alan c. hindmarsh, odepack, a systematized collection of ode
22     c solvers, in scientific computing, r. s. stepleman et al. (eds.),
23     c north-holland, amsterdam, 1983, pp. 55-64.
24     c-----------------------------------------------------------------------
25     c author and contact.. alan c. hindmarsh,
26     c computing and mathematics research div., l-316
27     c lawrence livermore national laboratory
28     c livermore, ca 94550.
29     c-----------------------------------------------------------------------
30     c summary of usage.
31     c
32     c communication between the user and the lsode package, for normal
33     c situations, is summarized here. this summary describes only a subset
34     c of the full set of options available. see the full description for
35     c details, including optional communication, nonstandard options,
36     c and instructions for special situations. see also the example
37     c problem (with program and output) following this summary.
38     c
39     c a. first provide a subroutine of the form..
40     c subroutine f (neq, t, y, ydot)
41     c dimension y(neq), ydot(neq)
42     c which supplies the vector function f by loading ydot(i) with f(i).
43     c
44     c b. next determine (or guess) whether or not the problem is stiff.
45     c stiffness occurs when the jacobian matrix df/dy has an eigenvalue
46     c whose real part is negative and large in magnitude, compared to the
47     c reciprocal of the t span of interest. if the problem is nonstiff,
48     c use a method flag mf = 10. if it is stiff, there are four standard
49     c choices for mf, and lsode requires the jacobian matrix in some form.
50     c this matrix is regarded either as full (mf = 21 or 22),
51     c or banded (mf = 24 or 25). in the banded case, lsode requires two
52     c half-bandwidth parameters ml and mu. these are, respectively, the
53     c widths of the lower and upper parts of the band, excluding the main
54     c diagonal. thus the band consists of the locations (i,j) with
55     c i-ml .le. j .le. i+mu, and the full bandwidth is ml+mu+1.
56     c
57     c c. if the problem is stiff, you are encouraged to supply the jacobian
58     c directly (mf = 21 or 24), but if this is not feasible, lsode will
59     c compute it internally by difference quotients (mf = 22 or 25).
60     c if you are supplying the jacobian, provide a subroutine of the form..
61     c subroutine jac (neq, t, y, ml, mu, pd, nrowpd)
62     c dimension y(neq), pd(nrowpd,neq)
63     c which supplies df/dy by loading pd as follows..
64     c for a full jacobian (mf = 21), load pd(i,j) with df(i)/dy(j),
65     c the partial derivative of f(i) with respect to y(j). (ignore the
66     c ml and mu arguments in this case.)
67     c for a banded jacobian (mf = 24), load pd(i-j+mu+1,j) with
68     c df(i)/dy(j), i.e. load the diagonal lines of df/dy into the rows of
69     c pd from the top down.
70     c in either case, only nonzero elements need be loaded.
71     c
72     c d. write a main program which calls subroutine lsode once for
73     c each point at which answers are desired. this should also provide
74     c for possible use of logical unit 6 for output of error messages
75     c by lsode. on the first call to lsode, supply arguments as follows..
76     c f = name of subroutine for right-hand side vector f.
77     c this name must be declared external in calling program.
78     c neq = number of first order ode-s.
79     c y = array of initial values, of length neq.
80     c t = the initial value of the independent variable.
81     c tout = first point where output is desired (.ne. t).
82     c itol = 1 or 2 according as atol (below) is a scalar or array.
83     c rtol = relative tolerance parameter (scalar).
84     c atol = absolute tolerance parameter (scalar or array).
85     c the estimated local error in y(i) will be controlled so as
86     c to be roughly less (in magnitude) than
87     c ewt(i) = rtol*abs(y(i)) + atol if itol = 1, or
88     c ewt(i) = rtol*abs(y(i)) + atol(i) if itol = 2.
89     c thus the local error test passes if, in each component,
90     c either the absolute error is less than atol (or atol(i)),
91     c or the relative error is less than rtol.
92     c use rtol = 0.0 for pure absolute error control, and
93     c use atol = 0.0 (or atol(i) = 0.0) for pure relative error
94     c control. caution.. actual (global) errors may exceed these
95     c local tolerances, so choose them conservatively.
96     c itask = 1 for normal computation of output values of y at t = tout.
97     c istate = integer flag (input and output). set istate = 1.
98     c iopt = 0 to indicate no optional inputs used.
99     c rwork = real work array of length at least..
100     c 20 + 16*neq for mf = 10,
101     c 22 + 9*neq + neq**2 for mf = 21 or 22,
102     c 22 + 10*neq + (2*ml + mu)*neq for mf = 24 or 25.
103     c lrw = declared length of rwork (in user-s dimension).
104     c iwork = integer work array of length at least..
105     c 20 for mf = 10,
106     c 20 + neq for mf = 21, 22, 24, or 25.
107     c if mf = 24 or 25, input in iwork(1),iwork(2) the lower
108     c and upper half-bandwidths ml,mu.
109     c liw = declared length of iwork (in user-s dimension).
110     c jac = name of subroutine for jacobian matrix (mf = 21 or 24).
111     c if used, this name must be declared external in calling
112     c program. if not used, pass a dummy name.
113     c mf = method flag. standard values are..
114     c 10 for nonstiff (adams) method, no jacobian used.
115     c 21 for stiff (bdf) method, user-supplied full jacobian.
116     c 22 for stiff method, internally generated full jacobian.
117     c 24 for stiff method, user-supplied banded jacobian.
118     c 25 for stiff method, internally generated banded jacobian.
119     c note that the main program must declare arrays y, rwork, iwork,
120     c and possibly atol.
121     c
122     c e. the output from the first call (or any call) is..
123     c y = array of computed values of y(t) vector.
124     c t = corresponding value of independent variable (normally tout).
125     c istate = 2 if lsode was successful, negative otherwise.
126     c -1 means excess work done on this call (perhaps wrong mf).
127     c -2 means excess accuracy requested (tolerances too small).
128     c -3 means illegal input detected (see printed message).
129     c -4 means repeated error test failures (check all inputs).
130     c -5 means repeated convergence failures (perhaps bad jacobian
131     c supplied or wrong choice of mf or tolerances).
132     c -6 means error weight became zero during problem. (solution
133     c component i vanished, and atol or atol(i) = 0.)
134     c
135     c f. to continue the integration after a successful return, simply
136     c reset tout and call lsode again. no other parameters need be reset.
137     c
138     c-----------------------------------------------------------------------
139     c example problem.
140     c
141     c the following is a simple example problem, with the coding
142     c needed for its solution by lsode. the problem is from chemical
143     c kinetics, and consists of the following three rate equations..
144     c dy1/dt = -.04*y1 + 1.e4*y2*y3
145     c dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
146     c dy3/dt = 3.e7*y2**2
147     c on the interval from t = 0.0 to t = 4.e10, with initial conditions
148     c y1 = 1.0, y2 = y3 = 0. the problem is stiff.
149     c
150     c the following coding solves this problem with lsode, using mf = 21
151     c and printing results at t = .4, 4., ..., 4.e10. it uses
152     c itol = 2 and atol much smaller for y2 than y1 or y3 because
153     c y2 has much smaller values.
154     c at the end of the run, statistical quantities of interest are
155     c printed (see optional outputs in the full description below).
156     c
157     c external fex, jex
158     c double precision atol, rtol, rwork, t, tout, y
159     c dimension y(3), atol(3), rwork(58), iwork(23)
160     c neq = 3
161     c y(1) = 1.d0
162     c y(2) = 0.d0
163     c y(3) = 0.d0
164     c t = 0.d0
165     c tout = .4d0
166     c itol = 2
167     c rtol = 1.d-4
168     c atol(1) = 1.d-6
169     c atol(2) = 1.d-10
170     c atol(3) = 1.d-6
171     c itask = 1
172     c istate = 1
173     c iopt = 0
174     c lrw = 58
175     c liw = 23
176     c mf = 21
177     c do 40 iout = 1,12
178     c call lsode(fex,neq,y,t,tout,itol,rtol,atol,itask,istate,
179     c 1 iopt,rwork,lrw,iwork,liw,jex,mf)
180     c write(6,20)t,y(1),y(2),y(3)
181     c 20 format(7h at t =,e12.4,6h y =,3e14.6)
182     c if (istate .lt. 0) go to 80
183     c 40 tout = tout*10.d0
184     c write(6,60)iwork(11),iwork(12),iwork(13)
185     c 60 format(/12h no. steps =,i4,11h no. f-s =,i4,11h no. j-s =,i4)
186     c stop
187     c 80 write(6,90)istate
188     c 90 format(///22h error halt.. istate =,i3)
189     c stop
190     c end
191     c
192     c subroutine fex (neq, t, y, ydot)
193     c double precision t, y, ydot
194     c dimension y(3), ydot(3)
195     c ydot(1) = -.04d0*y(1) + 1.d4*y(2)*y(3)
196     c ydot(3) = 3.d7*y(2)*y(2)
197     c ydot(2) = -ydot(1) - ydot(3)
198     c return
199     c end
200     c
201     c subroutine jex (neq, t, y, ml, mu, pd, nrpd)
202     c double precision pd, t, y
203     c dimension y(3), pd(nrpd,3)
204     c pd(1,1) = -.04d0
205     c pd(1,2) = 1.d4*y(3)
206     c pd(1,3) = 1.d4*y(2)
207     c pd(2,1) = .04d0
208     c pd(2,3) = -pd(1,3)
209     c pd(3,2) = 6.d7*y(2)
210     c pd(2,2) = -pd(1,2) - pd(3,2)
211     c return
212     c end
213     c
214     c the output of this program (on a cdc-7600 in single precision)
215     c is as follows..
216     c
217     c at t = 4.0000e-01 y = 9.851726e-01 3.386406e-05 1.479357e-02
218     c at t = 4.0000e+00 y = 9.055142e-01 2.240418e-05 9.446344e-02
219     c at t = 4.0000e+01 y = 7.158050e-01 9.184616e-06 2.841858e-01
220     c at t = 4.0000e+02 y = 4.504846e-01 3.222434e-06 5.495122e-01
221     c at t = 4.0000e+03 y = 1.831701e-01 8.940379e-07 8.168290e-01
222     c at t = 4.0000e+04 y = 3.897016e-02 1.621193e-07 9.610297e-01
223     c at t = 4.0000e+05 y = 4.935213e-03 1.983756e-08 9.950648e-01
224     c at t = 4.0000e+06 y = 5.159269e-04 2.064759e-09 9.994841e-01
225     c at t = 4.0000e+07 y = 5.306413e-05 2.122677e-10 9.999469e-01
226     c at t = 4.0000e+08 y = 5.494529e-06 2.197824e-11 9.999945e-01
227     c at t = 4.0000e+09 y = 5.129458e-07 2.051784e-12 9.999995e-01
228     c at t = 4.0000e+10 y = -7.170586e-08 -2.868234e-13 1.000000e+00
229     c
230     c no. steps = 330 no. f-s = 405 no. j-s = 69
231     c-----------------------------------------------------------------------
232     c full description of user interface to lsode.
233     c
234     c the user interface to lsode consists of the following parts.
235     c
236     c i. the call sequence to subroutine lsode, which is a driver
237     c routine for the solver. this includes descriptions of both
238     c the call sequence arguments and of user-supplied routines.
239     c following these descriptions is a description of
240     c optional inputs available through the call sequence, and then
241     c a description of optional outputs (in the work arrays).
242     c
243     c ii. descriptions of other routines in the lsode package that may be
244     c (optionally) called by the user. these provide the ability to
245     c alter error message handling, save and restore the internal
246     c common, and obtain specified derivatives of the solution y(t).
247     c
248     c iii. descriptions of common blocks to be declared in overlay
249     c or similar environments, or to be saved when doing an interrupt
250     c of the problem and continued solution later.
251     c
252     c iv. description of two routines in the lsode package, either of
253     c which the user may replace with his own version, if desired.
254     c these relate to the measurement of errors.
255     c
256     c-----------------------------------------------------------------------
257     c part i. call sequence.
258     c
259     c the call sequence parameters used for input only are
260     c f, neq, tout, itol, rtol, atol, itask, iopt, lrw, liw, jac, mf,
261     c and those used for both input and output are
262     c y, t, istate.
263     c the work arrays rwork and iwork are also used for conditional and
264     c optional inputs and optional outputs. (the term output here refers
265     c to the return from subroutine lsode to the user-s calling program.)
266     c
267     c the legality of input parameters will be thoroughly checked on the
268     c initial call for the problem, but not checked thereafter unless a
269     c change in input parameters is flagged by istate = 3 on input.
270     c
271     c the descriptions of the call arguments are as follows.
272     c
273     c f = the name of the user-supplied subroutine defining the
274     c ode system. the system must be put in the first-order
275     c form dy/dt = f(t,y), where f is a vector-valued function
276     c of the scalar t and the vector y. subroutine f is to
277     c compute the function f. it is to have the form
278     c subroutine f (neq, t, y, ydot)
279     c dimension y(1), ydot(1)
280     c where neq, t, and y are input, and the array ydot = f(t,y)
281     c is output. y and ydot are arrays of length neq.
282     c (in the dimension statement above, 1 is a dummy
283     c dimension.. it can be replaced by any value.)
284     c subroutine f should not alter y(1),...,y(neq).
285     c f must be declared external in the calling program.
286     c
287     c subroutine f may access user-defined quantities in
288     c neq(2),... and/or in y(neq(1)+1),... if neq is an array
289     c (dimensioned in f) and/or y has length exceeding neq(1).
290     c see the descriptions of neq and y below.
291     c
292     c if quantities computed in the f routine are needed
293     c externally to lsode, an extra call to f should be made
294     c for this purpose, for consistent and accurate results.
295     c if only the derivative dy/dt is needed, use intdy instead.
296     c
297     c neq = the size of the ode system (number of first order
298     c ordinary differential equations). used only for input.
299     c neq may be decreased, but not increased, during the problem.
300     c if neq is decreased (with istate = 3 on input), the
301     c remaining components of y should be left undisturbed, if
302     c these are to be accessed in f and/or jac.
303     c
304     c normally, neq is a scalar, and it is generally referred to
305     c as a scalar in this user interface description. however,
306     c neq may be an array, with neq(1) set to the system size.
307     c (the lsode package accesses only neq(1).) in either case,
308     c this parameter is passed as the neq argument in all calls
309     c to f and jac. hence, if it is an array, locations
310     c neq(2),... may be used to store other integer data and pass
311     c it to f and/or jac. subroutines f and/or jac must include
312     c neq in a dimension statement in that case.
313     c
314     c y = a real array for the vector of dependent variables, of
315     c length neq or more. used for both input and output on the
316     c first call (istate = 1), and only for output on other calls.
317     c on the first call, y must contain the vector of initial
318     c values. on output, y contains the computed solution vector,
319     c evaluated at t. if desired, the y array may be used
320     c for other purposes between calls to the solver.
321     c
322     c this array is passed as the y argument in all calls to
323     c f and jac. hence its length may exceed neq, and locations
324     c y(neq+1),... may be used to store other real data and
325     c pass it to f and/or jac. (the lsode package accesses only
326     c y(1),...,y(neq).)
327     c
328     c t = the independent variable. on input, t is used only on the
329     c first call, as the initial point of the integration.
330     c on output, after each call, t is the value at which a
331     c computed solution y is evaluated (usually the same as tout).
332     c on an error return, t is the farthest point reached.
333     c
334     c tout = the next value of t at which a computed solution is desired.
335     c used only for input.
336     c
337     c when starting the problem (istate = 1), tout may be equal
338     c to t for one call, then should .ne. t for the next call.
339     c for the initial t, an input value of tout .ne. t is used
340     c in order to determine the direction of the integration
341     c (i.e. the algebraic sign of the step sizes) and the rough
342     c scale of the problem. integration in either direction
343     c (forward or backward in t) is permitted.
344     c
345     c if itask = 2 or 5 (one-step modes), tout is ignored after
346     c the first call (i.e. the first call with tout .ne. t).
347     c otherwise, tout is required on every call.
348     c
349     c if itask = 1, 3, or 4, the values of tout need not be
350     c monotone, but a value of tout which backs up is limited
351     c to the current internal t interval, whose endpoints are
352     c tcur - hu and tcur (see optional outputs, below, for
353     c tcur and hu).
354     c
355     c itol = an indicator for the type of error control. see
356     c description below under atol. used only for input.
357     c
358     c rtol = a relative error tolerance parameter, either a scalar or
359     c an array of length neq. see description below under atol.
360     c input only.
361     c
362     c atol = an absolute error tolerance parameter, either a scalar or
363     c an array of length neq. input only.
364     c
365     c the input parameters itol, rtol, and atol determine
366     c the error control performed by the solver. the solver will
367     c control the vector e = (e(i)) of estimated local errors
368     c in y, according to an inequality of the form
369     c rms-norm of ( e(i)/ewt(i) ) .le. 1,
370     c where ewt(i) = rtol(i)*abs(y(i)) + atol(i),
371     c and the rms-norm (root-mean-square norm) here is
372     c rms-norm(v) = sqrt(sum v(i)**2 / neq). here ewt = (ewt(i))
373     c is a vector of weights which must always be positive, and
374     c the values of rtol and atol should all be non-negative.
375     c the following table gives the types (scalar/array) of
376     c rtol and atol, and the corresponding form of ewt(i).
377     c
378     c itol rtol atol ewt(i)
379     c 1 scalar scalar rtol*abs(y(i)) + atol
380     c 2 scalar array rtol*abs(y(i)) + atol(i)
381     c 3 array scalar rtol(i)*abs(y(i)) + atol
382     c 4 array array rtol(i)*abs(y(i)) + atol(i)
383     c
384     c when either of these parameters is a scalar, it need not
385     c be dimensioned in the user-s calling program.
386     c
387     c if none of the above choices (with itol, rtol, and atol
388     c fixed throughout the problem) is suitable, more general
389     c error controls can be obtained by substituting
390     c user-supplied routines for the setting of ewt and/or for
391     c the norm calculation. see part iv below.
392     c
393     c if global errors are to be estimated by making a repeated
394     c run on the same problem with smaller tolerances, then all
395     c components of rtol and atol (i.e. of ewt) should be scaled
396     c down uniformly.
397     c
398     c itask = an index specifying the task to be performed.
399     c input only. itask has the following values and meanings.
400     c 1 means normal computation of output values of y(t) at
401     c t = tout (by overshooting and interpolating).
402     c 2 means take one step only and return.
403     c 3 means stop at the first internal mesh point at or
404     c beyond t = tout and return.
405     c 4 means normal computation of output values of y(t) at
406     c t = tout but without overshooting t = tcrit.
407     c tcrit must be input as rwork(1). tcrit may be equal to
408     c or beyond tout, but not behind it in the direction of
409     c integration. this option is useful if the problem
410     c has a singularity at or beyond t = tcrit.
411     c 5 means take one step, without passing tcrit, and return.
412     c tcrit must be input as rwork(1).
413     c
414     c note.. if itask = 4 or 5 and the solver reaches tcrit
415     c (within roundoff), it will return t = tcrit (exactly) to
416     c indicate this (unless itask = 4 and tout comes before tcrit,
417     c in which case answers at t = tout are returned first).
418     c
419     c istate = an index used for input and output to specify the
420     c the state of the calculation.
421     c
422     c on input, the values of istate are as follows.
423     c 1 means this is the first call for the problem
424     c (initializations will be done). see note below.
425     c 2 means this is not the first call, and the calculation
426     c is to continue normally, with no change in any input
427     c parameters except possibly tout and itask.
428     c (if itol, rtol, and/or atol are changed between calls
429     c with istate = 2, the new values will be used but not
430     c tested for legality.)
431     c 3 means this is not the first call, and the
432     c calculation is to continue normally, but with
433     c a change in input parameters other than
434     c tout and itask. changes are allowed in
435     c neq, itol, rtol, atol, iopt, lrw, liw, mf, ml, mu,
436     c and any of the optional inputs except h0.
437     c (see iwork description for ml and mu.)
438     c note.. a preliminary call with tout = t is not counted
439     c as a first call here, as no initialization or checking of
440     c input is done. (such a call is sometimes useful for the
441     c purpose of outputting the initial conditions.)
442     c thus the first call for which tout .ne. t requires
443     c istate = 1 on input.
444     c
445     c on output, istate has the following values and meanings.
446     c 1 means nothing was done, as tout was equal to t with
447     c istate = 1 on input. (however, an internal counter was
448     c set to detect and prevent repeated calls of this type.)
449     c 2 means the integration was performed successfully.
450     c -1 means an excessive amount of work (more than mxstep
451     c steps) was done on this call, before completing the
452     c requested task, but the integration was otherwise
453     c successful as far as t. (mxstep is an optional input
454     c and is normally 500.) to continue, the user may
455     c simply reset istate to a value .gt. 1 and call again
456     c (the excess work step counter will be reset to 0).
457     c in addition, the user may increase mxstep to avoid
458     c this error return (see below on optional inputs).
459     c -2 means too much accuracy was requested for the precision
460     c of the machine being used. this was detected before
461     c completing the requested task, but the integration
462     c was successful as far as t. to continue, the tolerance
463     c parameters must be reset, and istate must be set
464     c to 3. the optional output tolsf may be used for this
465     c purpose. (note.. if this condition is detected before
466     c taking any steps, then an illegal input return
467     c (istate = -3) occurs instead.)
468     c -3 means illegal input was detected, before taking any
469     c integration steps. see written message for details.
470     c note.. if the solver detects an infinite loop of calls
471     c to the solver with illegal input, it will cause
472     c the run to stop.
473     c -4 means there were repeated error test failures on
474     c one attempted step, before completing the requested
475     c task, but the integration was successful as far as t.
476     c the problem may have a singularity, or the input
477     c may be inappropriate.
478     c -5 means there were repeated convergence test failures on
479     c one attempted step, before completing the requested
480     c task, but the integration was successful as far as t.
481     c this may be caused by an inaccurate jacobian matrix,
482     c if one is being used.
483     c -6 means ewt(i) became zero for some i during the
484     c integration. pure relative error control (atol(i)=0.0)
485     c was requested on a variable which has now vanished.
486     c the integration was successful as far as t.
487     c
488     c note.. since the normal output value of istate is 2,
489     c it does not need to be reset for normal continuation.
490     c also, since a negative input value of istate will be
491     c regarded as illegal, a negative output value requires the
492     c user to change it, and possibly other inputs, before
493     c calling the solver again.
494     c
495     c iopt = an integer flag to specify whether or not any optional
496     c inputs are being used on this call. input only.
497     c the optional inputs are listed separately below.
498     c iopt = 0 means no optional inputs are being used.
499     c default values will be used in all cases.
500     c iopt = 1 means one or more optional inputs are being used.
501     c
502     c rwork = a real working array (double precision).
503     c the length of rwork must be at least
504     c 20 + nyh*(maxord + 1) + 3*neq + lwm where
505     c nyh = the initial value of neq,
506     c maxord = 12 (if meth = 1) or 5 (if meth = 2) (unless a
507     c smaller value is given as an optional input),
508     c lwm = 0 if miter = 0,
509     c lwm = neq**2 + 2 if miter is 1 or 2,
510     c lwm = neq + 2 if miter = 3, and
511     c lwm = (2*ml+mu+1)*neq + 2 if miter is 4 or 5.
512     c (see the mf description for meth and miter.)
513     c thus if maxord has its default value and neq is constant,
514     c this length is..
515     c 20 + 16*neq for mf = 10,
516     c 22 + 16*neq + neq**2 for mf = 11 or 12,
517     c 22 + 17*neq for mf = 13,
518     c 22 + 17*neq + (2*ml+mu)*neq for mf = 14 or 15,
519     c 20 + 9*neq for mf = 20,
520     c 22 + 9*neq + neq**2 for mf = 21 or 22,
521     c 22 + 10*neq for mf = 23,
522     c 22 + 10*neq + (2*ml+mu)*neq for mf = 24 or 25.
523     c the first 20 words of rwork are reserved for conditional
524     c and optional inputs and optional outputs.
525     c
526     c the following word in rwork is a conditional input..
527     c rwork(1) = tcrit = critical value of t which the solver
528     c is not to overshoot. required if itask is
529     c 4 or 5, and ignored otherwise. (see itask.)
530     c
531     c lrw = the length of the array rwork, as declared by the user.
532     c (this will be checked by the solver.)
533     c
534     c iwork = an integer work array. the length of iwork must be at least
535     c 20 if miter = 0 or 3 (mf = 10, 13, 20, 23), or
536     c 20 + neq otherwise (mf = 11, 12, 14, 15, 21, 22, 24, 25).
537     c the first few words of iwork are used for conditional and
538     c optional inputs and optional outputs.
539     c
540     c the following 2 words in iwork are conditional inputs..
541     c iwork(1) = ml these are the lower and upper
542     c iwork(2) = mu half-bandwidths, respectively, of the
543     c banded jacobian, excluding the main diagonal.
544     c the band is defined by the matrix locations
545     c (i,j) with i-ml .le. j .le. i+mu. ml and mu
546     c must satisfy 0 .le. ml,mu .le. neq-1.
547     c these are required if miter is 4 or 5, and
548     c ignored otherwise. ml and mu may in fact be
549     c the band parameters for a matrix to which
550     c df/dy is only approximately equal.
551     c
552     c liw = the length of the array iwork, as declared by the user.
553     c (this will be checked by the solver.)
554     c
555     c note.. the work arrays must not be altered between calls to lsode
556     c for the same problem, except possibly for the conditional and
557     c optional inputs, and except for the last 3*neq words of rwork.
558     c the latter space is used for internal scratch space, and so is
559     c available for use by the user outside lsode between calls, if
560     c desired (but not for use by f or jac).
561     c
562     c jac = the name of the user-supplied routine (miter = 1 or 4) to
563     c compute the jacobian matrix, df/dy, as a function of
564     c the scalar t and the vector y. it is to have the form
565     c subroutine jac (neq, t, y, ml, mu, pd, nrowpd)
566     c dimension y(1), pd(nrowpd,1)
567     c where neq, t, y, ml, mu, and nrowpd are input and the array
568     c pd is to be loaded with partial derivatives (elements of
569     c the jacobian matrix) on output. pd must be given a first
570     c dimension of nrowpd. t and y have the same meaning as in
571     c subroutine f. (in the dimension statement above, 1 is a
572     c dummy dimension.. it can be replaced by any value.)
573     c in the full matrix case (miter = 1), ml and mu are
574     c ignored, and the jacobian is to be loaded into pd in
575     c columnwise manner, with df(i)/dy(j) loaded into pd(i,j).
576     c in the band matrix case (miter = 4), the elements
577     c within the band are to be loaded into pd in columnwise
578     c manner, with diagonal lines of df/dy loaded into the rows
579     c of pd. thus df(i)/dy(j) is to be loaded into pd(i-j+mu+1,j).
580     c ml and mu are the half-bandwidth parameters (see iwork).
581     c the locations in pd in the two triangular areas which
582     c correspond to nonexistent matrix elements can be ignored
583     c or loaded arbitrarily, as they are overwritten by lsode.
584     c jac need not provide df/dy exactly. a crude
585     c approximation (possibly with a smaller bandwidth) will do.
586     c in either case, pd is preset to zero by the solver,
587     c so that only the nonzero elements need be loaded by jac.
588     c each call to jac is preceded by a call to f with the same
589     c arguments neq, t, and y. thus to gain some efficiency,
590     c intermediate quantities shared by both calculations may be
591     c saved in a user common block by f and not recomputed by jac,
592     c if desired. also, jac may alter the y array, if desired.
593     c jac must be declared external in the calling program.
594     c subroutine jac may access user-defined quantities in
595     c neq(2),... and/or in y(neq(1)+1),... if neq is an array
596     c (dimensioned in jac) and/or y has length exceeding neq(1).
597     c see the descriptions of neq and y above.
598     c
599     c mf = the method flag. used only for input. the legal values of
600     c mf are 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24, and 25.
601     c mf has decimal digits meth and miter.. mf = 10*meth + miter.
602     c meth indicates the basic linear multistep method..
603     c meth = 1 means the implicit adams method.
604     c meth = 2 means the method based on backward
605     c differentiation formulas (bdf-s).
606     c miter indicates the corrector iteration method..
607     c miter = 0 means functional iteration (no jacobian matrix
608     c is involved).
609     c miter = 1 means chord iteration with a user-supplied
610     c full (neq by neq) jacobian.
611     c miter = 2 means chord iteration with an internally
612     c generated (difference quotient) full jacobian
613     c (using neq extra calls to f per df/dy value).
614     c miter = 3 means chord iteration with an internally
615     c generated diagonal jacobian approximation.
616     c (using 1 extra call to f per df/dy evaluation).
617     c miter = 4 means chord iteration with a user-supplied
618     c banded jacobian.
619     c miter = 5 means chord iteration with an internally
620     c generated banded jacobian (using ml+mu+1 extra
621     c calls to f per df/dy evaluation).
622     c if miter = 1 or 4, the user must supply a subroutine jac
623     c (the name is arbitrary) as described above under jac.
624     c for other values of miter, a dummy argument can be used.
625     c-----------------------------------------------------------------------
626     c optional inputs.
627     c
628     c the following is a list of the optional inputs provided for in the
629     c call sequence. (see also part ii.) for each such input variable,
630     c this table lists its name as used in this documentation, its
631     c location in the call sequence, its meaning, and the default value.
632     c the use of any of these inputs requires iopt = 1, and in that
633     c case all of these inputs are examined. a value of zero for any
634     c of these optional inputs will cause the default value to be used.
635     c thus to use a subset of the optional inputs, simply preload
636     c locations 5 to 10 in rwork and iwork to 0.0 and 0 respectively, and
637     c then set those of interest to nonzero values.
638     c
639     c name location meaning and default value
640     c
641     c h0 rwork(5) the step size to be attempted on the first step.
642     c the default value is determined by the solver.
643     c
644     c hmax rwork(6) the maximum absolute step size allowed.
645     c the default value is infinite.
646     c
647     c hmin rwork(7) the minimum absolute step size allowed.
648     c the default value is 0. (this lower bound is not
649     c enforced on the final step before reaching tcrit
650     c when itask = 4 or 5.)
651     c
652     c maxord iwork(5) the maximum order to be allowed. the default
653     c value is 12 if meth = 1, and 5 if meth = 2.
654     c if maxord exceeds the default value, it will
655     c be reduced to the default value.
656     c if maxord is changed during the problem, it may
657     c cause the current order to be reduced.
658     c
659     c mxstep iwork(6) maximum number of (internally defined) steps
660     c allowed during one call to the solver.
661     c the default value is 500.
662     c
663     c mxhnil iwork(7) maximum number of messages printed (per problem)
664     c warning that t + h = t on a step (h = step size).
665     c this must be positive to result in a non-default
666     c value. the default value is 10.
667     c-----------------------------------------------------------------------
668     c optional outputs.
669     c
670     c as optional additional output from lsode, the variables listed
671     c below are quantities related to the performance of lsode
672     c which are available to the user. these are communicated by way of
673     c the work arrays, but also have internal mnemonic names as shown.
674     c except where stated otherwise, all of these outputs are defined
675     c on any successful return from lsode, and on any return with
676     c istate = -1, -2, -4, -5, or -6. on an illegal input return
677     c (istate = -3), they will be unchanged from their existing values
678     c (if any), except possibly for tolsf, lenrw, and leniw.
679     c on any error return, outputs relevant to the error will be defined,
680     c as noted below.
681     c
682     c name location meaning
683     c
684     c hu rwork(11) the step size in t last used (successfully).
685     c
686     c hcur rwork(12) the step size to be attempted on the next step.
687     c
688     c tcur rwork(13) the current value of the independent variable
689     c which the solver has actually reached, i.e. the
690     c current internal mesh point in t. on output, tcur
691     c will always be at least as far as the argument
692     c t, but may be farther (if interpolation was done).
693     c
694     c tolsf rwork(14) a tolerance scale factor, greater than 1.0,
695     c computed when a request for too much accuracy was
696     c detected (istate = -3 if detected at the start of
697     c the problem, istate = -2 otherwise). if itol is
698     c left unaltered but rtol and atol are uniformly
699     c scaled up by a factor of tolsf for the next call,
700     c then the solver is deemed likely to succeed.
701     c (the user may also ignore tolsf and alter the
702     c tolerance parameters in any other way appropriate.)
703     c
704     c nst iwork(11) the number of steps taken for the problem so far.
705     c
706     c nfe iwork(12) the number of f evaluations for the problem so far.
707     c
708     c nje iwork(13) the number of jacobian evaluations (and of matrix
709     c lu decompositions) for the problem so far.
710     c
711     c nqu iwork(14) the method order last used (successfully).
712     c
713     c nqcur iwork(15) the order to be attempted on the next step.
714     c
715     c imxer iwork(16) the index of the component of largest magnitude in
716     c the weighted local error vector ( e(i)/ewt(i) ),
717     c on an error return with istate = -4 or -5.
718     c
719     c lenrw iwork(17) the length of rwork actually required.
720     c this is defined on normal returns and on an illegal
721     c input return for insufficient storage.
722     c
723     c leniw iwork(18) the length of iwork actually required.
724     c this is defined on normal returns and on an illegal
725     c input return for insufficient storage.
726     c
727     c the following two arrays are segments of the rwork array which
728     c may also be of interest to the user as optional outputs.
729     c for each array, the table below gives its internal name,
730     c its base address in rwork, and its description.
731     c
732     c name base address description
733     c
734     c yh 21 the nordsieck history array, of size nyh by
735     c (nqcur + 1), where nyh is the initial value
736     c of neq. for j = 0,1,...,nqcur, column j+1
737     c of yh contains hcur**j/factorial(j) times
738     c the j-th derivative of the interpolating
739     c polynomial currently representing the solution,
740     c evaluated at t = tcur.
741     c
742     c acor lenrw-neq+1 array of size neq used for the accumulated
743     c corrections on each step, scaled on output
744     c to represent the estimated local error in y
745     c on the last step. this is the vector e in
746     c the description of the error control. it is
747     c defined only on a successful return from lsode.
748     c
749     c-----------------------------------------------------------------------
750     c part ii. other routines callable.
751     c
752     c the following are optional calls which the user may make to
753     c gain additional capabilities in conjunction with lsode.
754     c (the routines xsetun and xsetf are designed to conform to the
755     c slatec error handling package.)
756     c
757     c form of call function
758     c call xsetun(lun) set the logical unit number, lun, for
759     c output of messages from lsode, if
760     c the default is not desired.
761     c the default value of lun is 6.
762     c
763     c call xsetf(mflag) set a flag to control the printing of
764     c messages by lsode.
765     c mflag = 0 means do not print. (danger..
766     c this risks losing valuable information.)
767     c mflag = 1 means print (the default).
768     c
769     c either of the above calls may be made at
770     c any time and will take effect immediately.
771     c
772     c call srcom(rsav,isav,job) saves and restores the contents of
773     c the internal common blocks used by
774     c lsode (see part iii below).
775     c rsav must be a real array of length 218
776     c or more, and isav must be an integer
777     c array of length 41 or more.
778     c job=1 means save common into rsav/isav.
779     c job=2 means restore common from rsav/isav.
780     c srcom is useful if one is
781     c interrupting a run and restarting
782     c later, or alternating between two or
783     c more problems solved with lsode.
784     c
785     c call intdy(,,,,,) provide derivatives of y, of various
786     c (see below) orders, at a specified point t, if
787     c desired. it may be called only after
788     c a successful return from lsode.
789     c
790     c the detailed instructions for using intdy are as follows.
791     c the form of the call is..
792     c
793     c call intdy (t, k, rwork(21), nyh, dky, iflag)
794     c
795     c the input parameters are..
796     c
797     c t = value of independent variable where answers are desired
798     c (normally the same as the t last returned by lsode).
799     c for valid results, t must lie between tcur - hu and tcur.
800     c (see optional outputs for tcur and hu.)
801     c k = integer order of the derivative desired. k must satisfy
802     c 0 .le. k .le. nqcur, where nqcur is the current order
803     c (see optional outputs). the capability corresponding
804     c to k = 0, i.e. computing y(t), is already provided
805     c by lsode directly. since nqcur .ge. 1, the first
806     c derivative dy/dt is always available with intdy.
807     c rwork(21) = the base address of the history array yh.
808     c nyh = column length of yh, equal to the initial value of neq.
809     c
810     c the output parameters are..
811     c
812     c dky = a real array of length neq containing the computed value
813     c of the k-th derivative of y(t).
814     c iflag = integer flag, returned as 0 if k and t were legal,
815     c -1 if k was illegal, and -2 if t was illegal.
816     c on an error return, a message is also written.
817     c-----------------------------------------------------------------------
818     c part iii. common blocks.
819     c
820     c if lsode is to be used in an overlay situation, the user
821     c must declare, in the primary overlay, the variables in..
822     c (1) the call sequence to lsode,
823     c (2) the two internal common blocks
824     c /ls0001/ of length 257 (218 double precision words
825     c followed by 39 integer words),
826     c /eh0001/ of length 2 (integer words).
827     c
828     c if lsode is used on a system in which the contents of internal
829     c common blocks are not preserved between calls, the user should
830     c declare the above two common blocks in his main program to insure
831     c that their contents are preserved.
832     c
833     c if the solution of a given problem by lsode is to be interrupted
834     c and then later continued, such as when restarting an interrupted run
835     c or alternating between two or more problems, the user should save,
836     c following the return from the last lsode call prior to the
837     c interruption, the contents of the call sequence variables and the
838     c internal common blocks, and later restore these values before the
839     c next lsode call for that problem. to save and restore the common
840     c blocks, use subroutine srcom (see part ii above).
841     c
842     c-----------------------------------------------------------------------
843     c part iv. optionally replaceable solver routines.
844     c
845     c below are descriptions of two routines in the lsode package which
846     c relate to the measurement of errors. either routine can be
847     c replaced by a user-supplied version, if desired. however, since such
848     c a replacement may have a major impact on performance, it should be
849     c done only when absolutely necessary, and only with great caution.
850     c (note.. the means by which the package version of a routine is
851     c superseded by the user-s version may be system-dependent.)
852     c
853     c (a) ewset.
854     c the following subroutine is called just before each internal
855     c integration step, and sets the array of error weights, ewt, as
856     c described under itol/rtol/atol above..
857     c subroutine ewset (neq, itol, rtol, atol, ycur, ewt)
858     c where neq, itol, rtol, and atol are as in the lsode call sequence,
859     c ycur contains the current dependent variable vector, and
860     c ewt is the array of weights set by ewset.
861     c
862     c if the user supplies this subroutine, it must return in ewt(i)
863     c (i = 1,...,neq) a positive quantity suitable for comparing errors
864     c in y(i) to. the ewt array returned by ewset is passed to the
865     c vnorm routine (see below), and also used by lsode in the computation
866     c of the optional output imxer, the diagonal jacobian approximation,
867     c and the increments for difference quotient jacobians.
868     c
869     c in the user-supplied version of ewset, it may be desirable to use
870     c the current values of derivatives of y. derivatives up to order nq
871     c are available from the history array yh, described above under
872     c optional outputs. in ewset, yh is identical to the ycur array,
873     c extended to nq + 1 columns with a column length of nyh and scale
874     c factors of h**j/factorial(j). on the first call for the problem,
875     c given by nst = 0, nq is 1 and h is temporarily set to 1.0.
876     c the quantities nq, nyh, h, and nst can be obtained by including
877     c in ewset the statements..
878     c double precision h, rls
879     c common /ls0001/ rls(218),ils(39)
880     c nq = ils(35)
881     c nyh = ils(14)
882     c nst = ils(36)
883     c h = rls(212)
884     c thus, for example, the current value of dy/dt can be obtained as
885     c ycur(nyh+i)/h (i=1,...,neq) (and the division by h is
886     c unnecessary when nst = 0).
887     c
888     c (b) vnorm.
889     c the following is a real function routine which computes the weighted
890     c root-mean-square norm of a vector v..
891     c d = vnorm (n, v, w)
892     c where..
893     c n = the length of the vector,
894     c v = real array of length n containing the vector,
895     c w = real array of length n containing weights,
896     c d = sqrt( (1/n) * sum(v(i)*w(i))**2 ).
897     c vnorm is called with n = neq and with w(i) = 1.0/ewt(i), where
898     c ewt is as set by subroutine ewset.
899     c
900     c if the user supplies this function, it should return a non-negative
901     c value of vnorm suitable for use in the error control in lsode.
902     c none of the arguments should be altered by vnorm.
903     c for example, a user-supplied vnorm routine might..
904     c -substitute a max-norm of (v(i)*w(i)) for the rms-norm, or
905     c -ignore some components of v in the norm, with the effect of
906     c suppressing the error control on those components of y.
907     c-----------------------------------------------------------------------
908     c-----------------------------------------------------------------------
909     c other routines in the lsode package.
910     c
911     c in addition to subroutine lsode, the lsode package includes the
912     c following subroutines and function routines..
913     c intdy computes an interpolated value of the y vector at t = tout.
914     c stode is the core integrator, which does one step of the
915     c integration and the associated error control.
916     c cfode sets all method coefficients and test constants.
917     c prepj computes and preprocesses the jacobian matrix j = df/dy
918     c and the newton iteration matrix p = i - h*l0*j.
919     c solsy manages solution of linear system in chord iteration.
920     c ewset sets the error weight vector ewt before each step.
921     c vnorm computes the weighted r.m.s. norm of a vector.
922     c srcom is a user-callable routine to save and restore
923     c the contents of the internal common blocks.
924     c dgefa and dgesl are routines from linpack for solving full
925     c systems of linear algebraic equations.
926     c dgbfa and dgbsl are routines from linpack for solving banded
927     c linear systems.
928     c daxpy, dscal, idamax, and ddot are basic linear algebra modules
929     c (blas) used by the above linpack routines.
930     c d1mach computes the unit roundoff in a machine-independent manner.
931     c XERRWV, xsetun, and xsetf handle the printing of all error
932     c messages and warnings. XERRWV is machine-dependent.
933     C ASCEND C XERRWV has been replaced by XASCWV from interface/Lsode.c.
934     c note.. vnorm, idamax, ddot, and d1mach are function routines.
935     c all the others are subroutines.
936     c
937     c the intrinsic and external routines used by lsode are..
938     c dabs, dmax1, dmin1, dfloat, max0, min0, mod, dsign, dsqrt, and write.
939     c
940     c a block data subprogram is also included with the package,
941     c for loading some of the variables in internal common.
942     c
943     c-----------------------------------------------------------------------
944     c the following card is for optimized compilation on llnl compilers.
945     clll. optimize
946     c-----------------------------------------------------------------------
947     cascend changes
948     Ckaa external aftime
949     external prepj, solsy
950     integer illin, init, lyh, lewt, lacor, lsavf, lwm, liwm,
951     1 mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns
952     integer icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
953     1 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
954     integer i, i1, i2, iflag, imxer, kgo, lf0,
955     1 leniw, lenrw, lenwm, ml, mord, mu, mxhnl0, mxstp0
956     double precision rowns,
957     1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
958     double precision atoli, ayi, big, ewti, h0, hmax, hmx, rh, rtoli,
959     1 tcrit, tdist, tnext, tol, tolsf, tp, size, sum, w0,
960     2 d1mach, vnorm
961     dimension mord(2)
962     logical ihit
963     c-----------------------------------------------------------------------
964     c the following internal common block contains
965     c (a) variables which are local to any subroutine but whose values must
966     c be preserved between calls to the routine (own variables), and
967     c (b) variables which are communicated between subroutines.
968     c the structure of the block is as follows.. all real variables are
969     c listed first, followed by all integers. within each type, the
970     c variables are grouped with those local to subroutine lsode first,
971     c then those local to subroutine stode, and finally those used
972     c for communication. the block is declared in subroutines
973     c lsode, intdy, stode, prepj, and solsy. groups of variables are
974     c replaced by dummy arrays in the common declarations in routines
975     c where those variables are not used.
976     c-----------------------------------------------------------------------
977     common /ls0001/ rowns(209),
978     1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
979     2 illin, init, lyh, lewt, lacor, lsavf, lwm, liwm,
980     3 mxstep, mxhnil, nhnil, ntrep, nslast, nyh, iowns(6),
981     4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
982     5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
983     c
984     data mord(1),mord(2)/12,5/, mxstp0/500/, mxhnl0/10/
985    
986     C ascend debugging options
987     C type *, 'iopt',iopt
988     C type *, 'h0',rwork(5)
989     C type *, 'hmax',rwork(6)
990     C type *, 'hmin',rwork(7)
991     C type *, 'maxs',iwork(6)
992     c-----------------------------------------------------------------------
993     c block a.
994     c this code block is executed on every call.
995     c it tests istate and itask for legality and branches appropriately.
996     c if istate .gt. 1 but the flag init shows that initialization has
997     c not yet been done, an error return occurs.
998     c if istate = 1 and tout = t, jump to block g and return immediately.
999     c-----------------------------------------------------------------------
1000     if (istate .lt. 1 .or. istate .gt. 3) go to 601
1001     if (itask .lt. 1 .or. itask .gt. 5) go to 602
1002     if (istate .eq. 1) go to 10
1003     if (init .eq. 0) go to 603
1004     if (istate .eq. 2) go to 200
1005     go to 20
1006     10 init = 0
1007     if (tout .eq. t) go to 430
1008     20 ntrep = 0
1009     c-----------------------------------------------------------------------
1010     c block b.
1011     c the next code block is executed for the initial call (istate = 1),
1012     c or for a continuation call with parameter changes (istate = 3).
1013     c it contains checking of all inputs and various initializations.
1014     c
1015     c first check legality of the non-optional inputs neq, itol, iopt,
1016     c mf, ml, and mu.
1017     c-----------------------------------------------------------------------
1018     if (neq(1) .le. 0) go to 604
1019     if (istate .eq. 1) go to 25
1020     if (neq(1) .gt. n) go to 605
1021     25 n = neq(1)
1022     if (itol .lt. 1 .or. itol .gt. 4) go to 606
1023     if (iopt .lt. 0 .or. iopt .gt. 1) go to 607
1024     meth = mf/10
1025     miter = mf - 10*meth
1026     if (meth .lt. 1 .or. meth .gt. 2) go to 608
1027     if (miter .lt. 0 .or. miter .gt. 5) go to 608
1028     if (miter .le. 3) go to 30
1029     ml = iwork(1)
1030     mu = iwork(2)
1031     if (ml .lt. 0 .or. ml .ge. n) go to 609
1032     if (mu .lt. 0 .or. mu .ge. n) go to 610
1033     30 continue
1034     c next process and check the optional inputs. --------------------------
1035     if (iopt .eq. 1) go to 40
1036     maxord = mord(meth)
1037     mxstep = mxstp0
1038     mxhnil = mxhnl0
1039     if (istate .eq. 1) h0 = 0.0d0
1040     hmxi = 0.0d0
1041     hmin = 0.0d0
1042     go to 60
1043     40 maxord = iwork(5)
1044     if (maxord .lt. 0) go to 611
1045     if (maxord .eq. 0) maxord = 100
1046     maxord = min0(maxord,mord(meth))
1047     mxstep = iwork(6)
1048     if (mxstep .lt. 0) go to 612
1049     if (mxstep .eq. 0) mxstep = mxstp0
1050     mxhnil = iwork(7)
1051     if (mxhnil .lt. 0) go to 613
1052     if (mxhnil .eq. 0) mxhnil = mxhnl0
1053     if (istate .ne. 1) go to 50
1054     h0 = rwork(5)
1055     if ((tout - t)*h0 .lt. 0.0d0) go to 614
1056     50 hmax = rwork(6)
1057     if (hmax .lt. 0.0d0) go to 615
1058     hmxi = 0.0d0
1059     if (hmax .gt. 0.0d0) hmxi = 1.0d0/hmax
1060     hmin = rwork(7)
1061     if (hmin .lt. 0.0d0) go to 616
1062     c-----------------------------------------------------------------------
1063     c set work array pointers and check lengths lrw and liw.
1064     c pointers to segments of rwork and iwork are named by prefixing l to
1065     c the name of the segment. e.g., the segment yh starts at rwork(lyh).
1066     c segments of rwork (in order) are denoted yh, wm, ewt, savf, acor.
1067     c-----------------------------------------------------------------------
1068     60 lyh = 21
1069     if (istate .eq. 1) nyh = n
1070     lwm = lyh + (maxord + 1)*nyh
1071     if (miter .eq. 0) lenwm = 0
1072     if (miter .eq. 1 .or. miter .eq. 2) lenwm = n*n + 2
1073     if (miter .eq. 3) lenwm = n + 2
1074     if (miter .ge. 4) lenwm = (2*ml + mu + 1)*n + 2
1075     lewt = lwm + lenwm
1076     lsavf = lewt + n
1077     lacor = lsavf + n
1078     lenrw = lacor + n - 1
1079     iwork(17) = lenrw
1080     liwm = 1
1081     leniw = 20 + n
1082     if (miter .eq. 0 .or. miter .eq. 3) leniw = 20
1083     iwork(18) = leniw
1084     if (lenrw .gt. lrw) go to 617
1085     if (leniw .gt. liw) go to 618
1086     c check rtol and atol for legality. ------------------------------------
1087     rtoli = rtol(1)
1088     atoli = atol(1)
1089     do 70 i = 1,n
1090     if (itol .ge. 3) rtoli = rtol(i)
1091     if (itol .eq. 2 .or. itol .eq. 4) atoli = atol(i)
1092     if (rtoli .lt. 0.0d0) go to 619
1093     if (atoli .lt. 0.0d0) go to 620
1094     70 continue
1095     if (istate .eq. 1) go to 100
1096     c if istate = 3, set flag to signal parameter changes to stode. --------
1097     jstart = -1
1098     if (nq .le. maxord) go to 90
1099     c maxord was reduced below nq. copy yh(*,maxord+2) into savf. ---------
1100     do 80 i = 1,n
1101     80 rwork(i+lsavf-1) = rwork(i+lwm-1)
1102     c reload wm(1) = rwork(lwm), since lwm may have changed. ---------------
1103     90 if (miter .gt. 0) rwork(lwm) = dsqrt(uround)
1104     if (n .eq. nyh) go to 200
1105     c neq was reduced. zero part of yh to avoid undefined references. -----
1106     i1 = lyh + l*nyh
1107     i2 = lyh + (maxord + 1)*nyh - 1
1108     if (i1 .gt. i2) go to 200
1109     do 95 i = i1,i2
1110     95 rwork(i) = 0.0d0
1111     go to 200
1112     c-----------------------------------------------------------------------
1113     c block c.
1114     c the next block is for the initial call only (istate = 1).
1115     c it contains all remaining initializations, the initial call to f,
1116     c and the calculation of the initial step size.
1117     c the error weights in ewt are inverted after being loaded.
1118     c-----------------------------------------------------------------------
1119     100 uround = d1mach(4)
1120     tn = t
1121     if (itask .ne. 4 .and. itask .ne. 5) go to 110
1122     tcrit = rwork(1)
1123     if ((tcrit - tout)*(tout - t) .lt. 0.0d0) go to 625
1124     if (h0 .ne. 0.0d0 .and. (t + h0 - tcrit)*h0 .gt. 0.0d0)
1125     1 h0 = tcrit - t
1126     110 jstart = 0
1127     if (miter .gt. 0) rwork(lwm) = dsqrt(uround)
1128     nhnil = 0
1129     nst = 0
1130     nje = 0
1131     nslast = 0
1132     hu = 0.0d0
1133     nqu = 0
1134     ccmax = 0.3d0
1135     maxcor = 3
1136     msbp = 20
1137     mxncf = 10
1138     c initial call to f. (lf0 points to yh(*,2).) -------------------------
1139     lf0 = lyh + nyh
1140     call f (neq, t, y, rwork(lf0))
1141     nfe = 1
1142     c load the initial value vector in yh. ---------------------------------
1143     do 115 i = 1,n
1144     115 rwork(i+lyh-1) = y(i)
1145     c load and invert the ewt array. (h is temporarily set to 1.0.) -------
1146     nq = 1
1147     h = 1.0d0
1148     call ewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
1149     do 120 i = 1,n
1150     if (rwork(i+lewt-1) .le. 0.0d0) go to 621
1151     120 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
1152     c-----------------------------------------------------------------------
1153     c the coding below computes the step size, h0, to be attempted on the
1154     c first step, unless the user has supplied a value for this.
1155     c first check that tout - t differs significantly from zero.
1156     c a scalar tolerance quantity tol is computed, as max(rtol(i))
1157     c if this is positive, or max(atol(i)/abs(y(i))) otherwise, adjusted
1158     c so as to be between 100*uround and 1.0e-3.
1159     c then the computed value h0 is given by..
1160     c neq
1161     c h0**2 = tol / ( w0**-2 + (1/neq) * sum ( f(i)/ywt(i) )**2 )
1162     c 1
1163     c where w0 = max ( abs(t), abs(tout) ),
1164     c f(i) = i-th component of initial value of f,
1165     c ywt(i) = ewt(i)/tol (a weight for y(i)).
1166     c the sign of h0 is inferred from the initial values of tout and t.
1167     c-----------------------------------------------------------------------
1168     if (h0 .ne. 0.0d0) go to 180
1169     tdist = dabs(tout - t)
1170     w0 = dmax1(dabs(t),dabs(tout))
1171     if (tdist .lt. 2.0d0*uround*w0) go to 622
1172     tol = rtol(1)
1173     if (itol .le. 2) go to 140
1174     do 130 i = 1,n
1175     130 tol = dmax1(tol,rtol(i))
1176     140 if (tol .gt. 0.0d0) go to 160
1177     atoli = atol(1)
1178     do 150 i = 1,n
1179     if (itol .eq. 2 .or. itol .eq. 4) atoli = atol(i)
1180     ayi = dabs(y(i))
1181     if (ayi .ne. 0.0d0) tol = dmax1(tol,atoli/ayi)
1182     150 continue
1183     160 tol = dmax1(tol,100.0d0*uround)
1184     tol = dmin1(tol,0.001d0)
1185     sum = vnorm (n, rwork(lf0), rwork(lewt))
1186     sum = 1.0d0/(tol*w0*w0) + tol*sum**2
1187     h0 = 1.0d0/dsqrt(sum)
1188     h0 = dmin1(h0,tdist)
1189     h0 = dsign(h0,tout-t)
1190     c adjust h0 if necessary to meet hmax bound. ---------------------------
1191     180 rh = dabs(h0)*hmxi
1192     if (rh .gt. 1.0d0) h0 = h0/rh
1193     c load h with h0 and scale yh(*,2) by h0. ------------------------------
1194     h = h0
1195     do 190 i = 1,n
1196     190 rwork(i+lf0-1) = h0*rwork(i+lf0-1)
1197     go to 270
1198     c-----------------------------------------------------------------------
1199     c block d.
1200     c the next code block is for continuation calls only (istate = 2 or 3)
1201     c and is to check stop conditions before taking a step.
1202     c-----------------------------------------------------------------------
1203     200 nslast = nst
1204     go to (210, 250, 220, 230, 240), itask
1205     210 if ((tn - tout)*h .lt. 0.0d0) go to 250
1206     call intdy (tout, 0, rwork(lyh), nyh, y, iflag)
1207     if (iflag .ne. 0) go to 627
1208     t = tout
1209     go to 420
1210     220 tp = tn - hu*(1.0d0 + 100.0d0*uround)
1211     if ((tp - tout)*h .gt. 0.0d0) go to 623
1212     if ((tn - tout)*h .lt. 0.0d0) go to 250
1213     go to 400
1214     230 tcrit = rwork(1)
1215     if ((tn - tcrit)*h .gt. 0.0d0) go to 624
1216     if ((tcrit - tout)*h .lt. 0.0d0) go to 625
1217     if ((tn - tout)*h .lt. 0.0d0) go to 245
1218     call intdy (tout, 0, rwork(lyh), nyh, y, iflag)
1219     if (iflag .ne. 0) go to 627
1220     t = tout
1221     go to 420
1222     240 tcrit = rwork(1)
1223     if ((tn - tcrit)*h .gt. 0.0d0) go to 624
1224     245 hmx = dabs(tn) + dabs(h)
1225     ihit = dabs(tn - tcrit) .le. 100.0d0*uround*hmx
1226     if (ihit) go to 400
1227     tnext = tn + h*(1.0d0 + 4.0d0*uround)
1228     if ((tnext - tcrit)*h .le. 0.0d0) go to 250
1229     h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
1230     if (istate .eq. 2) jstart = -2
1231     c-----------------------------------------------------------------------
1232     c block e.
1233     c the next block is normally executed for all calls and contains
1234     c the call to the one-step core integrator stode.
1235     c
1236     c this is a looping point for the integration steps.
1237     c
1238     c first check for too many steps being taken, update ewt (if not at
1239     c start of problem), check for too much accuracy being requested, and
1240     c check for h below the roundoff level in t.
1241     c-----------------------------------------------------------------------
1242     250 continue
1243     if ((nst-nslast) .ge. mxstep) go to 500
1244     call ewset (n, itol, rtol, atol, rwork(lyh), rwork(lewt))
1245     do 260 i = 1,n
1246     if (rwork(i+lewt-1) .le. 0.0d0) go to 510
1247     260 rwork(i+lewt-1) = 1.0d0/rwork(i+lewt-1)
1248     270 tolsf = uround*vnorm (n, rwork(lyh), rwork(lewt))
1249     if (tolsf .le. 1.0d0) go to 280
1250     tolsf = tolsf*2.0d0
1251     if (nst .eq. 0) go to 626
1252     go to 520
1253     280 if ((tn + h) .ne. tn) go to 290
1254     nhnil = nhnil + 1
1255     if (nhnil .gt. mxhnil) go to 290
1256     call xascwv(50hlsode-- warning..internal t (=r1) and h (=r2) are,
1257     1 50, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1258     call xascwv(
1259     1 60h such that in the machine, t + h = t on the next step ,
1260     1 60, 101, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1261     call xascwv(50h (h = step size). solver will continue anyway,
1262     1 50, 101, 0, 0, 0, 0, 2, tn, h)
1263     if (nhnil .lt. mxhnil) go to 290
1264     call xascwv(50hlsode-- above warning has been issued i1 times. ,
1265     1 50, 102, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1266     call xascwv(50h it will not be issued again for this problem,
1267     1 50, 102, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
1268     290 continue
1269     c-----------------------------------------------------------------------
1270     c call stode(neq,y,yh,nyh,yh,ewt,savf,acor,wm,iwm,f,jac,prepj,solsy)
1271     c-----------------------------------------------------------------------
1272     c boyd
1273     call stode (neq, y, rwork(lyh), nyh, rwork(lyh), rwork(lewt),
1274     1 rwork(lsavf), rwork(lacor), rwork(lwm), iwork(liwm),
1275     2 f, jac, prepj, solsy)
1276     kgo = 1 - kflag
1277     go to (300, 530, 540), kgo
1278     c-----------------------------------------------------------------------
1279     c block f.
1280     c the following block handles the case of a successful return from the
1281     c core integrator (kflag = 0). test for stop conditions.
1282     c-----------------------------------------------------------------------
1283     300 init = 1
1284     go to (310, 400, 330, 340, 350), itask
1285     c itask = 1. if tout has been reached, interpolate. -------------------
1286     310 if ((tn - tout)*h .lt. 0.0d0) go to 250
1287     call intdy (tout, 0, rwork(lyh), nyh, y, iflag)
1288     t = tout
1289     go to 420
1290     c itask = 3. jump to exit if tout was reached. ------------------------
1291     330 if ((tn - tout)*h .ge. 0.0d0) go to 400
1292     go to 250
1293     c itask = 4. see if tout or tcrit was reached. adjust h if necessary.
1294     340 if ((tn - tout)*h .lt. 0.0d0) go to 345
1295     call intdy (tout, 0, rwork(lyh), nyh, y, iflag)
1296     t = tout
1297     go to 420
1298     345 hmx = dabs(tn) + dabs(h)
1299     ihit = dabs(tn - tcrit) .le. 100.0d0*uround*hmx
1300     if (ihit) go to 400
1301     tnext = tn + h*(1.0d0 + 4.0d0*uround)
1302     if ((tnext - tcrit)*h .le. 0.0d0) go to 250
1303     h = (tcrit - tn)*(1.0d0 - 4.0d0*uround)
1304     jstart = -2
1305     go to 250
1306     c itask = 5. see if tcrit was reached and jump to exit. ---------------
1307     350 hmx = dabs(tn) + dabs(h)
1308     ihit = dabs(tn - tcrit) .le. 100.0d0*uround*hmx
1309     c-----------------------------------------------------------------------
1310     c block g.
1311     c the following block handles all successful returns from lsode.
1312     c if itask .ne. 1, y is loaded from yh and t is set accordingly.
1313     c istate is set to 2, the illegal input counter is zeroed, and the
1314     c optional outputs are loaded into the work arrays before returning.
1315     c if istate = 1 and tout = t, there is a return with no action taken,
1316     c except that if this has happened repeatedly, the run is terminated.
1317     c-----------------------------------------------------------------------
1318     400 do 410 i = 1,n
1319     410 y(i) = rwork(i+lyh-1)
1320     t = tn
1321     if (itask .ne. 4 .and. itask .ne. 5) go to 420
1322     if (ihit) t = tcrit
1323     420 istate = 2
1324     illin = 0
1325     rwork(11) = hu
1326     rwork(12) = h
1327     rwork(13) = tn
1328     iwork(11) = nst
1329     iwork(12) = nfe
1330     iwork(13) = nje
1331     iwork(14) = nqu
1332     iwork(15) = nq
1333     return
1334     c
1335     430 ntrep = ntrep + 1
1336     if (ntrep .lt. 5) return
1337     call xascwv(
1338     1 60hlsode-- repeated calls with istate = 1 and tout = t (=r1) ,
1339     1 60, 301, 0, 0, 0, 0, 1, t, 0.0d0)
1340     go to 800
1341     c-----------------------------------------------------------------------
1342     c block h.
1343     c the following block handles all unsuccessful returns other than
1344     c those for illegal input. first the error message routine is called.
1345     c if there was an error test or convergence test failure, imxer is set.
1346     c then y is loaded from yh, t is set to tn, and the illegal input
1347     c counter illin is set to 0. the optional outputs are loaded into
1348     c the work arrays before returning.
1349     c-----------------------------------------------------------------------
1350     c the maximum number of steps was taken before reaching tout. ----------
1351     500 call xascwv(50hlsode-- at current t (=r1), mxstep (=i1) steps ,
1352     1 50, 201, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1353     call xascwv(50h taken on this call before reaching tout ,
1354     1 50, 201, 0, 1, mxstep, 0, 1, tn, 0.0d0)
1355     istate = -1
1356     go to 580
1357     c ewt(i) .le. 0.0 for some i (not at start of problem). ----------------
1358     510 ewti = rwork(lewt+i-1)
1359     call xascwv(50hlsode-- at t (=r1), ewt(i1) has become r2 .le. 0.,
1360     1 50, 202, 0, 1, i, 0, 2, tn, ewti)
1361     istate = -6
1362     go to 580
1363     c too much accuracy requested for machine precision. -------------------
1364     520 call xascwv(50hlsode-- at t (=r1), too much accuracy requested ,
1365     1 50, 203, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1366     call xascwv(50h for precision of machine.. see tolsf (=r2) ,
1367     1 50, 203, 0, 0, 0, 0, 2, tn, tolsf)
1368     rwork(14) = tolsf
1369     istate = -2
1370     go to 580
1371     c kflag = -1. error test failed repeatedly or with abs(h) = hmin. -----
1372     530 call xascwv(50hlsode-- at t(=r1) and step size h(=r2), the error,
1373     1 50, 204, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1374     call xascwv(50h test failed repeatedly or with abs(h) = hmin,
1375     1 50, 204, 0, 0, 0, 0, 2, tn, h)
1376     istate = -4
1377     go to 560
1378     c kflag = -2. convergence failed repeatedly or with abs(h) = hmin. ----
1379     540 call xascwv(50hlsode-- at t (=r1) and step size h (=r2), the ,
1380     1 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1381     call xascwv(50h corrector convergence failed repeatedly ,
1382     1 50, 205, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1383     call xascwv(30h or with abs(h) = hmin ,
1384     1 30, 205, 0, 0, 0, 0, 2, tn, h)
1385     istate = -5
1386     c compute imxer if relevant. -------------------------------------------
1387     560 big = 0.0d0
1388     imxer = 1
1389     do 570 i = 1,n
1390     size = dabs(rwork(i+lacor-1)*rwork(i+lewt-1))
1391     if (big .ge. size) go to 570
1392     big = size
1393     imxer = i
1394     570 continue
1395     iwork(16) = imxer
1396     c set y vector, t, illin, and optional outputs. ------------------------
1397     580 do 590 i = 1,n
1398     590 y(i) = rwork(i+lyh-1)
1399     t = tn
1400     illin = 0
1401     rwork(11) = hu
1402     rwork(12) = h
1403     rwork(13) = tn
1404     iwork(11) = nst
1405     iwork(12) = nfe
1406     iwork(13) = nje
1407     iwork(14) = nqu
1408     iwork(15) = nq
1409     return
1410     c-----------------------------------------------------------------------
1411     c block i.
1412     c the following block handles all error returns due to illegal input
1413     c (istate = -3), as detected before calling the core integrator.
1414     c first the error message routine is called. then if there have been
1415     c 5 consecutive such returns just before this call to the solver,
1416     c the run is halted.
1417     c-----------------------------------------------------------------------
1418     601 call xascwv(30hlsode-- istate (=i1) illegal ,
1419     1 30, 1, 0, 1, istate, 0, 0, 0.0d0, 0.0d0)
1420     go to 700
1421     602 call xascwv(30hlsode-- itask (=i1) illegal ,
1422     1 30, 2, 0, 1, itask, 0, 0, 0.0d0, 0.0d0)
1423     go to 700
1424     603 call xascwv(50hlsode-- istate .gt. 1 but lsode not initialized ,
1425     1 50, 3, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1426     go to 700
1427     604 call xascwv(30hlsode-- neq (=i1) .lt. 1 ,
1428     1 30, 4, 0, 1, neq(1), 0, 0, 0.0d0, 0.0d0)
1429     go to 700
1430     605 call xascwv(50hlsode-- istate = 3 and neq increased (i1 to i2) ,
1431     1 50, 5, 0, 2, n, neq(1), 0, 0.0d0, 0.0d0)
1432     go to 700
1433     606 call xascwv(30hlsode-- itol (=i1) illegal ,
1434     1 30, 6, 0, 1, itol, 0, 0, 0.0d0, 0.0d0)
1435     go to 700
1436     607 call xascwv(30hlsode-- iopt (=i1) illegal ,
1437     1 30, 7, 0, 1, iopt, 0, 0, 0.0d0, 0.0d0)
1438     go to 700
1439     608 call xascwv(30hlsode-- mf (=i1) illegal ,
1440     1 30, 8, 0, 1, mf, 0, 0, 0.0d0, 0.0d0)
1441     go to 700
1442     609 call xascwv(50hlsode-- ml (=i1) illegal.. .lt.0 or .ge.neq (=i2),
1443     1 50, 9, 0, 2, ml, neq(1), 0, 0.0d0, 0.0d0)
1444     go to 700
1445     610 call xascwv(50hlsode-- mu (=i1) illegal.. .lt.0 or .ge.neq (=i2),
1446     1 50, 10, 0, 2, mu, neq(1), 0, 0.0d0, 0.0d0)
1447     go to 700
1448     611 call xascwv(30hlsode-- maxord (=i1) .lt. 0 ,
1449     1 30, 11, 0, 1, maxord, 0, 0, 0.0d0, 0.0d0)
1450     go to 700
1451     612 call xascwv(30hlsode-- mxstep (=i1) .lt. 0 ,
1452     1 30, 12, 0, 1, mxstep, 0, 0, 0.0d0, 0.0d0)
1453     go to 700
1454     613 call xascwv(30hlsode-- mxhnil (=i1) .lt. 0 ,
1455     1 30, 13, 0, 1, mxhnil, 0, 0, 0.0d0, 0.0d0)
1456     go to 700
1457     614 call xascwv(40hlsode-- tout (=r1) behind t (=r2) ,
1458     1 40, 14, 0, 0, 0, 0, 2, tout, t)
1459     call xascwv(50h integration direction is given by h0 (=r1) ,
1460     1 50, 14, 0, 0, 0, 0, 1, h0, 0.0d0)
1461     go to 700
1462     615 call xascwv(30hlsode-- hmax (=r1) .lt. 0.0 ,
1463     1 30, 15, 0, 0, 0, 0, 1, hmax, 0.0d0)
1464     go to 700
1465     616 call xascwv(30hlsode-- hmin (=r1) .lt. 0.0 ,
1466     1 30, 16, 0, 0, 0, 0, 1, hmin, 0.0d0)
1467     go to 700
1468     617 call xascwv(
1469     1 60hlsode-- rwork length needed, lenrw (=i1), exceeds lrw (=i2),
1470     1 60, 17, 0, 2, lenrw, lrw, 0, 0.0d0, 0.0d0)
1471     go to 700
1472     618 call xascwv(
1473     1 60hlsode-- iwork length needed, leniw (=i1), exceeds liw (=i2),
1474     1 60, 18, 0, 2, leniw, liw, 0, 0.0d0, 0.0d0)
1475     go to 700
1476     619 call xascwv(40hlsode-- rtol(i1) is r1 .lt. 0.0 ,
1477     1 40, 19, 0, 1, i, 0, 1, rtoli, 0.0d0)
1478     go to 700
1479     620 call xascwv(40hlsode-- atol(i1) is r1 .lt. 0.0 ,
1480     1 40, 20, 0, 1, i, 0, 1, atoli, 0.0d0)
1481     go to 700
1482     621 ewti = rwork(lewt+i-1)
1483     call xascwv(40hlsode-- ewt(i1) is r1 .le. 0.0 ,
1484     1 40, 21, 0, 1, i, 0, 1, ewti, 0.0d0)
1485     go to 700
1486     622 call xascwv(
1487     1 60hlsode-- tout (=r1) too close to t(=r2) to start integration,
1488     1 60, 22, 0, 0, 0, 0, 2, tout, t)
1489     go to 700
1490     623 call xascwv(
1491     1 60hlsode-- itask = i1 and tout (=r1) behind tcur - hu (= r2) ,
1492     1 60, 23, 0, 1, itask, 0, 2, tout, tp)
1493     go to 700
1494     624 call xascwv(
1495     1 60hlsode-- itask = 4 or 5 and tcrit (=r1) behind tcur (=r2) ,
1496     1 60, 24, 0, 0, 0, 0, 2, tcrit, tn)
1497     go to 700
1498     625 call xascwv(
1499     1 60hlsode-- itask = 4 or 5 and tcrit (=r1) behind tout (=r2) ,
1500     1 60, 25, 0, 0, 0, 0, 2, tcrit, tout)
1501     go to 700
1502     626 call xascwv(50hlsode-- at start of problem, too much accuracy ,
1503     1 50, 26, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1504     call xascwv(
1505     1 60h requested for precision of machine.. see tolsf (=r1) ,
1506     1 60, 26, 0, 0, 0, 0, 1, tolsf, 0.0d0)
1507     rwork(14) = tolsf
1508     go to 700
1509     627 call xascwv(50hlsode-- trouble from intdy. itask = i1, tout = r1,
1510     1 50, 27, 0, 1, itask, 0, 1, tout, 0.0d0)
1511     c
1512     700 if (illin .eq. 5) go to 710
1513     illin = illin + 1
1514     istate = -3
1515     return
1516     710 call xascwv(50hlsode-- repeated occurrences of illegal input ,
1517     1 50, 302, 0, 0, 0, 0, 0, 0.0d0, 0.0d0)
1518     c
1519     800 call xascwv(50hlsode-- run aborted.. apparent infinite loop ,
1520     1 50, 303, 2, 0, 0, 0, 0, 0.0d0, 0.0d0)
1521     return
1522     c----------------------- end of subroutine lsode -----------------------
1523     end
1524     subroutine solsy (wm, iwm, x, tem)
1525     clll. optimize
1526     integer iwm
1527     integer iownd, iowns,
1528     1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
1529     2 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1530     integer i, meband, ml, mu
1531     double precision wm, x, tem
1532     double precision rowns,
1533     1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
1534     double precision di, hl0, phl0, r
1535     dimension wm(1), iwm(1), x(1), tem(1)
1536     common /ls0001/ rowns(209),
1537     2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
1538     3 iownd(14), iowns(6),
1539     4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
1540     5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1541     c-----------------------------------------------------------------------
1542     c this routine manages the solution of the linear system arising from
1543     c a chord iteration. it is called if miter .ne. 0.
1544     c if miter is 1 or 2, it calls dgesl to accomplish this.
1545     c if miter = 3 it updates the coefficient h*el0 in the diagonal
1546     c matrix, and then computes the solution.
1547     c if miter is 4 or 5, it calls dgbsl.
1548     c communication with solsy uses the following variables..
1549     c wm = real work space containing the inverse diagonal matrix if
1550     c miter = 3 and the lu decomposition of the matrix otherwise.
1551     c storage of matrix elements starts at wm(3).
1552     c wm also contains the following matrix-related data..
1553     c wm(1) = sqrt(uround) (not used here),
1554     c wm(2) = hl0, the previous value of h*el0, used if miter = 3.
1555     c iwm = integer work space containing pivot information, starting at
1556     c iwm(21), if miter is 1, 2, 4, or 5. iwm also contains band
1557     c parameters ml = iwm(1) and mu = iwm(2) if miter is 4 or 5.
1558     c x = the right-hand side vector on input, and the solution vector
1559     c on output, of length n.
1560     c tem = vector of work space of length n, not used in this version.
1561     c iersl = output flag (in common). iersl = 0 if no trouble occurred.
1562     c iersl = 1 if a singular matrix arose with miter = 3.
1563     c this routine also uses the common variables el0, h, miter, and n.
1564     c-----------------------------------------------------------------------
1565     iersl = 0
1566     go to (100, 100, 300, 400, 400), miter
1567     100 call dgesl (wm(3), n, n, iwm(21), x, 0)
1568     return
1569     c
1570     300 phl0 = wm(2)
1571     hl0 = h*el0
1572     wm(2) = hl0
1573     if (hl0 .eq. phl0) go to 330
1574     r = hl0/phl0
1575     do 320 i = 1,n
1576     di = 1.0d0 - r*(1.0d0 - 1.0d0/wm(i+2))
1577     if (dabs(di) .eq. 0.0d0) go to 390
1578     320 wm(i+2) = 1.0d0/di
1579     330 do 340 i = 1,n
1580     340 x(i) = wm(i+2)*x(i)
1581     return
1582     390 iersl = 1
1583     return
1584     c
1585     400 ml = iwm(1)
1586     mu = iwm(2)
1587     meband = 2*ml + mu + 1
1588     call dgbsl (wm(3), meband, n, ml, mu, iwm(21), x, 0)
1589     return
1590     c----------------------- end of subroutine solsy -----------------------
1591     end
1592     subroutine prepj (neq, y, yh, nyh, ewt, ftem, savf, wm, iwm,
1593     1 f, jac)
1594     clll. optimize
1595     cascend changes
1596     external f, jac
1597     integer neq, nyh, iwm
1598     integer iownd, iowns,
1599     1 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
1600     2 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1601     integer i, i1, i2, ier, ii, j, j1, jj, lenp,
1602     1 mba, mband, meb1, meband, ml, ml3, mu, np1
1603     double precision y, yh, ewt, ftem, savf, wm
1604     double precision rowns,
1605     1 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround
1606     double precision con, di, fac, hl0, r, r0, srur, yi, yj, yjj,
1607     1 vnorm
1608     Ckaa double precision time1, time2
1609     dimension neq(1), y(1), yh(nyh,1), ewt(1), ftem(1), savf(1),
1610     1 wm(1), iwm(1)
1611     common /ls0001/ rowns(209),
1612     2 ccmax, el0, h, hmin, hmxi, hu, rc, tn, uround,
1613     3 iownd(14), iowns(6),
1614     4 icf, ierpj, iersl, jcur, jstart, kflag, l, meth, miter,
1615     5 maxord, maxcor, msbp, mxncf, n, nq, nst, nfe, nje, nqu
1616     c-----------------------------------------------------------------------
1617     c prepj is called by stode to compute and process the matrix
1618     c p = i - h*el(1)*j , where j is an approximation to the jacobian.
1619     c here j is computed by the user-supplied routine jac if
1620     c miter = 1 or 4, or by finite differencing if miter = 2, 3, or 5.
1621     c if miter = 3, a diagonal approximation to j is used.
1622     c j is stored in wm and replaced by p. if miter .ne. 3, p is then
1623     c subjected to lu decomposition in preparation for later solution
1624     c of linear systems with p as coefficient matrix. this is done
1625     c by dgefa if miter = 1 or 2, and by dgbfa if miter = 4 or 5.
1626     c
1627     c in addition to variables described previously, communication
1628     c with prepj uses the following..
1629     c y = array containing predicted values on entry.
1630     c ftem = work array of length n (acor in stode).
1631     c savf = array containing f evaluated at predicted y.
1632     c wm = real work space for matrices. on output it contains the
1633     c inverse diagonal matrix if miter = 3 and the lu decomposition
1634     c of p if miter is 1, 2 , 4, or 5.
1635     c storage of matrix elements starts at wm(3).
1636     c wm also contains the following matrix-related data..
1637     c wm(1) = sqrt(uround), used in numerical jacobian increments.
1638     c wm(2) = h*el0, saved for later use if miter = 3.
1639     c iwm = integer work space containing pivot information, starting at
1640     c iwm(21), if miter is 1, 2, 4, or 5. iwm also contains band
1641     c parameters ml = iwm(1) and mu = iwm(2) if miter is 4 or 5.
1642     c el0 = el(1) (input).
1643     c ierpj = output error flag, = 0 if no trouble, .gt. 0 if
1644     c p matrix found to be singular.
1645     c jcur = output flag = 1 to indicate that the jacobian matrix
1646     c (or approximation) is now current.
1647     c this routine also uses the common variables el0, h, tn, uround,
1648     c miter, n, nfe, and nje.
1649     c-----------------------------------------------------------------------
1650     nje = nje + 1
1651     ierpj = 0
1652     jcur = 1
1653     hl0 = h*el0