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