| blob | d7b0e39363f1a2e5d13428b71edb198d333b9951 |
1 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
2 ! LSODE - Stiff method based on backward differentiation formulas (BDF) !
3 ! By default the code employs the KPP sparse linear algebra routines !
4 ! Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK) !
5 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
6 ! A. Sandu - version of July 2005
8 MODULE KPP_ROOT_Integrator
10 USE KPP_ROOT_Precision
15 Set2zero, WLAMCH
17 IMPLICIT NONE
18 PUBLIC
19 SAVE
21 !~~~> Statistics on the work performed by the LSODE method
25 ! SDIRK method coefficients
29 ! mz_rs_20050717: TODO: use strings of IERR_NAMES for error messages
30 ! description of the error numbers IERR
43 CONTAINS
48 USE KPP_ROOT_Parameters
49 USE KPP_ROOT_Global
50 IMPLICIT NONE
54 ! Optional input parameters and statistics
63 !!$ INTEGER, SAVE :: Ntotal = 0
72 ! If optional parameters are given, and if they are >0,
73 ! then they overwrite default settings.
76 END IF
79 END IF
84 ! mz_rs_20050716: IERR and ISTATUS are returned to the user who then
85 ! decides what to do about it, i.e. either stop the run or ignore it.
86 !!$ IF (IERR < 0) THEN
87 !!$ PRINT*,'LSODE: Unsuccessful exit at T=',TIN,' (IERR=',IERR,')'
88 !!$ STOP
89 !!$ ENDIF
90 !!$ Ntotal = Ntotal + ISTATUS(3)
91 !!$ PRINT*,'Nsteps = ', ISTATUS(3),' (',Ntotal,')'
95 ! if optional parameters are given for output they to return information
100 END SUBROUTINE INTEGRATE
102 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
105 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
106 !
107 !~~~> INPUT PARAMETERS:
108 !
109 ! Note: For input parameters equal to zero the default values of the
110 ! corresponding variables are used.
111 !
112 ! Note: For input parameters equal to zero the default values of the
113 ! corresponding variables are used.
114 !~~~>
115 ! ICNTRL(1) = not used
116 !
117 ! ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors
118 ! = 1: AbsTol, RelTol are scalars
119 !
120 ! ICNTRL(3) = not used
121 !
122 ! ICNTRL(4) -> maximum number of integration steps
123 ! For ICNTRL(4)=0 the default value of 100000 is used
124 !
125 ! ICNTRL(5) -> maximum order of the integration formula allowed
126 !
127 !~~~> Real parameters
128 !
129 ! RCNTRL(1) -> Hmin, lower bound for the integration step size
130 ! It is strongly recommended to keep Hmin = ZERO
131 ! RCNTRL(2) -> Hmax, upper bound for the integration step size
132 ! RCNTRL(3) -> Hstart, starting value for the integration step size
133 !
134 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
135 !
136 !~~~> OUTPUT PARAMETERS:
137 !
138 ! Note: each call to Rosenbrock adds the current no. of fcn calls
139 ! to previous value of ISTATUS(1), and similar for the other params.
140 ! Set ISTATUS(1:10) = 0 before call to avoid this accumulation.
141 !
142 ! ISTATUS(1) = No. of function calls
143 ! ISTATUS(2) = No. of jacobian calls
144 ! ISTATUS(3) = No. of steps
145 !
146 ! RSTATUS(1) -> Texit, the time corresponding to the
147 ! computed Y upon return
148 ! RSTATUS(2) -> Hexit, last predicted step before exit
149 ! For multiple restarts, use Hexit as Hstart in the following run
150 !
151 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
153 IMPLICIT NONE
163 !~~~> NORMAL COMPUTATION
173 ELSE
175 END IF
195 END SUBROUTINE KppLsode
196 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
199 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
200 !DECK DLSODE
206 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
207 !***BEGIN PROLOGUE DLSODE
208 !***PURPOSE Livermore Solver for Ordinary Differential Equations.
209 ! DLSODE solves the initial-value problem for stiff or
210 ! nonstiff systems of first-order ODE's,
211 ! dy/dt = f(t,y), or, in component form,
212 ! dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(N)), i=1,...,N.
213 !***CATEGORY I1A
214 !***TYPE KPP_REAL (SLSODE-S, DLSODE-D)
215 !***KEYWORDS ORDINARY DIFFERENTIAL EQUATIONS, INITIAL VALUE PROBLEM,
216 ! STIFF, NONSTIFF
217 !***AUTHOR Hindmarsh, Alan C., (LLNL)
218 ! Center for Applied Scientific Computing, L-561
219 ! Lawrence Livermore National Laboratory
220 ! Livermore, CA 94551.
221 !***DESCRIPTION
222 !
223 ! NOTE: The "Usage" and "Arguments" sections treat only a subset of
224 ! available options, in condensed fashion. The options
225 ! covered and the information supplied will support most
226 ! standard uses of DLSODE.
227 !
228 ! For more sophisticated uses, full details on all options are
229 ! given in the concluding section, headed "Long Description."
230 ! A synopsis of the DLSODE Long Description is provided at the
231 ! beginning of that section; general topics covered are:
232 ! - Elements of the call sequence; optional input and output
233 ! - Optional supplemental routines in the DLSODE package
234 ! - internal COMMON block
235 !
236 ! *Usage:
237 ! Communication between the user and the DLSODE package, for normal
238 ! situations, is summarized here. This summary describes a subset
239 ! of the available options. See "Long Description" for complete
240 ! details, including optional communication, nonstandard options,
241 ! and instructions for special situations.
242 !
243 ! A sample program is given in the "Examples" section.
244 !
245 ! Refer to the argument descriptions for the definitions of the
246 ! quantities that appear in the following sample declarations.
247 !
248 ! For MF = 10,
249 ! PARAMETER (LRW = 20 + 16*NEQ, LIW = 20)
250 ! For MF = 21 or 22,
251 ! PARAMETER (LRW = 22 + 9*NEQ + NEQ**2, LIW = 20 + NEQ)
252 ! For MF = 24 or 25,
253 ! PARAMETER (LRW = 22 + 10*NEQ + (2*ML+MU)*NEQ,
254 ! * LIW = 20 + NEQ)
255 !
256 ! EXTERNAL F, JAC
257 ! INTEGER NEQ, ITOL, ITASK, ISTATE, IOPT, LRW, IWORK(LIW),
258 ! * LIW, MF
259 ! KPP_REAL Y(NEQ), T, TOUT, RelTol, AbsTol(ntol), RWORK(LRW)
260 !
261 ! CALL DLSODE (F, NEQ, Y, T, TOUT, ITOL, RelTol, AbsTol, ITASK,
262 ! * ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF)
263 !
264 ! *Arguments:
265 ! F :EXT Name of subroutine for right-hand-side vector f.
266 ! This name must be declared EXTERNAL in calling
267 ! program. The form of F must be:
268 !
269 ! SUBROUTINE F (NEQ, T, Y, YDOT)
270 ! INTEGER NEQ
271 ! KPP_REAL T, Y(*), YDOT(*)
272 !
273 ! The inputs are NEQ, T, Y. F is to set
274 !
275 ! YDOT(i) = f(i,T,Y(1),Y(2),...,Y(NEQ)),
276 ! i = 1, ..., NEQ .
277 !
278 ! NEQ :IN Number of first-order ODE's.
279 !
280 ! Y :INOUT Array of values of the y(t) vector, of length NEQ.
281 ! Input: For the first call, Y should contain the
282 ! values of y(t) at t = T. (Y is an input
283 ! variable only if ISTATE = 1.)
284 ! Output: On return, Y will contain the values at the
285 ! new t-value.
286 !
287 ! T :INOUT Value of the independent variable. On return it
288 ! will be the current value of t (normally TOUT).
289 !
290 ! TOUT :IN Next point where output is desired (.NE. T).
291 !
292 ! ITOL :IN 1 or 2 according as AbsTol (below) is a scalar or
293 ! an array.
294 !
295 ! RelTol :IN Relative tolerance parameter (scalar).
296 !
297 ! AbsTol :IN Absolute tolerance parameter (scalar or array).
298 ! If ITOL = 1, AbsTol need not be dimensioned.
299 ! If ITOL = 2, AbsTol must be dimensioned at least NEQ.
300 !
301 ! The estimated local error in Y(i) will be controlled
302 ! so as to be roughly less (in magnitude) than
303 !
304 ! EWT(i) = RelTol*ABS(Y(i)) + AbsTol if ITOL = 1, or
305 ! EWT(i) = RelTol*ABS(Y(i)) + AbsTol(i) if ITOL = 2.
306 !
307 ! Thus the local error test passes if, in each
308 ! component, either the absolute error is less than
309 ! AbsTol (or AbsTol(i)), or the relative error is less
310 ! than RelTol.
311 !
312 ! Use RelTol = 0.0 for pure absolute error control, and
313 ! use AbsTol = 0.0 (or AbsTol(i) = 0.0) for pure relative
314 ! error control. Caution: Actual (global) errors may
315 ! exceed these local tolerances, so choose them
316 ! conservatively.
317 !
318 ! ITASK :IN Flag indicating the task DLSODE is to perform.
319 ! Use ITASK = 1 for normal computation of output
320 ! values of y at t = TOUT.
321 !
322 ! ISTATE:INOUT Index used for input and output to specify the state
323 ! of the calculation.
324 ! Input:
325 ! 1 This is the first call for a problem.
326 ! 2 This is a subsequent call.
327 ! Output:
328 ! 1 Nothing was done, because TOUT was equal to T.
329 ! 2 DLSODE was successful (otherwise, negative).
330 ! Note that ISTATE need not be modified after a
331 ! successful return.
332 ! -1 Excess work done on this call (perhaps wrong
333 ! MF).
334 ! -2 Excess accuracy requested (tolerances too
335 ! small).
336 ! -3 Illegal input detected (see printed message).
337 ! -4 Repeated error test failures (check all
338 ! inputs).
339 ! -5 Repeated convergence failures (perhaps bad
340 ! Jacobian supplied or wrong choice of MF or
341 ! tolerances).
342 ! -6 Error weight became zero during problem
343 ! (solution component i vanished, and AbsTol or
344 ! AbsTol(i) = 0.).
345 !
346 ! IOPT :IN Flag indicating whether optional inputs are used:
347 ! 0 No.
348 ! 1 Yes. (See "Optional inputs" under "Long
349 ! Description," Part 1.)
350 !
351 ! RWORK :WORK Real work array of length at least:
352 ! 20 + 16*NEQ for MF = 10,
353 ! 22 + 9*NEQ + NEQ**2 for MF = 21 or 22,
354 ! 22 + 10*NEQ + (2*ML + MU)*NEQ for MF = 24 or 25.
355 !
356 ! LRW :IN Declared length of RWORK (in user's DIMENSION
357 ! statement).
358 !
359 ! IWORK :WORK Integer work array of length at least:
360 ! 20 for MF = 10,
361 ! 20 + NEQ for MF = 21, 22, 24, or 25.
362 !
363 ! If MF = 24 or 25, input in IWORK(1),IWORK(2) the
364 ! lower and upper Jacobian half-bandwidths ML,MU.
365 !
366 ! On return, IWORK contains information that may be
367 ! of interest to the user:
368 !
369 ! Name Location Meaning
370 ! ----- --------- -----------------------------------------
371 ! NST IWORK(11) Number of steps taken for the problem so
372 ! far.
373 ! NFE IWORK(12) Number of f evaluations for the problem
374 ! so far.
375 ! NJE IWORK(13) Number of Jacobian evaluations (and of
376 ! matrix LU decompositions) for the problem
377 ! so far.
378 ! NQU IWORK(14) Method order last used (successfully).
379 ! LENRW IWORK(17) Length of RWORK actually required. This
380 ! is defined on normal returns and on an
381 ! illegal input return for insufficient
382 ! storage.
383 ! LENIW IWORK(18) Length of IWORK actually required. This
384 ! is defined on normal returns and on an
385 ! illegal input return for insufficient
386 ! storage.
387 !
388 ! LIW :IN Declared length of IWORK (in user's DIMENSION
389 ! statement).
390 !
391 ! JAC :EXT Name of subroutine for Jacobian matrix (MF =
392 ! 21 or 24). If used, this name must be declared
393 ! EXTERNAL in calling program. If not used, pass a
394 ! dummy name. The form of JAC must be:
395 !
396 ! SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD)
397 ! INTEGER NEQ, ML, MU, NROWPD
398 ! KPP_REAL T, Y(*), PD(NROWPD,*)
399 !
400 ! See item c, under "Description" below for more
401 ! information about JAC.
402 !
403 ! MF :IN Method flag. Standard values are:
404 ! 10 Nonstiff (Adams) method, no Jacobian used.
405 ! 21 Stiff (BDF) method, user-supplied full Jacobian.
406 ! 22 Stiff method, internally generated full
407 ! Jacobian.
408 ! 24 Stiff method, user-supplied banded Jacobian.
409 ! 25 Stiff method, internally generated banded
410 ! Jacobian.
411 !
412 ! *Description:
413 ! DLSODE solves the initial value problem for stiff or nonstiff
414 ! systems of first-order ODE's,
415 !
416 ! dy/dt = f(t,y) ,
417 !
418 ! or, in component form,
419 !
420 ! dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(NEQ))
421 ! (i = 1, ..., NEQ) .
422 !
423 ! DLSODE is a package based on the GEAR and GEARB packages, and on
424 ! the October 23, 1978, version of the tentative ODEPACK user
425 ! interface standard, with minor modifications.
426 !
427 ! The steps in solving such a problem are as follows.
428 !
429 ! a. First write a subroutine of the form
430 !
431 ! SUBROUTINE F (NEQ, T, Y, YDOT)
432 ! INTEGER NEQ
433 ! KPP_REAL T, Y(*), YDOT(*)
434 !
435 ! which supplies the vector function f by loading YDOT(i) with
436 ! f(i).
437 !
438 ! b. Next determine (or guess) whether or not the problem is stiff.
439 ! Stiffness occurs when the Jacobian matrix df/dy has an
440 ! eigenvalue whose real part is negative and large in magnitude
441 ! compared to the reciprocal of the t span of interest. If the
442 ! problem is nonstiff, use method flag MF = 10. If it is stiff,
443 ! there are four standard choices for MF, and DLSODE requires the
444 ! Jacobian matrix in some form. This matrix is regarded either
445 ! as full (MF = 21 or 22), or banded (MF = 24 or 25). In the
446 ! banded case, DLSODE requires two half-bandwidth parameters ML
447 ! and MU. These are, respectively, the widths of the lower and
448 ! upper parts of the band, excluding the main diagonal. Thus the
449 ! band consists of the locations (i,j) with
450 !
451 ! i - ML <= j <= i + MU ,
452 !
453 ! and the full bandwidth is ML + MU + 1 .
454 !
455 ! c. If the problem is stiff, you are encouraged to supply the
456 ! Jacobian directly (MF = 21 or 24), but if this is not feasible,
457 ! DLSODE will compute it internally by difference quotients (MF =
458 ! 22 or 25). If you are supplying the Jacobian, write a
459 ! subroutine of the form
460 !
461 ! SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD)
462 ! INTEGER NEQ, ML, MU, NRWOPD
463 ! KPP_REAL T, Y(*), PD(NROWPD,*)
464 !
465 ! which provides df/dy by loading PD as follows:
466 ! - For a full Jacobian (MF = 21), load PD(i,j) with df(i)/dy(j),
467 ! the partial derivative of f(i) with respect to y(j). (Ignore
468 ! the ML and MU arguments in this case.)
469 ! - For a banded Jacobian (MF = 24), load PD(i-j+MU+1,j) with
470 ! df(i)/dy(j); i.e., load the diagonal lines of df/dy into the
471 ! rows of PD from the top down.
472 ! - In either case, only nonzero elements need be loaded.
473 !
474 ! d. Write a main program that calls subroutine DLSODE once for each
475 ! point at which answers are desired. This should also provide
476 ! for possible use of logical unit 6 for output of error messages
477 ! by DLSODE.
478 !
479 ! Before the first call to DLSODE, set ISTATE = 1, set Y and T to
480 ! the initial values, and set TOUT to the first output point. To
481 ! continue the integration after a successful return, simply
482 ! reset TOUT and call DLSODE again. No other parameters need be
483 ! reset.
484 !
485 ! *Examples:
486 ! The following is a simple example problem, with the coding needed
487 ! for its solution by DLSODE. The problem is from chemical kinetics,
488 ! and consists of the following three rate equations:
489 !
490 ! dy1/dt = -.04*y1 + 1.E4*y2*y3
491 ! dy2/dt = .04*y1 - 1.E4*y2*y3 - 3.E7*y2**2
492 ! dy3/dt = 3.E7*y2**2
493 !
494 ! on the interval from t = 0.0 to t = 4.E10, with initial conditions
495 ! y1 = 1.0, y2 = y3 = 0. The problem is stiff.
496 !
497 ! The following coding solves this problem with DLSODE, using
498 ! MF = 21 and printing results at t = .4, 4., ..., 4.E10. It uses
499 ! ITOL = 2 and AbsTol much smaller for y2 than for y1 or y3 because y2
500 ! has much smaller values. At the end of the run, statistical
501 ! quantities of interest are printed.
502 !
503 ! EXTERNAL FEX, JEX
504 ! INTEGER IOPT, IOUT, ISTATE, ITASK, ITOL, IWORK(23), LIW, LRW,
505 ! * MF, NEQ
506 ! KPP_REAL AbsTol(3), RelTol, RWORK(58), T, TOUT, Y(3)
507 ! NEQ = 3
508 ! Y(1) = 1.D0
509 ! Y(2) = 0.D0
510 ! Y(3) = 0.D0
511 ! T = 0.D0
512 ! TOUT = .4D0
513 ! ITOL = 2
514 ! RelTol = 1.D-4
515 ! AbsTol(1) = 1.D-6
516 ! AbsTol(2) = 1.D-10
517 ! AbsTol(3) = 1.D-6
518 ! ITASK = 1
519 ! ISTATE = 1
520 ! IOPT = 0
521 ! LRW = 58
522 ! LIW = 23
523 ! MF = 21
524 ! DO 40 IOUT = 1,12
525 ! CALL DLSODE (FEX, NEQ, Y, T, TOUT, ITOL, RelTol, AbsTol, ITASK,
526 ! * ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JEX, MF)
527 ! WRITE(6,20) T, Y(1), Y(2), Y(3)
528 ! 20 FORMAT(' At t =',D12.4,' y =',3D14.6)
529 ! IF (ISTATE .LT. 0) GO TO 80
530 ! 40 TOUT = TOUT*10.D0
531 ! WRITE(6,60) IWORK(11), IWORK(12), IWORK(13)
532 ! 60 FORMAT(/' No. steps =',i4,', No. f-s =',i4,', No. J-s =',i4)
533 ! STOP
534 ! 80 WRITE(6,90) ISTATE
535 ! 90 FORMAT(///' Error halt.. ISTATE =',I3)
536 ! STOP
537 ! END
538 !
539 ! SUBROUTINE FEX (NEQ, T, Y, YDOT)
540 ! INTEGER NEQ
541 ! KPP_REAL T, Y(3), YDOT(3)
542 ! YDOT(1) = -.04D0*Y(1) + 1.D4*Y(2)*Y(3)
543 ! YDOT(3) = 3.D7*Y(2)*Y(2)
544 ! YDOT(2) = -YDOT(1) - YDOT(3)
545 ! RETURN
546 ! END
547 !
548 ! SUBROUTINE JEX (NEQ, T, Y, ML, MU, PD, NRPD)
549 ! INTEGER NEQ, ML, MU, NRPD
550 ! KPP_REAL T, Y(3), PD(NRPD,3)
551 ! PD(1,1) = -.04D0
552 ! PD(1,2) = 1.D4*Y(3)
553 ! PD(1,3) = 1.D4*Y(2)
554 ! PD(2,1) = .04D0
555 ! PD(2,3) = -PD(1,3)
556 ! PD(3,2) = 6.D7*Y(2)
557 ! PD(2,2) = -PD(1,2) - PD(3,2)
558 ! RETURN
559 ! END
560 !
561 ! The output from this program (on a Cray-1 in single precision)
562 ! is as follows.
563 !
564 ! At t = 4.0000e-01 y = 9.851726e-01 3.386406e-05 1.479357e-02
565 ! At t = 4.0000e+00 y = 9.055142e-01 2.240418e-05 9.446344e-02
566 ! At t = 4.0000e+01 y = 7.158050e-01 9.184616e-06 2.841858e-01
567 ! At t = 4.0000e+02 y = 4.504846e-01 3.222434e-06 5.495122e-01
568 ! At t = 4.0000e+03 y = 1.831701e-01 8.940379e-07 8.168290e-01
569 ! At t = 4.0000e+04 y = 3.897016e-02 1.621193e-07 9.610297e-01
570 ! At t = 4.0000e+05 y = 4.935213e-03 1.983756e-08 9.950648e-01
571 ! At t = 4.0000e+06 y = 5.159269e-04 2.064759e-09 9.994841e-01
572 ! At t = 4.0000e+07 y = 5.306413e-05 2.122677e-10 9.999469e-01
573 ! At t = 4.0000e+08 y = 5.494530e-06 2.197825e-11 9.999945e-01
574 ! At t = 4.0000e+09 y = 5.129458e-07 2.051784e-12 9.999995e-01
575 ! At t = 4.0000e+10 y = -7.170603e-08 -2.868241e-13 1.000000e+00
576 !
577 ! No. steps = 330, No. f-s = 405, No. J-s = 69
578 !
579 ! *Accuracy:
580 ! The accuracy of the solution depends on the choice of tolerances
581 ! RelTol and AbsTol. Actual (global) errors may exceed these local
582 ! tolerances, so choose them conservatively.
583 !
584 ! *Cautions:
585 ! The work arrays should not be altered between calls to DLSODE for
586 ! the same problem, except possibly for the conditional and optional
587 ! inputs.
588 !
589 ! *Portability:
590 ! Since NEQ is dimensioned inside DLSODE, some compilers may object
591 ! to a call to DLSODE with NEQ a scalar variable. In this event,
592 ! use DIMENSION NEQ. Similar remarks apply to RelTol and AbsTol.
593 !
594 ! Note to Cray users:
595 ! For maximum efficiency, use the CFT77 compiler. Appropriate
596 ! compiler optimization directives have been inserted for CFT77.
597 !
598 ! *Reference:
599 ! Alan C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE
600 ! Solvers," in Scientific Computing, R. S. Stepleman, et al., Eds.
601 ! (North-Holland, Amsterdam, 1983), pp. 55-64.
602 !
603 ! *Long Description:
604 ! The following complete description of the user interface to
605 ! DLSODE consists of four parts:
606 !
607 ! 1. The call sequence to subroutine DLSODE, which is a driver
608 ! routine for the solver. This includes descriptions of both
609 ! the call sequence arguments and user-supplied routines.
610 ! Following these descriptions is a description of optional
611 ! inputs available through the call sequence, and then a
612 ! description of optional outputs in the work arrays.
613 !
614 ! 2. Descriptions of other routines in the DLSODE package that may
615 ! be (optionally) called by the user. These provide the ability
616 ! to alter error message handling, save and restore the internal
617 ! COMMON, and obtain specified derivatives of the solution y(t).
618 !
619 ! 3. Descriptions of COMMON block to be declared in overlay or
620 ! similar environments, or to be saved when doing an interrupt
621 ! of the problem and continued solution later.
622 !
623 ! 4. Description of two routines in the DLSODE package, either of
624 ! which the user may replace with his own version, if desired.
625 ! These relate to the measurement of errors.
626 !
627 !
628 ! Part 1. Call Sequence
629 ! ----------------------
630 !
631 ! Arguments
632 ! ---------
633 ! The call sequence parameters used for input only are
634 !
635 ! F, NEQ, TOUT, ITOL, RelTol, AbsTol, ITASK, IOPT, LRW, LIW, JAC, MF,
636 !
637 ! and those used for both input and output are
638 !
639 ! Y, T, ISTATE.
640 !
641 ! The work arrays RWORK and IWORK are also used for conditional and
642 ! optional inputs and optional outputs. (The term output here
643 ! refers to the return from subroutine DLSODE to the user's calling
644 ! program.)
645 !
646 ! The legality of input parameters will be thoroughly checked on the
647 ! initial call for the problem, but not checked thereafter unless a
648 ! change in input parameters is flagged by ISTATE = 3 on input.
649 !
650 ! The descriptions of the call arguments are as follows.
651 !
652 ! F The name of the user-supplied subroutine defining the ODE
653 ! system. The system must be put in the first-order form
654 ! dy/dt = f(t,y), where f is a vector-valued function of
655 ! the scalar t and the vector y. Subroutine F is to compute
656 ! the function f. It is to have the form
657 !
658 ! SUBROUTINE F (NEQ, T, Y, YDOT)
659 ! KPP_REAL T, Y(*), YDOT(*)
660 !
661 ! where NEQ, T, and Y are input, and the array YDOT =
662 ! f(T,Y) is output. Y and YDOT are arrays of length NEQ.
663 ! Subroutine F should not alter Y(1),...,Y(NEQ). F must be
664 ! declared EXTERNAL in the calling program.
665 !
666 ! Subroutine F may access user-defined quantities in
667 ! NEQ(2),... and/or in Y(NEQ+1),..., if NEQ is an array
668 ! (dimensioned in F) and/or Y has length exceeding NEQ.
669 ! See the descriptions of NEQ and Y below.
670 !
671 ! If quantities computed in the F routine are needed
672 ! externally to DLSODE, an extra call to F should be made
673 ! for this purpose, for consistent and accurate results.
674 ! If only the derivative dy/dt is needed, use DINTDY
675 ! instead.
676 !
677 ! NEQ The size of the ODE system (number of first-order
678 ! ordinary differential equations). Used only for input.
679 ! NEQ may be decreased, but not increased, during the
680 ! problem. If NEQ is decreased (with ISTATE = 3 on input),
681 ! the remaining components of Y should be left undisturbed,
682 ! if these are to be accessed in F and/or JAC.
683 !
684 ! Normally, NEQ is a scalar, and it is generally referred
685 ! to as a scalar in this user interface description.
686 ! However, NEQ may be an array, with NEQ set to the
687 ! system size. (The DLSODE package accesses only NEQ.)
688 ! In either case, this parameter is passed as the NEQ
689 ! argument in all calls to F and JAC. Hence, if it is an
690 ! array, locations NEQ(2),... may be used to store other
691 ! integer data and pass it to F and/or JAC. Subroutines
692 ! F and/or JAC must include NEQ in a DIMENSION statement
693 ! in that case.
694 !
695 ! Y A real array for the vector of dependent variables, of
696 ! length NEQ or more. Used for both input and output on
697 ! the first call (ISTATE = 1), and only for output on
698 ! other calls. On the first call, Y must contain the
699 ! vector of initial values. On output, Y contains the
700 ! computed solution vector, evaluated at T. If desired,
701 ! the Y array may be used for other purposes between
702 ! calls to the solver.
703 !
704 ! This array is passed as the Y argument in all calls to F
705 ! and JAC. Hence its length may exceed NEQ, and locations
706 ! Y(NEQ+1),... may be used to store other real data and
707 ! pass it to F and/or JAC. (The DLSODE package accesses
708 ! only Y(1),...,Y(NEQ).)
709 !
710 ! T The independent variable. On input, T is used only on
711 ! the first call, as the initial point of the integration.
712 ! On output, after each call, T is the value at which a
713 ! computed solution Y is evaluated (usually the same as
714 ! TOUT). On an error return, T is the farthest point
715 ! reached.
716 !
717 ! TOUT The next value of T at which a computed solution is
718 ! desired. Used only for input.
719 !
720 ! When starting the problem (ISTATE = 1), TOUT may be equal
721 ! to T for one call, then should not equal T for the next
722 ! call. For the initial T, an input value of TOUT .NE. T
723 ! is used in order to determine the direction of the
724 ! integration (i.e., the algebraic sign of the step sizes)
725 ! and the rough scale of the problem. Integration in
726 ! either direction (forward or backward in T) is permitted.
727 !
728 ! If ITASK = 2 or 5 (one-step modes), TOUT is ignored
729 ! after the first call (i.e., the first call with
730 ! TOUT .NE. T). Otherwise, TOUT is required on every call.
731 !
732 ! If ITASK = 1, 3, or 4, the values of TOUT need not be
733 ! monotone, but a value of TOUT which backs up is limited
734 ! to the current internal T interval, whose endpoints are
735 ! TCUR - HU and TCUR. (See "Optional Outputs" below for
736 ! TCUR and HU.)
737 !
738 !
739 ! ITOL An indicator for the type of error control. See
740 ! description below under AbsTol. Used only for input.
741 !
742 ! RelTol A relative error tolerance parameter, either a scalar or
743 ! an array of length NEQ. See description below under
744 ! AbsTol. Input only.
745 !
746 ! AbsTol An absolute error tolerance parameter, either a scalar or
747 ! an array of length NEQ. Input only.
748 !
749 ! The input parameters ITOL, RelTol, and AbsTol determine the
750 ! error control performed by the solver. The solver will
751 ! control the vector e = (e(i)) of estimated local errors
752 ! in Y, according to an inequality of the form
753 !
754 ! rms-norm of ( e(i)/EWT(i) ) <= 1,
755 !
756 ! where
757 !
758 ! EWT(i) = RelTol(i)*ABS(Y(i)) + AbsTol(i),
759 !
760 ! and the rms-norm (root-mean-square norm) here is
761 !
762 ! rms-norm(v) = SQRT(sum v(i)**2 / NEQ).
763 !
764 ! Here EWT = (EWT(i)) is a vector of weights which must
765 ! always be positive, and the values of RelTol and AbsTol
766 ! should all be nonnegative. The following table gives the
767 ! types (scalar/array) of RelTol and AbsTol, and the
768 ! corresponding form of EWT(i).
769 !
770 ! ITOL RelTol AbsTol EWT(i)
771 ! ---- ------ ------ -----------------------------
772 ! 1 scalar scalar RelTol*ABS(Y(i)) + AbsTol
773 ! 2 scalar array RelTol*ABS(Y(i)) + AbsTol(i)
774 ! 3 array scalar RelTol(i)*ABS(Y(i)) + AbsTol
775 ! 4 array array RelTol(i)*ABS(Y(i)) + AbsTol(i)
776 !
777 ! When either of these parameters is a scalar, it need not
778 ! be dimensioned in the user's calling program.
779 !
780 ! If none of the above choices (with ITOL, RelTol, and AbsTol
781 ! fixed throughout the problem) is suitable, more general
782 ! error controls can be obtained by substituting
783 ! user-supplied routines for the setting of EWT and/or for
784 ! the norm calculation. See Part 4 below.
785 !
786 ! If global errors are to be estimated by making a repeated
787 ! run on the same problem with smaller tolerances, then all
788 ! components of RelTol and AbsTol (i.e., of EWT) should be
789 ! scaled down uniformly.
790 !
791 ! ITASK An index specifying the task to be performed. Input
792 ! only. ITASK has the following values and meanings:
793 ! 1 Normal computation of output values of y(t) at
794 ! t = TOUT (by overshooting and interpolating).
795 ! 2 Take one step only and return.
796 ! 3 Stop at the first internal mesh point at or beyond
797 ! t = TOUT and return.
798 ! 4 Normal computation of output values of y(t) at
799 ! t = TOUT but without overshooting t = TCRIT. TCRIT
800 ! must be input as RWORK(1). TCRIT may be equal to or
801 ! beyond TOUT, but not behind it in the direction of
802 ! integration. This option is useful if the problem
803 ! has a singularity at or beyond t = TCRIT.
804 ! 5 Take one step, without passing TCRIT, and return.
805 ! TCRIT must be input as RWORK(1).
806 !
807 ! Note: If ITASK = 4 or 5 and the solver reaches TCRIT
808 ! (within roundoff), it will return T = TCRIT (exactly) to
809 ! indicate this (unless ITASK = 4 and TOUT comes before
810 ! TCRIT, in which case answers at T = TOUT are returned
811 ! first).
812 !
813 ! ISTATE An index used for input and output to specify the state
814 ! of the calculation.
815 !
816 ! On input, the values of ISTATE are as follows:
817 ! 1 This is the first call for the problem
818 ! (initializations will be done). See "Note" below.
819 ! 2 This is not the first call, and the calculation is to
820 ! continue normally, with no change in any input
821 ! parameters except possibly TOUT and ITASK. (If ITOL,
822 ! RelTol, and/or AbsTol are changed between calls with
823 ! ISTATE = 2, the new values will be used but not
824 ! tested for legality.)
825 ! 3 This is not the first call, and the calculation is to
826 ! continue normally, but with a change in input
827 ! parameters other than TOUT and ITASK. Changes are
828 ! allowed in NEQ, ITOL, RelTol, AbsTol, IOPT, LRW, LIW, MF,
829 ! ML, MU, and any of the optional inputs except H0.
830 ! (See IWORK description for ML and MU.)
831 !
832 ! Note: A preliminary call with TOUT = T is not counted as
833 ! a first call here, as no initialization or checking of
834 ! input is done. (Such a call is sometimes useful for the
835 ! purpose of outputting the initial conditions.) Thus the
836 ! first call for which TOUT .NE. T requires ISTATE = 1 on
837 ! input.
838 !
839 ! On output, ISTATE has the following values and meanings:
840 ! 1 Nothing was done, as TOUT was equal to T with
841 ! ISTATE = 1 on input.
842 ! 2 The integration was performed successfully.
843 ! -1 An excessive amount of work (more than MXSTEP steps)
844 ! was done on this call, before completing the
845 ! requested task, but the integration was otherwise
846 ! successful as far as T. (MXSTEP is an optional input
847 ! and is normally 500.) To continue, the user may
848 ! simply reset ISTATE to a value >1 and call again (the
849 ! excess work step counter will be reset to 0). In
850 ! addition, the user may increase MXSTEP to avoid this
851 ! error return; see "Optional Inputs" below.
852 ! -2 Too much accuracy was requested for the precision of
853 ! the machine being used. This was detected before
854 ! completing the requested task, but the integration
855 ! was successful as far as T. To continue, the
856 ! tolerance parameters must be reset, and ISTATE must
857 ! be set to 3. The optional output TOLSF may be used
858 ! for this purpose. (Note: If this condition is
859 ! detected before taking any steps, then an illegal
860 ! input return (ISTATE = -3) occurs instead.)
861 ! -3 Illegal input was detected, before taking any
862 ! integration steps. See written message for details.
863 ! (Note: If the solver detects an infinite loop of
864 ! calls to the solver with illegal input, it will cause
865 ! the run to stop.)
866 ! -4 There were repeated error-test failures on one
867 ! attempted step, before completing the requested task,
868 ! but the integration was successful as far as T. The
869 ! problem may have a singularity, or the input may be
870 ! inappropriate.
871 ! -5 There were repeated convergence-test failures on one
872 ! attempted step, before completing the requested task,
873 ! but the integration was successful as far as T. This
874 ! may be caused by an inaccurate Jacobian matrix, if
875 ! one is being used.
876 ! -6 EWT(i) became zero for some i during the integration.
877 ! Pure relative error control (AbsTol(i)=0.0) was
878 ! requested on a variable which has now vanished. The
879 ! integration was successful as far as T.
880 !
881 ! Note: Since the normal output value of ISTATE is 2, it
882 ! does not need to be reset for normal continuation. Also,
883 ! since a negative input value of ISTATE will be regarded
884 ! as illegal, a negative output value requires the user to
885 ! change it, and possibly other inputs, before calling the
886 ! solver again.
887 !
888 ! IOPT An integer flag to specify whether any optional inputs
889 ! are being used on this call. Input only. The optional
890 ! inputs are listed under a separate heading below.
891 ! 0 No optional inputs are being used. Default values
892 ! will be used in all cases.
893 ! 1 One or more optional inputs are being used.
894 !
895 ! RWORK A real working array (double precision). The length of
896 ! RWORK must be at least
897 !
898 ! 20 + NYH*(MAXORD + 1) + 3*NEQ + LWM
899 !
900 ! where
901 ! NYH = the initial value of NEQ,
902 ! MAXORD = 12 (if METH = 1) or 5 (if METH = 2) (unless a
903 ! smaller value is given as an optional input),
904 ! LWM = 0 if MITER = 0,
905 ! LWM = NEQ**2 + 2 if MITER = 1 or 2,
906 ! LWM = NEQ + 2 if MITER = 3, and
907 ! LWM = (2*ML + MU + 1)*NEQ + 2
908 ! if MITER = 4 or 5.
909 ! (See the MF description below for METH and MITER.)
910 !
911 ! Thus if MAXORD has its default value and NEQ is constant,
912 ! this length is:
913 ! 20 + 16*NEQ for MF = 10,
914 ! 22 + 16*NEQ + NEQ**2 for MF = 11 or 12,
915 ! 22 + 17*NEQ for MF = 13,
916 ! 22 + 17*NEQ + (2*ML + MU)*NEQ for MF = 14 or 15,
917 ! 20 + 9*NEQ for MF = 20,
918 ! 22 + 9*NEQ + NEQ**2 for MF = 21 or 22,
919 ! 22 + 10*NEQ for MF = 23,
920 ! 22 + 10*NEQ + (2*ML + MU)*NEQ for MF = 24 or 25.
921 !
922 ! The first 20 words of RWORK are reserved for conditional
923 ! and optional inputs and optional outputs.
924 !
925 ! The following word in RWORK is a conditional input:
926 ! RWORK(1) = TCRIT, the critical value of t which the
927 ! solver is not to overshoot. Required if ITASK
928 ! is 4 or 5, and ignored otherwise. See ITASK.
929 !
930 ! LRW The length of the array RWORK, as declared by the user.
931 ! (This will be checked by the solver.)
932 !
933 ! IWORK An integer work array. Its length must be at least
934 ! 20 if MITER = 0 or 3 (MF = 10, 13, 20, 23), or
935 ! 20 + NEQ otherwise (MF = 11, 12, 14, 15, 21, 22, 24, 25).
936 ! (See the MF description below for MITER.) The first few
937 ! words of IWORK are used for conditional and optional
938 ! inputs and optional outputs.
939 !
940 ! The following two words in IWORK are conditional inputs:
941 ! IWORK(1) = ML These are the lower and upper half-
942 ! IWORK(2) = MU bandwidths, respectively, of the banded
943 ! Jacobian, excluding the main diagonal.
944 ! The band is defined by the matrix locations
945 ! (i,j) with i - ML <= j <= i + MU. ML and MU
946 ! must satisfy 0 <= ML,MU <= NEQ - 1. These are
947 ! required if MITER is 4 or 5, and ignored
948 ! otherwise. ML and MU may in fact be the band
949 ! parameters for a matrix to which df/dy is only
950 ! approximately equal.
951 !
952 ! LIW The length of the array IWORK, as declared by the user.
953 ! (This will be checked by the solver.)
954 !
955 ! Note: The work arrays must not be altered between calls to DLSODE
956 ! for the same problem, except possibly for the conditional and
957 ! optional inputs, and except for the last 3*NEQ words of RWORK.
958 ! The latter space is used for internal scratch space, and so is
959 ! available for use by the user outside DLSODE between calls, if
960 ! desired (but not for use by F or JAC).
961 !
962 ! JAC The name of the user-supplied routine (MITER = 1 or 4) to
963 ! compute the Jacobian matrix, df/dy, as a function of the
964 ! scalar t and the vector y. (See the MF description below
965 ! for MITER.) It is to have the form
966 !
967 ! SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD)
968 ! KPP_REAL T, Y(*), PD(NROWPD,*)
969 !
970 ! where NEQ, T, Y, ML, MU, and NROWPD are input and the
971 ! array PD is to be loaded with partial derivatives
972 ! (elements of the Jacobian matrix) on output. PD must be
973 ! given a first dimension of NROWPD. T and Y have the same
974 ! meaning as in subroutine F.
975 !
976 ! In the full matrix case (MITER = 1), ML and MU are
977 ! ignored, and the Jacobian is to be loaded into PD in
978 ! columnwise manner, with df(i)/dy(j) loaded into PD(i,j).
979 !
980 ! In the band matrix case (MITER = 4), the elements within
981 ! the band are to be loaded into PD in columnwise manner,
982 ! with diagonal lines of df/dy loaded into the rows of PD.
983 ! Thus df(i)/dy(j) is to be loaded into PD(i-j+MU+1,j). ML
984 ! and MU are the half-bandwidth parameters (see IWORK).
985 ! The locations in PD in the two triangular areas which
986 ! correspond to nonexistent matrix elements can be ignored
987 ! or loaded arbitrarily, as they are overwritten by DLSODE.
988 !
989 ! JAC need not provide df/dy exactly. A crude approximation
990 ! (possibly with a smaller bandwidth) will do.
991 !
992 ! In either case, PD is preset to zero by the solver, so
993 ! that only the nonzero elements need be loaded by JAC.
994 ! Each call to JAC is preceded by a call to F with the same
995 ! arguments NEQ, T, and Y. Thus to gain some efficiency,
996 ! intermediate quantities shared by both calculations may
997 ! be saved in a user COMMON block by F and not recomputed
998 ! by JAC, if desired. Also, JAC may alter the Y array, if
999 ! desired. JAC must be declared EXTERNAL in the calling
1000 ! program.
1001 !
1002 ! Subroutine JAC may access user-defined quantities in
1003 ! NEQ(2),... and/or in Y(NEQ+1),... if NEQ is an array
1004 ! (dimensioned in JAC) and/or Y has length exceeding
1005 ! NEQ. See the descriptions of NEQ and Y above.
1006 !
1007 ! MF The method flag. Used only for input. The legal values
1008 ! of MF are 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24,
1009 ! and 25. MF has decimal digits METH and MITER:
1010 ! MF = 10*METH + MITER .
1011 !
1012 ! METH indicates the basic linear multistep method:
1013 ! 1 Implicit Adams method.
1014 ! 2 Method based on backward differentiation formulas
1015 ! (BDF's).
1016 !
1017 ! MITER indicates the corrector iteration method:
1018 ! 0 Functional iteration (no Jacobian matrix is
1019 ! involved).
1020 ! 1 Chord iteration with a user-supplied full (NEQ by
1021 ! NEQ) Jacobian.
1022 ! 2 Chord iteration with an internally generated
1023 ! (difference quotient) full Jacobian (using NEQ
1024 ! extra calls to F per df/dy value).
1025 ! 3 Chord iteration with an internally generated
1026 ! diagonal Jacobian approximation (using one extra call
1027 ! to F per df/dy evaluation).
1028 ! 4 Chord iteration with a user-supplied banded Jacobian.
1029 ! 5 Chord iteration with an internally generated banded
1030 ! Jacobian (using ML + MU + 1 extra calls to F per
1031 ! df/dy evaluation).
1032 !
1033 ! If MITER = 1 or 4, the user must supply a subroutine JAC
1034 ! (the name is arbitrary) as described above under JAC.
1035 ! For other values of MITER, a dummy argument can be used.
1036 !
1037 ! Optional Inputs
1038 ! ---------------
1039 ! The following is a list of the optional inputs provided for in the
1040 ! call sequence. (See also Part 2.) For each such input variable,
1041 ! this table lists its name as used in this documentation, its
1042 ! location in the call sequence, its meaning, and the default value.
1043 ! The use of any of these inputs requires IOPT = 1, and in that case
1044 ! all of these inputs are examined. A value of zero for any of
1045 ! these optional inputs will cause the default value to be used.
1046 ! Thus to use a subset of the optional inputs, simply preload
1047 ! locations 5 to 10 in RWORK and IWORK to 0.0 and 0 respectively,
1048 ! and then set those of interest to nonzero values.
1049 !
1050 ! Name Location Meaning and default value
1051 ! ------ --------- -----------------------------------------------
1052 ! H0 RWORK(5) Step size to be attempted on the first step.
1053 ! The default value is determined by the solver.
1054 ! HMAX RWORK(6) Maximum absolute step size allowed. The
1055 ! default value is infinite.
1056 ! HMIN RWORK(7) Minimum absolute step size allowed. The
1057 ! default value is 0. (This lower bound is not
1058 ! enforced on the final step before reaching
1059 ! TCRIT when ITASK = 4 or 5.)
1060 ! MAXORD IWORK(5) Maximum order to be allowed. The default value
1061 ! is 12 if METH = 1, and 5 if METH = 2. (See the
1062 ! MF description above for METH.) If MAXORD
1063 ! exceeds the default value, it will be reduced
1064 ! to the default value. If MAXORD is changed
1065 ! during the problem, it may cause the current
1066 ! order to be reduced.
1067 ! MXSTEP IWORK(6) Maximum number of (internally defined) steps
1068 ! allowed during one call to the solver. The
1069 ! default value is 500.
1070 ! MXHNIL IWORK(7) Maximum number of messages printed (per
1071 ! problem) warning that T + H = T on a step
1072 ! (H = step size). This must be positive to
1073 ! result in a nondefault value. The default
1074 ! value is 10.
1075 !
1076 ! Optional Outputs
1077 ! ----------------
1078 ! As optional additional output from DLSODE, the variables listed
1079 ! below are quantities related to the performance of DLSODE which
1080 ! are available to the user. These are communicated by way of the
1081 ! work arrays, but also have internal mnemonic names as shown.
1082 ! Except where stated otherwise, all of these outputs are defined on
1083 ! any successful return from DLSODE, and on any return with ISTATE =
1084 ! -1, -2, -4, -5, or -6. On an illegal input return (ISTATE = -3),
1085 ! they will be unchanged from their existing values (if any), except
1086 ! possibly for TOLSF, LENRW, and LENIW. On any error return,
1087 ! outputs relevant to the error will be defined, as noted below.
1088 !
1089 ! Name Location Meaning
1090 ! ----- --------- ------------------------------------------------
1091 ! HU RWORK(11) Step size in t last used (successfully).
1092 ! HCUR RWORK(12) Step size to be attempted on the next step.
1093 ! TCUR RWORK(13) Current value of the independent variable which
1094 ! the solver has actually reached, i.e., the
1095 ! current internal mesh point in t. On output,
1096 ! TCUR will always be at least as far as the
1097 ! argument T, but may be farther (if interpolation
1098 ! was done).
1099 ! TOLSF RWORK(14) Tolerance scale factor, greater than 1.0,
1100 ! computed when a request for too much accuracy
1101 ! was detected (ISTATE = -3 if detected at the
1102 ! start of the problem, ISTATE = -2 otherwise).
1103 ! If ITOL is left unaltered but RelTol and AbsTol are
1104 ! uniformly scaled up by a factor of TOLSF for the
1105 ! next call, then the solver is deemed likely to
1106 ! succeed. (The user may also ignore TOLSF and
1107 ! alter the tolerance parameters in any other way
1108 ! appropriate.)
1109 ! NST IWORK(11) Number of steps taken for the problem so far.
1110 ! NFE IWORK(12) Number of F evaluations for the problem so far.
1111 ! NJE IWORK(13) Number of Jacobian evaluations (and of matrix LU
1112 ! decompositions) for the problem so far.
1113 ! NQU IWORK(14) Method order last used (successfully).
1114 ! NQCUR IWORK(15) Order to be attempted on the next step.
1115 ! IMXER IWORK(16) Index of the component of largest magnitude in
1116 ! the weighted local error vector ( e(i)/EWT(i) ),
1117 ! on an error return with ISTATE = -4 or -5.
1118 ! LENRW IWORK(17) Length of RWORK actually required. This is
1119 ! defined on normal returns and on an illegal
1120 ! input return for insufficient storage.
1121 ! LENIW IWORK(18) Length of IWORK actually required. This is
1122 ! defined on normal returns and on an illegal
1123 ! input return for insufficient storage.
1124 !
1125 ! The following two arrays are segments of the RWORK array which may
1126 ! also be of interest to the user as optional outputs. For each
1127 ! array, the table below gives its internal name, its base address
1128 ! in RWORK, and its description.
1129 !
1130 ! Name Base address Description
1131 ! ---- ------------ ----------------------------------------------
1132 ! YH 21 The Nordsieck history array, of size NYH by
1133 ! (NQCUR + 1), where NYH is the initial value of
1134 ! NEQ. For j = 0,1,...,NQCUR, column j + 1 of
1135 ! YH contains HCUR**j/factorial(j) times the jth
1136 ! derivative of the interpolating polynomial
1137 ! currently representing the solution, evaluated
1138 ! at t = TCUR.
1139 ! ACOR LENRW-NEQ+1 Array of size NEQ used for the accumulated
1140 ! corrections on each step, scaled on output to
1141 ! represent the estimated local error in Y on
1142 ! the last step. This is the vector e in the
1143 ! description of the error control. It is
1144 ! defined only on successful return from DLSODE.
1145 !
1146 !
1147 ! Part 2. Other Callable Routines
1148 ! --------------------------------
1149 !
1150 ! The following are optional calls which the user may make to gain
1151 ! additional capabilities in conjunction with DLSODE.
1152 !
1153 ! Form of call Function
1154 ! ------------------------ ----------------------------------------
1155 ! CALL XSETUN(LUN) Set the logical unit number, LUN, for
1156 ! output of messages from DLSODE, if the
1157 ! default is not desired. The default
1158 ! value of LUN is 6. This call may be made
1159 ! at any time and will take effect
1160 ! immediately.
1161 ! CALL XSETF(MFLAG) Set a flag to control the printing of
1162 ! messages by DLSODE. MFLAG = 0 means do
1163 ! not print. (Danger: this risks losing
1164 ! valuable information.) MFLAG = 1 means
1165 ! print (the default). This call may be
1166 ! made at any time and will take effect
1167 ! immediately.
1168 ! CALL DSRCOM(RSAV,ISAV,JOB) Saves and restores the contents of the
1169 ! internal COMMON blocks used by DLSODE
1170 ! (see Part 3 below). RSAV must be a
1171 ! real array of length 218 or more, and
1172 ! ISAV must be an integer array of length
1173 ! 37 or more. JOB = 1 means save COMMON
1174 ! into RSAV/ISAV. JOB = 2 means restore
1175 ! COMMON from same. DSRCOM is useful if
1176 ! one is interrupting a run and restarting
1177 ! later, or alternating between two or
1178 ! more problems solved with DLSODE.
1179 ! CALL DINTDY(,,,,,) Provide derivatives of y, of various
1180 ! (see below) orders, at a specified point t, if
1181 ! desired. It may be called only after a
1182 ! successful return from DLSODE. Detailed
1183 ! instructions follow.
1184 !
1185 ! Detailed instructions for using DINTDY
1186 ! --------------------------------------
1187 ! The form of the CALL is:
1188 !
1189 ! CALL DINTDY (T, K, RWORK(21), NYH, DKY, IFLAG)
1190 !
1191 ! The input parameters are:
1192 !
1193 ! T Value of independent variable where answers are
1194 ! desired (normally the same as the T last returned by
1195 ! DLSODE). For valid results, T must lie between
1196 ! TCUR - HU and TCUR. (See "Optional Outputs" above
1197 ! for TCUR and HU.)
1198 ! K Integer order of the derivative desired. K must
1199 ! satisfy 0 <= K <= NQCUR, where NQCUR is the current
1200 ! order (see "Optional Outputs"). The capability
1201 ! corresponding to K = 0, i.e., computing y(t), is
1202 ! already provided by DLSODE directly. Since
1203 ! NQCUR >= 1, the first derivative dy/dt is always
1204 ! available with DINTDY.
1205 ! RWORK(21) The base address of the history array YH.
1206 ! NYH Column length of YH, equal to the initial value of NEQ.
1207 !
1208 ! The output parameters are:
1209 !
1210 ! DKY Real array of length NEQ containing the computed value
1211 ! of the Kth derivative of y(t).
1212 ! IFLAG Integer flag, returned as 0 if K and T were legal,
1213 ! -1 if K was illegal, and -2 if T was illegal.
1214 ! On an error return, a message is also written.
1215 !
1216 !
1217 ! Part 3. Common Blocks
1218 ! ----------------------
1219 !
1220 ! If DLSODE is to be used in an overlay situation, the user must
1221 ! declare, in the primary overlay, the variables in:
1222 ! (1) the call sequence to DLSODE,
1223 ! (2) the internal COMMON block /DLS001/, of length 255
1224 ! (218 double precision words followed by 37 integer words).
1225 !
1226 ! If DLSODE is used on a system in which the contents of internal
1227 ! COMMON blocks are not preserved between calls, the user should
1228 ! declare the above COMMON block in his main program to insure that
1229 ! its contents are preserved.
1230 !
1231 ! If the solution of a given problem by DLSODE is to be interrupted
1232 ! and then later continued, as when restarting an interrupted run or
1233 ! alternating between two or more problems, the user should save,
1234 ! following the return from the last DLSODE call prior to the
1235 ! interruption, the contents of the call sequence variables and the
1236 ! internal COMMON block, and later restore these values before the
1237 ! next DLSODE call for that problem. In addition, if XSETUN and/or
1238 ! XSETF was called for non-default handling of error messages, then
1239 ! these calls must be repeated. To save and restore the COMMON
1240 ! block, use subroutine DSRCOM (see Part 2 above).
1241 !
1242 !
1243 ! Part 4. Optionally Replaceable Solver Routines
1244 ! -----------------------------------------------
1245 !
1246 ! Below are descriptions of two routines in the DLSODE package which
1247 ! relate to the measurement of errors. Either routine can be
1248 ! replaced by a user-supplied version, if desired. However, since
1249 ! such a replacement may have a major impact on performance, it
1250 ! should be done only when absolutely necessary, and only with great
1251 ! caution. (Note: The means by which the package version of a
1252 ! routine is superseded by the user's version may be system-
1253 ! dependent.)
1254 !
1255 ! DEWSET
1256 ! ------
1257 ! The following subroutine is called just before each internal
1258 ! integration step, and sets the array of error weights, EWT, as
1259 ! described under ITOL/RelTol/AbsTol above:
1260 !
1261 ! SUBROUTINE DEWSET (NEQ, ITOL, RelTol, AbsTol, YCUR, EWT)
1262 !
1263 ! where NEQ, ITOL, RelTol, and AbsTol are as in the DLSODE call
1264 ! sequence, YCUR contains the current dependent variable vector,
1265 ! and EWT is the array of weights set by DEWSET.
1266 !
1267 ! If the user supplies this subroutine, it must return in EWT(i)
1268 ! (i = 1,...,NEQ) a positive quantity suitable for comparing errors
1269 ! in Y(i) to. The EWT array returned by DEWSET is passed to the
1270 ! DVNORM routine (see below), and also used by DLSODE in the
1271 ! computation of the optional output IMXER, the diagonal Jacobian
1272 ! approximation, and the increments for difference quotient
1273 ! Jacobians.
1274 !
1275 ! In the user-supplied version of DEWSET, it may be desirable to use
1276 ! the current values of derivatives of y. Derivatives up to order NQ
1277 ! are available from the history array YH, described above under
1278 ! optional outputs. In DEWSET, YH is identical to the YCUR array,
1279 ! extended to NQ + 1 columns with a column length of NYH and scale
1280 ! factors of H**j/factorial(j). On the first call for the problem,
1281 ! given by NST = 0, NQ is 1 and H is temporarily set to 1.0.
1282 ! NYH is the initial value of NEQ. The quantities NQ, H, and NST
1283 ! can be obtained by including in SEWSET the statements:
1284 ! KPP_REAL RLS
1285 ! COMMON /DLS001/ RLS(218),ILS(37)
1286 ! NQ = ILS(33)
1287 ! NST = ILS(34)
1288 ! H = RLS(212)
1289 ! Thus, for example, the current value of dy/dt can be obtained as
1290 ! YCUR(NYH+i)/H (i=1,...,NEQ) (and the division by H is unnecessary
1291 ! when NST = 0).
1292 !
1293 ! DVNORM
1294 ! ------
1295 ! DVNORM is a real function routine which computes the weighted
1296 ! root-mean-square norm of a vector v:
1297 !
1298 ! d = DVNORM (n, v, w)
1299 !
1300 ! where:
1301 ! n = the length of the vector,
1302 ! v = real array of length n containing the vector,
1303 ! w = real array of length n containing weights,
1304 ! d = SQRT( (1/n) * sum(v(i)*w(i))**2 ).
1305 !
1306 ! DVNORM is called with n = NEQ and with w(i) = 1.0/EWT(i), where
1307 ! EWT is as set by subroutine DEWSET.
1308 !
1309 ! If the user supplies this function, it should return a nonnegative
1310 ! value of DVNORM suitable for use in the error control in DLSODE.
1311 ! None of the arguments should be altered by DVNORM. For example, a
1312 ! user-supplied DVNORM routine might:
1313 ! - Substitute a max-norm of (v(i)*w(i)) for the rms-norm, or
1314 ! - Ignore some components of v in the norm, with the effect of
1315 ! suppressing the error control on those components of Y.
1316 ! ---------------------------------------------------------------------
1317 !***ROUTINES CALLED DEWSET, DINTDY, DUMACH, DSTODE, DVNORM, XERRWD
1318 !***COMMON BLOCKS DLS001
1319 !***REVISION HISTORY (YYYYMMDD)
1320 ! 19791129 DATE WRITTEN
1321 ! 19791213 Minor changes to declarations; DELP init. in STODE.
1322 ! 19800118 Treat NEQ as array; integer declarations added throughout;
1323 ! minor changes to prologue.
1324 ! 19800306 Corrected TESCO(1,NQP1) setting in CFODE.
1325 ! 19800519 Corrected access of YH on forced order reduction;
1326 ! numerous corrections to prologues and other comments.
1327 ! 19800617 In main driver, added loading of SQRT(UROUND) in RWORK;
1328 ! minor corrections to main prologue.
1329 ! 19800923 Added zero initialization of HU and NQU.
1330 ! 19801218 Revised XERRWD routine; minor corrections to main prologue.
1331 ! 19810401 Minor changes to comments and an error message.
1332 ! 19810814 Numerous revisions: replaced EWT by 1/EWT; used flags
1333 ! JCUR, ICF, IERPJ, IERSL between STODE and subordinates;
1334 ! added tuning parameters CCMAX, MAXCOR, MSBP, MXNCF;
1335 ! reorganized returns from STODE; reorganized type decls.;
1336 ! fixed message length in XERRWD; changed default LUNIT to 6;
1337 ! changed Common lengths; changed comments throughout.
1338 ! 19870330 Major update by ACH: corrected comments throughout;
1339 ! removed TRET from Common; rewrote EWSET with 4 loops;
1340 ! fixed t test in INTDY; added Cray directives in STODE;
1341 ! in STODE, fixed DELP init. and logic around PJAC call;
1342 ! combined routines to save/restore Common;
1343 ! passed LEVEL = 0 in error message calls (except run abort).
1344 ! 19890426 Modified prologue to SLATEC/LDOC format. (FNF)
1345 ! 19890501 Many improvements to prologue. (FNF)
1346 ! 19890503 A few final corrections to prologue. (FNF)
1347 ! 19890504 Minor cosmetic changes. (FNF)
1348 ! 19890510 Corrected description of Y in Arguments section. (FNF)
1349 ! 19890517 Minor corrections to prologue. (FNF)
1350 ! 19920514 Updated with prologue edited 891025 by G. Shaw for manual.
1351 ! 19920515 Converted source lines to upper case. (FNF)
1352 ! 19920603 Revised XERRWD calls using mixed upper-lower case. (ACH)
1353 ! 19920616 Revised prologue comment regarding CFT. (ACH)
1354 ! 19921116 Revised prologue comments regarding Common. (ACH).
1355 ! 19930326 Added comment about non-reentrancy. (FNF)
1356 ! 19930723 Changed D1MACH to DUMACH. (FNF)
1357 ! 19930801 Removed ILLIN and NTREP from Common (affects driver logic);
1358 ! minor changes to prologue and internal comments;
1359 ! changed Hollerith strings to quoted strings;
1360 ! changed internal comments to mixed case;
1361 ! replaced XERRWD with new version using character type;
1362 ! changed dummy dimensions from 1 to *. (ACH)
1363 ! 19930809 Changed to generic intrinsic names; changed names of
1364 ! subprograms and Common blocks to DLSODE etc. (ACH)
1365 ! 19930929 Eliminated use of REAL intrinsic; other minor changes. (ACH)
1366 ! 20010412 Removed all 'own' variables from Common block /DLS001/
1367 ! (affects declarations in 6 routines). (ACH)
1368 ! 20010509 Minor corrections to prologue. (ACH)
1369 ! 20031105 Restored 'own' variables to Common block /DLS001/, to
1370 ! enable interrupt/restart feature. (ACH)
1371 ! 20031112 Added SAVE statements for data-loaded constants.
1372 !
1373 !***END PROLOGUE DLSODE
1374 !
1375 !*Internal Notes:
1376 !
1377 ! Other Routines in the DLSODE Package.
1378 !
1379 ! In addition to Subroutine DLSODE, the DLSODE package includes the
1380 ! following subroutines and function routines:
1381 ! DINTDY computes an interpolated value of the y vector at t = TOUT.
1382 ! DSTODE is the core integrator, which does one step of the
1383 ! integration and the associated error control.
1384 ! DCFODE sets all method coefficients and test constants.
1385 ! DPREPJ computes and preprocesses the Jacobian matrix J = df/dy
1386 ! and the Newton iteration matrix P = I - h*l0*J.
1387 ! DSOLSY manages solution of linear system in chord iteration.
1388 ! DEWSET sets the error weight vector EWT before each step.
1389 ! DVNORM computes the weighted R.M.S. norm of a vector.
1390 ! DSRCOM is a user-callable routine to save and restore
1391 ! the contents of the internal Common block.
1392 ! DGEFA and DGESL are routines from LINPACK for solving full
1393 ! systems of linear algebraic equations.
1394 ! DGBFA and DGBSL are routines from LINPACK for solving banded
1395 ! linear systems.
1396 ! DUMACH computes the unit roundoff in a machine-independent manner.
1397 ! XERRWD, XSETUN, XSETF, IXSAV, IUMACH handle the printing of all
1398 ! error messages and warnings. XERRWD is machine-dependent.
1399 ! Note: DVNORM, DUMACH, IXSAV, and IUMACH are function routines.
1400 ! All the others are subroutines.
1401 !
1402 !**End
1403 !
1404 ! Declare externals.
1405 ! Note: they are now internal
1406 !EXTERNAL DPREPJ, DSOLSY
1407 !KPP_REAL DUMACH, DVNORM
1408 !
1409 ! Declare all other variables.
1416 KPP_REAL ROWNS, &
1421 LOGICAL IHIT
1424 !-----------------------------------------------------------------------
1425 ! The following internal Common block contains
1426 ! (a) variables which are local to any subroutine but whose values must
1427 ! be preserved between calls to the routine ("own" variables), and
1428 ! (b) variables which are communicated between subroutines.
1429 ! The block DLS001 is declared in subroutines DLSODE, DINTDY, DSTODE,
1430 ! DPREPJ, and DSOLSY.
1431 ! Groups of variables are replaced by dummy arrays in the Common
1432 ! declarations in routines where those variables are not used.
1433 !-----------------------------------------------------------------------
1440 !
1442 !-----------------------------------------------------------------------
1443 ! Block A.
1444 ! This code block is executed on every call.
1445 ! It tests ISTATE and ITASK for legality and branches appropriately.
1446 ! If ISTATE .GT. 1 but the flag INIT shows that initialization has
1447 ! not yet been done, an error return occurs.
1448 ! If ISTATE = 1 and TOUT = T, return immediately.
1449 !-----------------------------------------------------------------------
1450 !
1451 !***FIRST EXECUTABLE STATEMENT DLSODE
1460 !-----------------------------------------------------------------------
1461 ! Block B.
1462 ! The next code block is executed for the initial call (ISTATE = 1),
1463 ! or for a continuation call with parameter changes (ISTATE = 3).
1464 ! It contains checking of all inputs and various initializations.
1465 !
1466 ! First check legality of the non-optional inputs NEQ, ITOL, IOPT,
1467 ! MF, ML, and MU.
1468 !-----------------------------------------------------------------------
1485 ! Next process and check the optional inputs. --------------------------
1488 MXSTEP = MXSTP0
1489 MXHNIL = MXHNL0
1513 !-----------------------------------------------------------------------
1514 ! Set work array pointers and check lengths LRW and LIW.
1515 ! Pointers to segments of RWORK and IWORK are named by prefixing L to
1516 ! the name of the segment. E.g., the segment YH starts at RWORK(LYH).
1517 ! Segments of RWORK (in order) are denoted YH, WM, EWT, SAVF, ACOR.
1518 !-----------------------------------------------------------------------
1537 ! Check RelTol and AbsTol for legality. ------------------------------------
1547 ! If ISTATE = 3, set flag to signal parameter changes to DSTODE. -------
1550 ! MAXORD was reduced below NQ. Copy YH(*,MAXORD+2) into SAVF. ---------
1553 ! Reload WM(1) = RWORK(LWM), since LWM may have changed. ---------------
1556 ! NEQ was reduced. Zero part of YH to avoid undefined references. -----
1563 !-----------------------------------------------------------------------
1564 ! Block C.
1565 ! The next block is for the initial call only (ISTATE = 1).
1566 ! It contains all remaining initializations, the initial call to F,
1567 ! and the calculation of the initial step size.
1568 ! The error weights in EWT are inverted after being loaded.
1569 !-----------------------------------------------------------------------
1571 TN = T
1589 ! Initial call to F. (LF0 points to YH(*,2).) -------------------------
1593 ! Load the initial value vector in YH. ---------------------------------
1596 ! Load and invert the EWT array. (H is temporarily set to 1.0.) -------
1603 !-----------------------------------------------------------------------
1604 ! The coding below computes the step size, H0, to be attempted on the
1605 ! first step, unless the user has supplied a value for this.
1606 ! First check that TOUT - T differs significantly from zero.
1607 ! A scalar tolerance quantity TOL is computed, as MAX(RelTol(I))
1608 ! if this is positive, or MAX(AbsTol(I)/ABS(Y(I))) otherwise, adjusted
1609 ! so as to be between 100*UROUND and 1.0E-3.
1610 ! Then the computed value H0 is given by..
1611 ! NEQ
1612 ! H0**2 = TOL / ( w0**-2 + (1/NEQ) * SUM ( f(i)/ywt(i) )**2 )
1613 ! 1
1614 ! where w0 = MAX ( ABS(T), ABS(TOUT) ),
1615 ! f(i) = i-th component of initial value of f,
1616 ! ywt(i) = EWT(i)/TOL (a weight for y(i)).
1617 ! The sign of H0 is inferred from the initial values of TOUT and T.
1618 !-----------------------------------------------------------------------
1641 ! Adjust H0 if necessary to meet HMAX bound. ---------------------------
1644 ! Load H with H0 and scale YH(*,2) by H0. ------------------------------
1645 H = H0
1649 !-----------------------------------------------------------------------
1650 ! Block D.
1651 ! The next code block is for continuation calls only (ISTATE = 2 or 3)
1652 ! and is to check stop conditions before taking a step.
1653 !-----------------------------------------------------------------------
1659 T = TOUT
1671 T = TOUT
1682 !-----------------------------------------------------------------------
1683 ! Block E.
1684 ! The next block is normally executed for all calls and contains
1685 ! the call to the one-step core integrator DSTODE.
1686 !
1687 ! This is a looping point for the integration steps.
1688 !
1689 ! First check for too many steps being taken, update EWT (if not at
1690 ! start of problem), check for too much accuracy being requested, and
1691 ! check for H below the roundoff level in T.
1692 !-----------------------------------------------------------------------
1719 !-----------------------------------------------------------------------
1720 ! CALL DSTODE(NEQ,Y,YH,NYH,YH,EWT,SAVF,ACOR,WM,IWM,F,JAC,DPREPJ,DSOLSY)
1721 !-----------------------------------------------------------------------
1725 !F, JAC, DPREPJ, DSOLSY)
1728 !-----------------------------------------------------------------------
1729 ! Block F.
1730 ! The following block handles the case of a successful return from the
1731 ! core integrator (KFLAG = 0). Test for stop conditions.
1732 !-----------------------------------------------------------------------
1735 ! ITASK = 1. If TOUT has been reached, interpolate. -------------------
1738 T = TOUT
1740 ! ITASK = 3. Jump to exit if TOUT was reached. ------------------------
1743 ! ITASK = 4. See if TOUT or TCRIT was reached. Adjust H if necessary.
1746 T = TOUT
1756 ! ITASK = 5. See if TCRIT was reached and jump to exit. ---------------
1759 !-----------------------------------------------------------------------
1760 ! Block G.
1761 ! The following block handles all successful returns from DLSODE.
1762 ! If ITASK .NE. 1, Y is loaded from YH and T is set accordingly.
1763 ! ISTATE is set to 2, and the optional outputs are loaded into the
1764 ! work arrays before returning.
1765 !-----------------------------------------------------------------------
1768 T = TN
1780 RETURN
1781 !-----------------------------------------------------------------------
1782 ! Block H.
1783 ! The following block handles all unsuccessful returns other than
1784 ! those for illegal input. First the error message routine is called.
1785 ! If there was an error test or convergence test failure, IMXER is set.
1786 ! Then Y is loaded from YH and T is set to TN. The optional outputs
1787 ! are loaded into the work arrays before returning.
1788 !-----------------------------------------------------------------------
1789 ! The maximum number of steps was taken before reaching TOUT. ----------
1796 ! EWT(I) .LE. 0.0 for some I (not at start of problem). ----------------
1802 ! Too much accuracy requested for machine precision. -------------------
1810 ! KFLAG = -1. Error test failed repeatedly or with ABS(H) = HMIN. -----
1817 ! KFLAG = -2. Convergence failed repeatedly or with ABS(H) = HMIN. ----
1825 ! Compute IMXER if relevant. -------------------------------------------
1831 BIG = SIZE
1832 IMXER = I
1835 ! Set Y vector, T, and optional outputs. -------------------------------
1838 T = TN
1847 RETURN
1848 !-----------------------------------------------------------------------
1849 ! Block I.
1850 ! The following block handles all error returns due to illegal input
1851 ! (ISTATE = -3), as detected before calling the core integrator.
1852 ! First the error message routine is called. If the illegal input
1853 ! is a negative ISTATE, the run is aborted (apparent infinite loop).
1854 !-----------------------------------------------------------------------
1948 !
1950 RETURN
1951 !
1954 RETURN
1955 !----------------------- END OF SUBROUTINE DLSODE ----------------------
1956 !END SUBROUTINE DLSODE
1957 CONTAINS
1960 !DECK DUMACH
1962 !***BEGIN PROLOGUE DUMACH
1963 !***PURPOSE Compute the unit roundoff of the machine.
1964 !***CATEGORY R1
1965 !***TYPE KPP_REAL (RUMACH-S, DUMACH-D)
1966 !***KEYWORDS MACHINE CONSTANTS
1967 !***AUTHOR Hindmarsh, Alan C., (LLNL)
1968 !***DESCRIPTION
1969 ! *Usage:
1970 ! KPP_REAL A, DUMACH
1971 ! A = DUMACH()
1972 !
1973 ! *Function Return Values:
1974 ! A : the unit roundoff of the machine.
1975 !
1976 ! *Description:
1977 ! The unit roundoff is defined as the smallest positive machine
1978 ! number u such that 1.0 + u .ne. 1.0. This is computed by DUMACH
1979 ! in a machine-independent manner.
1980 !
1981 !***REFERENCES (NONE)
1982 !***ROUTINES CALLED DUMSUM
1983 !***REVISION HISTORY (YYYYMMDD)
1984 ! 19930216 DATE WRITTEN
1985 ! 19930818 Added SLATEC-format prologue. (FNF)
1986 ! 20030707 Added DUMSUM to force normal storage of COMP. (ACH)
1987 !***END PROLOGUE DUMACH
1988 !
1989 KPP_REAL U, COMP
1990 !***FIRST EXECUTABLE STATEMENT DUMACH
1996 RETURN
1997 !----------------------- End of Function DUMACH ------------------------
1998 END FUNCTION DUMACH
2001 ! Routine to force normal storing of A + B, for DUMACH.
2004 RETURN
2005 END SUBROUTINE DUMSUM
2006 !DECK DCFODE
2008 !***BEGIN PROLOGUE DCFODE
2009 !***SUBSIDIARY
2010 !***PURPOSE Set ODE integrator coefficients.
2011 !***TYPE KPP_REAL (SCFODE-S, DCFODE-D)
2012 !***AUTHOR Hindmarsh, Alan C., (LLNL)
2013 !***DESCRIPTION
2014 !
2015 ! DCFODE is called by the integrator routine to set coefficients
2016 ! needed there. The coefficients for the current method, as
2017 ! given by the value of METH, are set for all orders and saved.
2018 ! The maximum order assumed here is 12 if METH = 1 and 5 if METH = 2.
2019 ! (A smaller value of the maximum order is also allowed.)
2020 ! DCFODE is called once at the beginning of the problem,
2021 ! and is not called again unless and until METH is changed.
2022 !
2023 ! The ELCO array contains the basic method coefficients.
2024 ! The coefficients el(i), 1 .le. i .le. nq+1, for the method of
2025 ! order nq are stored in ELCO(i,nq). They are given by a genetrating
2026 ! polynomial, i.e.,
2027 ! l(x) = el(1) + el(2)*x + ... + el(nq+1)*x**nq.
2028 ! For the implicit Adams methods, l(x) is given by
2029 ! dl/dx = (x+1)*(x+2)*...*(x+nq-1)/factorial(nq-1), l(-1) = 0.
2030 ! For the BDF methods, l(x) is given by
2031 ! l(x) = (x+1)*(x+2)* ... *(x+nq)/K,
2032 ! where K = factorial(nq)*(1 + 1/2 + ... + 1/nq).
2033 !
2034 ! The TESCO array contains test constants used for the
2035 ! local error test and the selection of step size and/or order.
2036 ! At order nq, TESCO(k,nq) is used for the selection of step
2037 ! size at order nq - 1 if k = 1, at order nq if k = 2, and at order
2038 ! nq + 1 if k = 3.
2039 !
2040 !***SEE ALSO DLSODE
2041 !***ROUTINES CALLED (NONE)
2042 !***REVISION HISTORY (YYMMDD)
2043 ! 791129 DATE WRITTEN
2044 ! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
2045 ! 890503 Minor cosmetic changes. (FNF)
2046 ! 930809 Renamed to allow single/double precision versions. (ACH)
2047 !***END PROLOGUE DCFODE
2048 !**End
2049 INTEGER METH
2053 !
2054 !***FIRST EXECUTABLE STATEMENT DCFODE
2056 !
2066 !-----------------------------------------------------------------------
2067 ! The PC array will contain the coefficients of the polynomial
2068 ! p(x) = (x+1)*(x+2)*...*(x+nq-1).
2069 ! Initially, p(x) = 1.
2070 !-----------------------------------------------------------------------
2071 RQ1FAC = RQFAC
2074 FNQM1 = NQM1
2076 ! Form coefficients of p(x)*(x+nq-1). ----------------------------------
2082 ! Compute integral, -1 to 0, of p(x) and x*p(x). -----------------------
2087 TSIGN = -TSIGN
2090 ! Store coefficients in ELCO and TESCO. --------------------------------
2101 RETURN
2102 !
2106 !-----------------------------------------------------------------------
2107 ! The PC array will contain the coefficients of the polynomial
2108 ! p(x) = (x+1)*(x+2)*...*(x+nq).
2109 ! Initially, p(x) = 1.
2110 !-----------------------------------------------------------------------
2111 FNQ = NQ
2113 ! Form coefficients of p(x)*(x+nq). ------------------------------------
2119 ! Store coefficients in ELCO and TESCO. --------------------------------
2128 RETURN
2129 !----------------------- END OF SUBROUTINE DCFODE ----------------------
2130 END SUBROUTINE DCFODE
2131 !DECK DINTDY
2133 !***BEGIN PROLOGUE DINTDY
2134 !***SUBSIDIARY
2135 !***PURPOSE Interpolate solution derivatives.
2136 !***TYPE KPP_REAL (SINTDY-S, DINTDY-D)
2137 !***AUTHOR Hindmarsh, Alan C., (LLNL)
2138 !***DESCRIPTION
2139 !
2140 ! DINTDY computes interpolated values of the K-th derivative of the
2141 ! dependent variable vector y, and stores it in DKY. This routine
2142 ! is called within the package with K = 0 and T = TOUT, but may
2143 ! also be called by the user for any K up to the current order.
2144 ! (See detailed instructions in the usage documentation.)
2145 !
2146 ! The computed values in DKY are gotten by interpolation using the
2147 ! Nordsieck history array YH. This array corresponds uniquely to a
2148 ! vector-valued polynomial of degree NQCUR or less, and DKY is set
2149 ! to the K-th derivative of this polynomial at T.
2150 ! The formula for DKY is:
2151 ! q
2152 ! DKY(i) = sum c(j,K) * (T - tn)**(j-K) * h**(-j) * YH(i,j+1)
2153 ! j=K
2154 ! where c(j,K) = j*(j-1)*...*(j-K+1), q = NQCUR, tn = TCUR, h = HCUR.
2155 ! The quantities nq = NQCUR, l = nq+1, N = NEQ, tn, and h are
2156 ! communicated by COMMON. The above sum is done in reverse order.
2157 ! IFLAG is returned negative if either K or T is out of bounds.
2158 !
2159 !***SEE ALSO DLSODE
2160 !***ROUTINES CALLED XERRWD
2161 !***COMMON BLOCKS DLS001
2162 !***REVISION HISTORY (YYMMDD)
2163 ! 791129 DATE WRITTEN
2164 ! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
2165 ! 890503 Minor cosmetic changes. (FNF)
2166 ! 930809 Renamed to allow single/double precision versions. (ACH)
2167 ! 010418 Reduced size of Common block /DLS001/. (ACH)
2168 ! 031105 Restored 'own' variables to Common block /DLS001/, to
2169 ! enable interrupt/restart feature. (ACH)
2170 ! 050427 Corrected roundoff decrement in TP. (ACH)
2171 !***END PROLOGUE DINTDY
2172 !**End
2179 KPP_REAL ROWNS, &
2190 !
2191 !***FIRST EXECUTABLE STATEMENT DINTDY
2196 !
2224 RETURN
2225 !
2229 RETURN
2235 RETURN
2236 !----------------------- END OF SUBROUTINE DINTDY ----------------------
2237 END SUBROUTINE DINTDY
2238 !DECK DPREPJ
2240 !***BEGIN PROLOGUE DPREPJ
2241 !***SUBSIDIARY
2242 !***PURPOSE Compute and process Newton iteration matrix.
2243 !***TYPE KPP_REAL (SPREPJ-S, DPREPJ-D)
2244 !***AUTHOR Hindmarsh, Alan C., (LLNL)
2245 !***DESCRIPTION
2246 !
2247 ! DPREPJ is called by DSTODE to compute and process the matrix
2248 ! P = I - h*el(1)*J , where J is an approximation to the Jacobian.
2249 ! Here J is computed by the user-supplied routine JAC if
2250 ! MITER = 1 or 4, or by finite differencing if MITER = 2, 3, or 5.
2251 ! If MITER = 3, a diagonal approximation to J is used.
2252 ! J is stored in WM and replaced by P. If MITER .ne. 3, P is then
2253 ! subjected to LU decomposition in preparation for later solution
2254 ! of linear systems with P as coefficient matrix. This is done
2255 ! by DGEFA if MITER = 1 or 2, and by DGBFA if MITER = 4 or 5.
2256 !
2257 ! In addition to variables described in DSTODE and DLSODE prologues,
2258 ! communication with DPREPJ uses the following:
2259 ! Y = array containing predicted values on entry.
2260 ! FTEM = work array of length N (ACOR in DSTODE).
2261 ! SAVF = array containing f evaluated at predicted y.
2262 ! WM = real work space for matrices. On output it contains the
2263 ! inverse diagonal matrix if MITER = 3 and the LU decomposition
2264 ! of P if MITER is 1, 2 , 4, or 5.
2265 ! Storage of matrix elements starts at WM(3).
2266 ! WM also contains the following matrix-related data:
2267 ! WM(1) = SQRT(UROUND), used in numerical Jacobian increments.
2268 ! WM(2) = H*EL0, saved for later use if MITER = 3.
2269 ! IWM = integer work space containing pivot information, starting at
2270 ! IWM(21), if MITER is 1, 2, 4, or 5. IWM also contains band
2271 ! parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5.
2272 ! EL0 = EL(1) (input).
2273 ! IERPJ = output error flag, = 0 if no trouble, .gt. 0 if
2274 ! P matrix found to be singular.
2275 ! JCUR = output flag = 1 to indicate that the Jacobian matrix
2276 ! (or approximation) is now current.
2277 ! This routine also uses the COMMON variables EL0, H, TN, UROUND,
2278 ! MITER, N, NFE, and NJE.
2279 !
2280 !***SEE ALSO DLSODE
2281 !***ROUTINES CALLED DGBFA, DGEFA, DVNORM
2282 !***COMMON BLOCKS DLS001
2283 !***REVISION HISTORY (YYMMDD)
2284 ! 791129 DATE WRITTEN
2285 ! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
2286 ! 890504 Minor cosmetic changes. (FNF)
2287 ! 930809 Renamed to allow single/double precision versions. (ACH)
2288 ! 010418 Reduced size of Common block /DLS001/. (ACH)
2289 ! 031105 Restored 'own' variables to Common block /DLS001/, to
2290 ! enable interrupt/restart feature. (ACH)
2291 !***END PROLOGUE DPREPJ
2292 !**End
2300 KPP_REAL ROWNS, &
2311 !KPP_REAL DVNORM
2312 !
2313 !***FIRST EXECUTABLE STATEMENT DPREPJ
2318 CON = -HL0
2320 #ifdef FULL_ALGEBRA
2324 END DO
2328 END DO
2329 ! Add identity matrix
2335 END DO
2336 ! Do LU decomposition on P
2338 #else
2342 END DO
2343 ! Add identity matrix
2347 END DO
2348 ! Do LU decomposition on P
2350 #endif
2352 RETURN
2353 !----------------------- END OF SUBROUTINE DPREPJ ----------------------
2354 END SUBROUTINE DPREPJ
2355 !DECK DSOLSY
2357 !***BEGIN PROLOGUE DSOLSY
2358 !***SUBSIDIARY
2359 !***PURPOSE ODEPACK linear system solver.
2360 !***TYPE KPP_REAL (SSOLSY-S, DSOLSY-D)
2361 !***AUTHOR Hindmarsh, Alan C., (LLNL)
2362 !***DESCRIPTION
2363 !
2364 ! This routine manages the solution of the linear system arising from
2365 ! a chord iteration. It is called if MITER .ne. 0.
2366 ! If MITER is 1 or 2, it calls DGESL to accomplish this.
2367 ! If MITER = 3 it updates the coefficient h*EL0 in the diagonal
2368 ! matrix, and then computes the solution.
2369 ! If MITER is 4 or 5, it calls DGBSL.
2370 ! Communication with DSOLSY uses the following variables:
2371 ! WM = real work space containing the inverse diagonal matrix if
2372 ! MITER = 3 and the LU decomposition of the matrix otherwise.
2373 ! Storage of matrix elements starts at WM(3).
2374 ! WM also contains the following matrix-related data:
2375 ! WM(1) = SQRT(UROUND) (not used here),
2376 ! WM(2) = HL0, the previous value of h*EL0, used if MITER = 3.
2377 ! IWM = integer work space containing pivot information, starting at
2378 ! IWM(21), if MITER is 1, 2, 4, or 5. IWM also contains band
2379 ! parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5.
2380 ! X = the right-hand side vector on input, and the solution vector
2381 ! on output, of length N.
2382 ! TEM = vector of work space of length N, not used in this version.
2383 ! IERSL = output flag (in COMMON). IERSL = 0 if no trouble occurred.
2384 ! IERSL = 1 if a singular matrix arose with MITER = 3.
2385 ! This routine also uses the COMMON variables EL0, H, MITER, and N.
2386 !
2387 !***SEE ALSO DLSODE
2388 !***ROUTINES CALLED DGBSL, DGESL
2389 !***COMMON BLOCKS DLS001
2390 !***REVISION HISTORY (YYMMDD)
2391 ! 791129 DATE WRITTEN
2392 ! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
2393 ! 890503 Minor cosmetic changes. (FNF)
2394 ! 930809 Renamed to allow single/double precision versions. (ACH)
2395 ! 010418 Reduced size of Common block /DLS001/. (ACH)
2396 ! 031105 Restored 'own' variables to Common block /DLS001/, to
2397 ! enable interrupt/restart feature. (ACH)
2398 !***END PROLOGUE DSOLSY
2399 !**End
2406 KPP_REAL ROWNS, &
2416 !
2417 !***FIRST EXECUTABLE STATEMENT DSOLSY
2419 #ifdef FULL_ALGEBRA
2421 #else
2423 #endif
2424 RETURN
2425 !----------------------- END OF SUBROUTINE DSOLSY ----------------------
2426 END SUBROUTINE DSOLSY
2427 !DECK DSRCOM
2429 !***BEGIN PROLOGUE DSRCOM
2430 !***SUBSIDIARY
2431 !***PURPOSE Save/restore ODEPACK COMMON blocks.
2432 !***TYPE KPP_REAL (SSRCOM-S, DSRCOM-D)
2433 !***AUTHOR Hindmarsh, Alan C., (LLNL)
2434 !***DESCRIPTION
2435 !
2436 ! This routine saves or restores (depending on JOB) the contents of
2437 ! the COMMON block DLS001, which is used internally
2438 ! by one or more ODEPACK solvers.
2439 !
2440 ! RSAV = real array of length 218 or more.
2441 ! ISAV = integer array of length 37 or more.
2442 ! JOB = flag indicating to save or restore the COMMON blocks:
2443 ! JOB = 1 if COMMON is to be saved (written to RSAV/ISAV)
2444 ! JOB = 2 if COMMON is to be restored (read from RSAV/ISAV)
2445 ! A call with JOB = 2 presumes a prior call with JOB = 1.
2446 !
2447 !***SEE ALSO DLSODE
2448 !***ROUTINES CALLED (NONE)
2449 !***COMMON BLOCKS DLS001
2450 !***REVISION HISTORY (YYMMDD)
2451 ! 791129 DATE WRITTEN
2452 ! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
2453 ! 890503 Minor cosmetic changes. (FNF)
2454 ! 921116 Deleted treatment of block /EH0001/. (ACH)
2455 ! 930801 Reduced Common block length by 2. (ACH)
2456 ! 930809 Renamed to allow single/double precision versions. (ACH)
2457 ! 010418 Reduced Common block length by 209+12. (ACH)
2458 ! 031105 Restored 'own' variables to Common block /DLS001/, to
2459 ! enable interrupt/restart feature. (ACH)
2460 ! 031112 Added SAVE statement for data-loaded constants.
2461 !***END PROLOGUE DSRCOM
2462 !**End
2464 INTEGER ILS
2470 !
2471 !***FIRST EXECUTABLE STATEMENT DSRCOM
2473 !
2478 RETURN
2479 !
2485 RETURN
2486 !----------------------- END OF SUBROUTINE DSRCOM ----------------------
2487 END SUBROUTINE DSRCOM
2488 !DECK DSTODE
2491 !WM, IWM, F, JAC, PJAC, SLVS)
2492 !***BEGIN PROLOGUE DSTODE
2493 !***SUBSIDIARY
2494 !***PURPOSE Performs one step of an ODEPACK integration.
2495 !***TYPE KPP_REAL (SSTODE-S, DSTODE-D)
2496 !***AUTHOR Hindmarsh, Alan C., (LLNL)
2497 !***DESCRIPTION
2498 !
2499 ! DSTODE performs one step of the integration of an initial value
2500 ! problem for a system of ordinary differential equations.
2501 ! Note: DSTODE is independent of the value of the iteration method
2502 ! indicator MITER, when this is .ne. 0, and hence is independent
2503 ! of the type of chord method used, or the Jacobian structure.
2504 ! Communication with DSTODE is done with the following variables:
2505 !
2506 ! NEQ = integer array containing problem size in NEQ, and
2507 ! passed as the NEQ argument in all calls to F and JAC.
2508 ! Y = an array of length .ge. N used as the Y argument in
2509 ! all calls to F and JAC.
2510 ! YH = an NYH by LMAX array containing the dependent variables
2511 ! and their approximate scaled derivatives, where
2512 ! LMAX = MAXORD + 1. YH(i,j+1) contains the approximate
2513 ! j-th derivative of y(i), scaled by h**j/factorial(j)
2514 ! (j = 0,1,...,NQ). on entry for the first step, the first
2515 ! two columns of YH must be set from the initial values.
2516 ! NYH = a constant integer .ge. N, the first dimension of YH.
2517 ! YH1 = a one-dimensional array occupying the same space as YH.
2518 ! EWT = an array of length N containing multiplicative weights
2519 ! for local error measurements. Local errors in Y(i) are
2520 ! compared to 1.0/EWT(i) in various error tests.
2521 ! SAVF = an array of working storage, of length N.
2522 ! Also used for input of YH(*,MAXORD+2) when JSTART = -1
2523 ! and MAXORD .lt. the current order NQ.
2524 ! ACOR = a work array of length N, used for the accumulated
2525 ! corrections. On a successful return, ACOR(i) contains
2526 ! the estimated one-step local error in Y(i).
2527 ! WM,IWM = real and integer work arrays associated with matrix
2528 ! operations in chord iteration (MITER .ne. 0).
2529 ! PJAC = name of routine to evaluate and preprocess Jacobian matrix
2530 ! and P = I - h*el0*JAC, if a chord method is being used.
2531 ! SLVS = name of routine to solve linear system in chord iteration.
2532 ! CCMAX = maximum relative change in h*el0 before PJAC is called.
2533 ! H = the step size to be attempted on the next step.
2534 ! H is altered by the error control algorithm during the
2535 ! problem. H can be either positive or negative, but its
2536 ! sign must remain constant throughout the problem.
2537 ! HMIN = the minimum absolute value of the step size h to be used.
2538 ! HMXI = inverse of the maximum absolute value of h to be used.
2539 ! HMXI = 0.0 is allowed and corresponds to an infinite hmax.
2540 ! HMIN and HMXI may be changed at any time, but will not
2541 ! take effect until the next change of h is considered.
2542 ! TN = the independent variable. TN is updated on each step taken.
2543 ! JSTART = an integer used for input only, with the following
2544 ! values and meanings:
2545 ! 0 perform the first step.
2546 ! .gt.0 take a new step continuing from the last.
2547 ! -1 take the next step with a new value of H, MAXORD,
2548 ! N, METH, MITER, and/or matrix parameters.
2549 ! -2 take the next step with a new value of H,
2550 ! but with other inputs unchanged.
2551 ! On return, JSTART is set to 1 to facilitate continuation.
2552 ! KFLAG = a completion code with the following meanings:
2553 ! 0 the step was succesful.
2554 ! -1 the requested error could not be achieved.
2555 ! -2 corrector convergence could not be achieved.
2556 ! -3 fatal error in PJAC or SLVS.
2557 ! A return with KFLAG = -1 or -2 means either
2558 ! abs(H) = HMIN or 10 consecutive failures occurred.
2559 ! On a return with KFLAG negative, the values of TN and
2560 ! the YH array are as of the beginning of the last
2561 ! step, and H is the last step size attempted.
2562 ! MAXORD = the maximum order of integration method to be allowed.
2563 ! MAXCOR = the maximum number of corrector iterations allowed.
2564 ! MSBP = maximum number of steps between PJAC calls (MITER .gt. 0).
2565 ! MXNCF = maximum number of convergence failures allowed.
2566 ! METH/MITER = the method flags. See description in driver.
2567 ! N = the number of first-order differential equations.
2568 ! The values of CCMAX, H, HMIN, HMXI, TN, JSTART, KFLAG, MAXORD,
2569 ! MAXCOR, MSBP, MXNCF, METH, MITER, and N are communicated via COMMON.
2570 !
2571 !***SEE ALSO DLSODE
2572 !***ROUTINES CALLED DCFODE, DVNORM
2573 !***COMMON BLOCKS DLS001
2574 !***REVISION HISTORY (YYMMDD)
2575 ! 791129 DATE WRITTEN
2576 ! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
2577 ! 890503 Minor cosmetic changes. (FNF)
2578 ! 930809 Renamed to allow single/double precision versions. (ACH)
2579 ! 010418 Reduced size of Common block /DLS001/. (ACH)
2580 ! 031105 Restored 'own' variables to Common block /DLS001/, to
2581 ! enable interrupt/restart feature. (ACH)
2582 !***END PROLOGUE DSTODE
2583 !**End
2597 !KPP_REAL DVNORM
2605 !
2606 !***FIRST EXECUTABLE STATEMENT DSTODE
2608 TOLD = TN
2618 !-----------------------------------------------------------------------
2619 ! On the first call, the order is set to 1, and other variables are
2620 ! initialized. RMAX is the maximum ratio by which H can be increased
2621 ! in a single step. It is initially 1.E4 to compensate for the small
2622 ! initial H, but then is normally equal to 10. If a failure
2623 ! occurs (in corrector convergence or error test), RMAX is set to 2
2624 ! for the next increase.
2625 !-----------------------------------------------------------------------
2634 HOLD = H
2635 MEO = METH
2637 IPUP = MITER
2640 !-----------------------------------------------------------------------
2641 ! The following block handles preliminaries needed when JSTART = -1.
2642 ! IPUP is set to MITER to force a matrix update.
2643 ! If an order increase is about to be considered (IALTH = 1),
2644 ! IALTH is reset to 2 to postpone consideration one more step.
2645 ! If the caller has changed METH, DCFODE is called to reset
2646 ! the coefficients of the method.
2647 ! If the caller has changed MAXORD to a value less than the current
2648 ! order NQ, NQ is reduced to MAXORD, and a new H chosen accordingly.
2649 ! If H is to be changed, YH must be rescaled.
2650 ! If H or METH is being changed, IALTH is reset to L = NQ + 1
2651 ! to prevent further changes in H for that many steps.
2652 !-----------------------------------------------------------------------
2658 MEO = METH
2660 IALTH = L
2665 L = LMAX
2679 H = HOLD
2681 !-----------------------------------------------------------------------
2682 ! DCFODE is called to get all the integration coefficients for the
2683 ! current METH. Then the EL vector and related constants are reset
2684 ! whenever the order NQ is changed, or at the start of the problem.
2685 !-----------------------------------------------------------------------
2694 !-----------------------------------------------------------------------
2695 ! If H is being changed, the H ratio RH is checked against
2696 ! RMAX, HMIN, and HMXI, and the YH array rescaled. IALTH is set to
2697 ! L = NQ + 1 to prevent a change of H for that many steps, unless
2698 ! forced by a convergence or error test failure.
2699 !-----------------------------------------------------------------------
2702 H = HOLD
2715 IALTH = L
2717 !-----------------------------------------------------------------------
2718 ! This section computes the predicted values by effectively
2719 ! multiplying the YH array by the Pascal Triangle matrix.
2720 ! RC is the ratio of new to old values of the coefficient H*EL(1).
2721 ! When RC differs from 1 by more than CCMAX, IPUP is set to MITER
2722 ! to force PJAC to be called, if a Jacobian is involved.
2723 ! In any case, PJAC is called at least every MSBP steps.
2724 !-----------------------------------------------------------------------
2731 !dir$ ivdep
2735 !-----------------------------------------------------------------------
2736 ! Up to MAXCOR corrector iterations are taken. A convergence test is
2737 ! made on the R.M.S. norm of each correction, weighted by the error
2738 ! weight vector EWT. The sum of the corrections is accumulated in the
2739 ! vector ACOR(i). The YH array is not altered in the corrector loop.
2740 !-----------------------------------------------------------------------
2747 !-----------------------------------------------------------------------
2748 ! If indicated, the matrix P = I - h*el(1)*J is reevaluated and
2749 ! preprocessed before starting the corrector iteration. IPUP is set
2750 ! to 0 as an indicator that this has been done.
2751 !-----------------------------------------------------------------------
2752 !CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC)
2756 NSLP = NST
2762 !-----------------------------------------------------------------------
2763 ! In the case of functional iteration, update Y directly from
2764 ! the result of the last function evaluation.
2765 !-----------------------------------------------------------------------
2774 !-----------------------------------------------------------------------
2775 ! In the case of the chord method, compute the corrector error,
2776 ! and solve the linear system with that as right-hand side and
2777 ! P as coefficient matrix.
2778 !-----------------------------------------------------------------------
2781 !CALL SLVS (WM, IWM, Y, SAVF)
2789 !-----------------------------------------------------------------------
2790 ! Test for convergence. If M.gt.0, an estimate of the convergence
2791 ! rate constant is stored in CRATE, and this is used in the test.
2792 !-----------------------------------------------------------------------
2799 DELP = DEL
2803 !-----------------------------------------------------------------------
2804 ! The corrector iteration failed to converge.
2805 ! If MITER .ne. 0 and the Jacobian is out of date, PJAC is called for
2806 ! the next try. Otherwise the YH array is retracted to its values
2807 ! before prediction, and H is reduced, if possible. If H cannot be
2808 ! reduced or MXNCF failures have occurred, exit with KFLAG = -2.
2809 !-----------------------------------------------------------------------
2812 IPUP = MITER
2817 TN = TOLD
2821 !dir$ ivdep
2829 IPUP = MITER
2832 !-----------------------------------------------------------------------
2833 ! The corrector has converged. JCUR is set to 0
2834 ! to signal that the Jacobian involved may need updating later.
2835 ! The local error test is made and control passes to statement 500
2836 ! if it fails.
2837 !-----------------------------------------------------------------------
2842 !-----------------------------------------------------------------------
2843 ! After a successful step, update the YH array.
2844 ! Consider changing H if IALTH = 1. Otherwise decrease IALTH by 1.
2845 ! If IALTH is then 1 and NQ .lt. MAXORD, then ACOR is saved for
2846 ! use in a possible order increase on the next step.
2847 ! If a change in H is considered, an increase or decrease in order
2848 ! by one is considered also. A change in H is made only if it is by a
2849 ! factor of at least 1.1. If not, IALTH is set to 3 to prevent
2850 ! testing for that many steps.
2851 !-----------------------------------------------------------------------
2855 HU = H
2856 NQU = NQ
2867 !-----------------------------------------------------------------------
2868 ! The error test failed. KFLAG keeps track of multiple failures.
2869 ! Restore TN and the YH array to their previous values, and prepare
2870 ! to try the step again. Compute the optimum step size for this or
2871 ! one lower order. After 2 or more failures, H is forced to decrease
2872 ! by a factor of 0.2 or less.
2873 !-----------------------------------------------------------------------
2875 TN = TOLD
2879 !dir$ ivdep
2889 !-----------------------------------------------------------------------
2890 ! Regardless of the success or failure of the step, factors
2891 ! RHDN, RHSM, and RHUP are computed, by which H could be multiplied
2892 ! at order NQ - 1, order NQ, or order NQ + 1, respectively.
2893 ! In the case of failure, RHUP = 0.0 to avoid an order increase.
2894 ! The largest of these is determined and the new order chosen
2895 ! accordingly. If the order is to be increased, we compute one
2896 ! additional scaled derivative.
2897 !-----------------------------------------------------------------------
2916 NEWQ = NQ
2917 RH = RHSM
2920 RH = RHDN
2924 RH = RHUP
2934 !-----------------------------------------------------------------------
2935 ! If there is a change of order, reset NQ, l, and the coefficients.
2936 ! In any case H is reset according to RH and the YH array is rescaled.
2937 ! Then exit from 690 if the step was OK, or redo the step otherwise.
2938 !-----------------------------------------------------------------------
2944 !-----------------------------------------------------------------------
2945 ! Control reaches this section if 3 or more failures have occured.
2946 ! If 10 failures have occurred, exit with KFLAG = -1.
2947 ! It is assumed that the derivatives that have accumulated in the
2948 ! YH array have errors of the wrong order. Hence the first
2949 ! derivative is recomputed, and the order is set to 1. Then
2950 ! H is reduced by a factor of 10, and the step is retried,
2951 ! until it succeeds or H reaches HMIN.
2952 !-----------------------------------------------------------------------
2963 IPUP = MITER
2970 !-----------------------------------------------------------------------
2971 ! All returns are made through this section. H is saved in HOLD
2972 ! to allow the caller to change H on the next step.
2973 !-----------------------------------------------------------------------
2986 RETURN
2987 !----------------------- END OF SUBROUTINE DSTODE ----------------------
2988 END SUBROUTINE DSTODE
2989 !DECK DEWSET
2991 !***BEGIN PROLOGUE DEWSET
2992 !***SUBSIDIARY
2993 !***PURPOSE Set error weight vector.
2994 !***TYPE KPP_REAL (SEWSET-S, DEWSET-D)
2995 !***AUTHOR Hindmarsh, Alan C., (LLNL)
2996 !***DESCRIPTION
2997 !
2998 ! This subroutine sets the error weight vector EWT according to
2999 ! EWT(i) = RelTol(i)*ABS(YCUR(i)) + AbsTol(i), i = 1,...,N,
3000 ! with the subscript on RelTol and/or AbsTol possibly replaced by 1 above,
3001 ! depending on the value of ITOL.
3002 !
3003 !***SEE ALSO DLSODE
3004 !***ROUTINES CALLED (NONE)
3005 !***REVISION HISTORY (YYMMDD)
3006 ! 791129 DATE WRITTEN
3007 ! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
3008 ! 890503 Minor cosmetic changes. (FNF)
3009 ! 930809 Renamed to allow single/double precision versions. (ACH)
3010 !***END PROLOGUE DEWSET
3011 !**End
3013 INTEGER I
3015 !
3016 !***FIRST EXECUTABLE STATEMENT DEWSET
3021 RETURN
3025 RETURN
3029 RETURN
3033 RETURN
3034 !----------------------- END OF SUBROUTINE DEWSET ----------------------
3035 END SUBROUTINE DEWSET
3036 !DECK DVNORM
3038 !***BEGIN PROLOGUE DVNORM
3039 !***SUBSIDIARY
3040 !***PURPOSE Weighted root-mean-square vector norm.
3041 !***TYPE KPP_REAL (SVNORM-S, DVNORM-D)
3042 !***AUTHOR Hindmarsh, Alan C., (LLNL)
3043 !***DESCRIPTION
3044 !
3045 ! This function routine computes the weighted root-mean-square norm
3046 ! of the vector of length N contained in the array V, with weights
3047 ! contained in the array W of length N:
3048 ! DVNORM = SQRT( (1/N) * SUM( V(i)*W(i) )**2 )
3049 !
3050 !***SEE ALSO DLSODE
3051 !***ROUTINES CALLED (NONE)
3052 !***REVISION HISTORY (YYMMDD)
3053 ! 791129 DATE WRITTEN
3054 ! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
3055 ! 890503 Minor cosmetic changes. (FNF)
3056 ! 930809 Renamed to allow single/double precision versions. (ACH)
3057 !***END PROLOGUE DVNORM
3058 !**End
3061 !
3062 !***FIRST EXECUTABLE STATEMENT DVNORM
3067 RETURN
3068 !----------------------- END OF FUNCTION DVNORM ------------------------
3069 END FUNCTION DVNORM
3070 !DECK XERRWD
3072 !***BEGIN PROLOGUE XERRWD
3073 !***SUBSIDIARY
3074 !***PURPOSE Write error message with values.
3075 !***CATEGORY R3C
3076 !***TYPE KPP_REAL (XERRWV-S, XERRWD-D)
3077 !***AUTHOR Hindmarsh, Alan C., (LLNL)
3078 !***DESCRIPTION
3079 !
3080 ! Subroutines XERRWD, XSETF, XSETUN, and the function routine IXSAV,
3081 ! as given here, constitute a simplified version of the SLATEC error
3082 ! handling package.
3083 !
3084 ! All arguments are input arguments.
3085 !
3086 ! MSG = The message (character array).
3087 ! NMES = The length of MSG (number of characters).
3088 ! NERR = The error number (not used).
3089 ! LEVEL = The error level..
3090 ! 0 or 1 means recoverable (control returns to caller).
3091 ! 2 means fatal (run is aborted--see note below).
3092 ! NI = Number of integers (0, 1, or 2) to be printed with message.
3093 ! I1,I2 = Integers to be printed, depending on NI.
3094 ! NR = Number of reals (0, 1, or 2) to be printed with message.
3095 ! R1,R2 = Reals to be printed, depending on NR.
3096 !
3097 ! Note.. this routine is machine-dependent and specialized for use
3098 ! in limited context, in the following ways..
3099 ! 1. The argument MSG is assumed to be of type CHARACTER, and
3100 ! the message is printed with a format of (1X,A).
3101 ! 2. The message is assumed to take only one line.
3102 ! Multi-line messages are generated by repeated calls.
3103 ! 3. If LEVEL = 2, control passes to the statement STOP
3104 ! to abort the run. This statement may be machine-dependent.
3105 ! 4. R1 and R2 are assumed to be in double precision and are printed
3106 ! in D21.13 format.
3107 !
3108 !***ROUTINES CALLED IXSAV
3109 !***REVISION HISTORY (YYMMDD)
3110 ! 920831 DATE WRITTEN
3111 ! 921118 Replaced MFLGSV/LUNSAV by IXSAV. (ACH)
3112 ! 930329 Modified prologue to SLATEC format. (FNF)
3113 ! 930407 Changed MSG from CHARACTER*1 array to variable. (FNF)
3114 ! 930922 Minor cosmetic change. (FNF)
3115 !***END PROLOGUE XERRWD
3116 !
3117 !*Internal Notes:
3118 !
3119 ! For a different default logical unit number, IXSAV (or a subsidiary
3120 ! routine that it calls) will need to be modified.
3121 ! For a different run-abort command, change the statement following
3122 ! statement 100 at the end.
3123 !-----------------------------------------------------------------------
3124 ! Subroutines called by XERRWD.. None
3125 ! Function routine called by XERRWD.. IXSAV
3126 !-----------------------------------------------------------------------
3127 !**End
3128 !
3129 ! Declare arguments.
3130 !
3131 KPP_REAL R1, R2
3134 !
3135 ! Declare local variables.
3136 !
3138 !
3139 ! Get logical unit number and message print flag.
3140 !
3141 !***FIRST EXECUTABLE STATEMENT XERRWD
3145 !
3146 ! Write the message.
3147 !
3158 !
3159 ! Abort the run if LEVEL = 2.
3160 !
3162 STOP
3163 !----------------------- End of Subroutine XERRWD ----------------------
3164 END SUBROUTINE XERRWD
3165 !DECK XSETF
3167 !***BEGIN PROLOGUE XSETF
3168 !***PURPOSE Reset the error print control flag.
3169 !***CATEGORY R3A
3170 !***TYPE ALL (XSETF-A)
3171 !***KEYWORDS ERROR CONTROL
3172 !***AUTHOR Hindmarsh, Alan C., (LLNL)
3173 !***DESCRIPTION
3174 !
3175 ! XSETF sets the error print control flag to MFLAG:
3176 ! MFLAG=1 means print all messages (the default).
3177 ! MFLAG=0 means no printing.
3178 !
3179 !***SEE ALSO XERRWD, XERRWV
3180 !***REFERENCES (NONE)
3181 !***ROUTINES CALLED IXSAV
3182 !***REVISION HISTORY (YYMMDD)
3183 ! 921118 DATE WRITTEN
3184 ! 930329 Added SLATEC format prologue. (FNF)
3185 ! 930407 Corrected SEE ALSO section. (FNF)
3186 ! 930922 Made user-callable, and other cosmetic changes. (FNF)
3187 !***END PROLOGUE XSETF
3188 !
3189 ! Subroutines called by XSETF.. None
3190 ! Function routine called by XSETF.. IXSAV
3191 !-----------------------------------------------------------------------
3192 !**End
3194 !
3195 !***FIRST EXECUTABLE STATEMENT XSETF
3197 RETURN
3198 !----------------------- End of Subroutine XSETF -----------------------
3199 END SUBROUTINE XSETF
3200 !DECK XSETUN
3202 !***BEGIN PROLOGUE XSETUN
3203 !***PURPOSE Reset the logical unit number for error messages.
3204 !***CATEGORY R3B
3205 !***TYPE ALL (XSETUN-A)
3206 !***KEYWORDS ERROR CONTROL
3207 !***DESCRIPTION
3208 !
3209 ! XSETUN sets the logical unit number for error messages to LUN.
3210 !
3211 !***AUTHOR Hindmarsh, Alan C., (LLNL)
3212 !***SEE ALSO XERRWD, XERRWV
3213 !***REFERENCES (NONE)
3214 !***ROUTINES CALLED IXSAV
3215 !***REVISION HISTORY (YYMMDD)
3216 ! 921118 DATE WRITTEN
3217 ! 930329 Added SLATEC format prologue. (FNF)
3218 ! 930407 Corrected SEE ALSO section. (FNF)
3219 ! 930922 Made user-callable, and other cosmetic changes. (FNF)
3220 !***END PROLOGUE XSETUN
3221 !
3222 ! Subroutines called by XSETUN.. None
3223 ! Function routine called by XSETUN.. IXSAV
3224 !-----------------------------------------------------------------------
3225 !**End
3227 !
3228 !***FIRST EXECUTABLE STATEMENT XSETUN
3230 RETURN
3231 !----------------------- End of Subroutine XSETUN ----------------------
3232 END SUBROUTINE XSETUN
3233 !DECK IXSAV
3235 !***BEGIN PROLOGUE IXSAV
3236 !***SUBSIDIARY
3237 !***PURPOSE Save and recall error message control parameters.
3238 !***CATEGORY R3C
3239 !***TYPE ALL (IXSAV-A)
3240 !***AUTHOR Hindmarsh, Alan C., (LLNL)
3241 !***DESCRIPTION
3242 !
3243 ! IXSAV saves and recalls one of two error message parameters:
3244 ! LUNIT, the logical unit number to which messages are printed, and
3245 ! MESFLG, the message print flag.
3246 ! This is a modification of the SLATEC library routine J4SAVE.
3247 !
3248 ! Saved local variables..
3249 ! LUNIT = Logical unit number for messages. The default is obtained
3250 ! by a call to IUMACH (may be machine-dependent).
3251 ! MESFLG = Print control flag..
3252 ! 1 means print all messages (the default).
3253 ! 0 means no printing.
3254 !
3255 ! On input..
3256 ! IPAR = Parameter indicator (1 for LUNIT, 2 for MESFLG).
3257 ! IVALUE = The value to be set for the parameter, if ISET = .TRUE.
3258 ! ISET = Logical flag to indicate whether to read or write.
3259 ! If ISET = .TRUE., the parameter will be given
3260 ! the value IVALUE. If ISET = .FALSE., the parameter
3261 ! will be unchanged, and IVALUE is a dummy argument.
3262 !
3263 ! On return..
3264 ! IXSAV = The (old) value of the parameter.
3265 !
3266 !***SEE ALSO XERRWD, XERRWV
3267 !***ROUTINES CALLED IUMACH
3268 !***REVISION HISTORY (YYMMDD)
3269 ! 921118 DATE WRITTEN
3270 ! 930329 Modified prologue to SLATEC format. (FNF)
3271 ! 930915 Added IUMACH call to get default output unit. (ACH)
3272 ! 930922 Minor cosmetic changes. (FNF)
3273 ! 010425 Type declaration for IUMACH added. (ACH)
3274 !***END PROLOGUE IXSAV
3275 !
3276 ! Subroutines called by IXSAV.. None
3277 ! Function routine called by IXSAV.. IUMACH
3278 !-----------------------------------------------------------------------
3279 !**End
3280 LOGICAL ISET
3282 !-----------------------------------------------------------------------
3284 !-----------------------------------------------------------------------
3285 ! The following Fortran-77 declaration is to cause the values of the
3286 ! listed (local) variables to be saved between calls to this routine.
3287 !-----------------------------------------------------------------------
3290 !
3291 !***FIRST EXECUTABLE STATEMENT IXSAV
3294 IXSAV = LUNIT
3296 ENDIF
3297 !
3299 IXSAV = MESFLG
3301 ENDIF
3302 !
3303 RETURN
3304 !----------------------- End of Function IXSAV -------------------------
3305 END FUNCTION IXSAV
3306 !DECK IUMACH
3308 !***BEGIN PROLOGUE IUMACH
3309 !***PURPOSE Provide standard output unit number.
3310 !***CATEGORY R1
3311 !***TYPE INTEGER (IUMACH-I)
3312 !***KEYWORDS MACHINE CONSTANTS
3313 !***AUTHOR Hindmarsh, Alan C., (LLNL)
3314 !***DESCRIPTION
3315 ! *Usage:
3316 ! INTEGER LOUT, IUMACH
3317 ! LOUT = IUMACH()
3318 !
3319 ! *Function Return Values:
3320 ! LOUT : the standard logical unit for Fortran output.
3321 !
3322 !***REFERENCES (NONE)
3323 !***ROUTINES CALLED (NONE)
3324 !***REVISION HISTORY (YYMMDD)
3325 ! 930915 DATE WRITTEN
3326 ! 930922 Made user-callable, and other cosmetic changes. (FNF)
3327 !***END PROLOGUE IUMACH
3328 !
3329 !*Internal Notes:
3330 ! The built-in value of 6 is standard on a wide range of Fortran
3331 ! systems. This may be machine-dependent.
3332 !**End
3333 !***FIRST EXECUTABLE STATEMENT IUMACH
3335 !
3336 RETURN
3337 !----------------------- End of Function IUMACH ------------------------
3338 END FUNCTION IUMACH
3340 !---- END OF SUBROUTINE DLSODE AND ITS INTERNAL PROCEDURES
3341 END SUBROUTINE DLSODE
3342 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3344 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3346 USE KPP_ROOT_Parameters
3347 USE KPP_ROOT_Global
3349 USE KPP_ROOT_Rates
3351 IMPLICIT NONE
3356 ! TOLD = TIME
3357 ! TIME = T
3358 ! CALL Update_SUN()
3359 ! CALL Update_RCONST()
3360 ! CALL Update_PHOTO()
3361 ! TIME = TOLD
3365 !Nfun=Nfun+1
3367 END SUBROUTINE FUN_CHEM
3370 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3372 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3374 USE KPP_ROOT_Parameters
3375 USE KPP_ROOT_Global
3376 USE KPP_ROOT_JacobianSP
3378 USE KPP_ROOT_Rates
3380 IMPLICIT NONE
3384 #ifdef FULL_ALGEBRA
3386 #else
3388 #endif
3390 ! TOLD = TIME
3391 ! TIME = T
3392 ! CALL Update_SUN()
3393 ! CALL Update_RCONST()
3394 ! CALL Update_PHOTO()
3395 ! TIME = TOLD
3397 #ifdef FULL_ALGEBRA
3402 END DO
3403 END DO
3406 END DO
3407 #else
3409 #endif
3410 !Njac=Njac+1
3412 END SUBROUTINE JAC_CHEM
3415 END MODULE KPP_ROOT_Integrator