Home

Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

kpp_lsode.f90 @ 4416

Last change on this file since 4416 was 2696, checked in by kanani, 7 years ago
Merge of branch palm4u into trunk
File size: 195.4 KB

Line
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
7
8	MODULE KPP_ROOT_Integrator
9
10	USE KPP_ROOT_Precision
11	USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME, ATOL, RTOL
12	USE KPP_ROOT_Parameters, ONLY: NVAR, NSPEC, NFIX, LU_NONZERO
13	USE KPP_ROOT_JacobianSP, ONLY: LU_DIAG
14	USE KPP_ROOT_LinearAlgebra, ONLY: KppDecomp, KppSolve, &
15	Set2zero, WLAMCH
16
17	IMPLICIT NONE
18	PUBLIC
19	SAVE
20
21	!~~~> Statistics on the work performed by the LSODE method
22	INTEGER :: Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng
23	INTEGER, PARAMETER :: ifun=1, ijac=2, istp=3, iacc=4, &
24	irej=5, idec=6, isol=7, isng=8, itexit=1, ihexit=2
25	! SDIRK method coefficients
26	KPP_REAL :: rkAlpha(5,4), rkBeta(5,4), rkD(4,5), &
27	rkGamma, rkA(5,5), rkB(5), rkC(5)
28
29	! mz_rs_20050717: TODO: use strings of IERR_NAMES for error messages
30	! description of the error numbers IERR
31	CHARACTER(LEN=50), PARAMETER, DIMENSION(-8:1) :: IERR_NAMES = (/ &
32	'Matrix is repeatedly singular ', & ! -8
33	'Step size too small ', & ! -7
34	'No of steps exceeds maximum bound ', & ! -6
35	'Improper tolerance values ', & ! -5
36	'FacMin/FacMax/FacRej must be positive ', & ! -4
37	'Hmin/Hmax/Hstart must be positive ', & ! -3
38	'Improper value for maximal no of Newton iterations', & ! -2
39	'Improper value for maximal no of steps ', & ! -1
40	' ', & ! 0 (not used)
41	'Success ' /) ! 1
42
43	CONTAINS
44
45	SUBROUTINE INTEGRATE( TIN, TOUT, &
46	ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U, IERR_U )
47
48	USE KPP_ROOT_Parameters
49	USE KPP_ROOT_Global
50	IMPLICIT NONE
51
52	KPP_REAL, INTENT(IN) :: TIN ! Start Time
53	KPP_REAL, INTENT(IN) :: TOUT ! End Time
54	! Optional input parameters and statistics
55	INTEGER, INTENT(IN), OPTIONAL :: ICNTRL_U(20)
56	KPP_REAL, INTENT(IN), OPTIONAL :: RCNTRL_U(20)
57	INTEGER, INTENT(OUT), OPTIONAL :: ISTATUS_U(20)
58	KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20)
59	INTEGER, INTENT(OUT), OPTIONAL :: IERR_U
60
61	KPP_REAL :: RCNTRL(20), RSTATUS(20)
62	INTEGER :: ICNTRL(20), ISTATUS(20), IERR
63	!!$ INTEGER, SAVE :: Ntotal = 0
64
65	ICNTRL(:) = 0
66	RCNTRL(:) = 0.0_dp
67	ISTATUS(:) = 0
68	RSTATUS(:) = 0.0_dp
69
70	ICNTRL(5) = 2 ! maximal order
71
72	! If optional parameters are given, and if they are >0,
73	! then they overwrite default settings.
74	IF (PRESENT(ICNTRL_U)) THEN
75	WHERE(ICNTRL_U(:) > 0) ICNTRL(:) = ICNTRL_U(:)
76	END IF
77	IF (PRESENT(RCNTRL_U)) THEN
78	WHERE(RCNTRL_U(:) > 0) RCNTRL(:) = RCNTRL_U(:)
79	END IF
80
81	CALL KppLsode( TIN,TOUT,VAR,RTOL,ATOL, &
82	RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR )
83
84	! INTEGER ITASK, ISTATE, NG, NEQ, IOUT, JROOT, ISTATS, &
85	! IERROR, I
86	! DOUBLE PRECISION ATOL, RTOL, RSTATS, T, TOUT, Y
87	! DIMENSION JROOT(2)
88	! TYPE(VODE_OPTS) :: OPTIONS
89	! IERROR = 0
90	! ITASK = 1
91	! ISTATE = 1
92	! NG = 2
93	! OPTIONS = SET_OPTS(DENSE_J=.TRUE.,RELERR=RTOL, &
94	! ABSERR_VECTOR=ATOL,NEVENTS=NG)
95	! CALL DVODE_F90(FEX,NEQ,Y,TIN,TOUT,ITASK,ISTATE,OPTIONS,G_FCN=GEX)
96	! CALL GET_STATS(RSTATS, ISTATS, NG, JROOT)
97
98
99	STEPMIN = RSTATUS(ihexit) ! Save last step
100
101	! if optional parameters are given for output they to return information
102	IF (PRESENT(ISTATUS_U)) ISTATUS_U(:) = ISTATUS(1:20)
103	IF (PRESENT(RSTATUS_U)) RSTATUS_U(:) = RSTATUS(1:20)
104	IF (PRESENT(IERR_U)) THEN
105	IF (IERR==2) THEN ! DLSODE returns "2" after successful completion
106	IERR_U = 1 ! IERR_U will return "1" for successful completion
107	ELSE
108	IERR_U = IERR
109	ENDIF
110	ENDIF
111
112	END SUBROUTINE INTEGRATE
113
114	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
115	SUBROUTINE KppLsode( TIN,TOUT,Y,RelTol,AbsTol, &
116	RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR )
117	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
118	!
119	!~~~> INPUT PARAMETERS:
120	!
121	! Note: For input parameters equal to zero the default values of the
122	! corresponding variables are used.
123	!
124	! Note: For input parameters equal to zero the default values of the
125	! corresponding variables are used.
126	!~~~>
127	! ICNTRL(1) = not used
128	!
129	! ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors
130	! = 1: AbsTol, RelTol are scalars
131	!
132	! ICNTRL(3) = not used
133	!
134	! ICNTRL(4) -> maximum number of integration steps
135	! For ICNTRL(4)=0 the default value of 100000 is used
136	!
137	! ICNTRL(5) -> maximum order of the integration formula allowed
138	!
139	!~~~> Real parameters
140	!
141	! RCNTRL(1) -> Hmin, lower bound for the integration step size
142	! It is strongly recommended to keep Hmin = ZERO
143	! RCNTRL(2) -> Hmax, upper bound for the integration step size
144	! RCNTRL(3) -> Hstart, starting value for the integration step size
145	!
146	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
147	!
148	!~~~> OUTPUT PARAMETERS:
149	!
150	! Note: each call to Rosenbrock adds the current no. of fcn calls
151	! to previous value of ISTATUS(1), and similar for the other params.
152	! Set ISTATUS(1:10) = 0 before call to avoid this accumulation.
153	!
154	! ISTATUS(1) = No. of function calls
155	! ISTATUS(2) = No. of jacobian calls
156	! ISTATUS(3) = No. of steps
157	!
158	! RSTATUS(1) -> Texit, the time corresponding to the
159	! computed Y upon return
160	! RSTATUS(2) -> Hexit, last predicted step before exit
161	! For multiple restarts, use Hexit as Hstart in the following run
162	!
163	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
164
165	IMPLICIT NONE
166	KPP_REAL :: Y(NVAR), AbsTol(NVAR), RelTol(NVAR), TIN, TOUT
167	KPP_REAL :: RCNTRL(20), RSTATUS(20)
168	INTEGER :: ICNTRL(20), ISTATUS(20)
169	INTEGER, PARAMETER :: LRW = 25 + 9NVAR+2NVAR*NVAR, &
170	LIW = 32 + NVAR
171	KPP_REAL :: RWORK(LRW), RPAR(1)
172	INTEGER :: IWORK(LIW), IPAR(1), ITOL, ITASK, &
173	IERR, IOPT, MF
174
175	!~~~> NORMAL COMPUTATION
176	ITASK = 1
177	IERR = 1
178	IOPT = 1 ! 0=no/1=use optional input
179
180	RWORK(1:30) = 0.0d0
181	IWORK(1:30) = 0
182
183	IF (ICNTRL(2)==0) THEN
184	ITOL = 4 ! Abs/RelTol are both vectors
185	ELSE
186	ITOL = 1 ! Abs/RelTol are both scalars
187	END IF
188	IWORK(6) = ICNTRL(4) ! max number of internal steps
189	IWORK(5) = ICNTRL(5) ! maximal order
190
191	MF = 21 !~~~> stiff case, analytic full Jacobian
192
193	RWORK(5) = RCNTRL(3) ! Hstart
194	RWORK(6) = RCNTRL(2) ! Hmax
195	RWORK(7) = RCNTRL(1) ! Hmin
196
197	CALL DLSODE ( FUN_CHEM, NVAR, Y, TIN, TOUT, ITOL, RelTol, AbsTol, ITASK,&
198	IERR, IOPT, RWORK, LRW, IWORK, LIW, JAC_CHEM, MF)
199
200	ISTATUS(1) = IWORK(12) ! Number of function evaluations
201	ISTATUS(2) = IWORK(13) ! Number of Jacobian evaluations
202	ISTATUS(3) = IWORK(11) ! Number of steps
203
204	RSTATUS(1) = TOUT ! mz_rs_20050717
205	RSTATUS(2) = RWORK(11) ! mz_rs_20050717
206
207	END SUBROUTINE KppLsode
208	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
209
210
211	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
212	!DECK DLSODE
213	SUBROUTINE DLSODE (F, NEQ, Y, T, TOUT, ITOL, RelTol, AbsTol, ITASK, &
214	ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF)
215	EXTERNAL F, JAC
216	INTEGER NEQ, ITOL, ITASK, ISTATE, IOPT, LRW, LIW, IWORK(LIW), MF
217	KPP_REAL Y(), T, TOUT, RelTol(), AbsTol(*), RWORK(LRW)
218	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
219	!***BEGIN PROLOGUE DLSODE
220	!***PURPOSE Livermore Solver for Ordinary Differential Equations.
221	! DLSODE solves the initial-value problem for stiff or
222	! nonstiff systems of first-order ODE's,
223	! dy/dt = f(t,y), or, in component form,
224	! dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(N)), i=1,...,N.
225	!***CATEGORY I1A
226	!***TYPE KPP_REAL (SLSODE-S, DLSODE-D)
227	!***KEYWORDS ORDINARY DIFFERENTIAL EQUATIONS, INITIAL VALUE PROBLEM,
228	! STIFF, NONSTIFF
229	!***AUTHOR Hindmarsh, Alan C., (LLNL)
230	! Center for Applied Scientific Computing, L-561
231	! Lawrence Livermore National Laboratory
232	! Livermore, CA 94551.
233	!***DESCRIPTION
234	!
235	! NOTE: The "Usage" and "Arguments" sections treat only a subset of
236	! available options, in condensed fashion. The options
237	! covered and the information supplied will support most
238	! standard uses of DLSODE.
239	!
240	! For more sophisticated uses, full details on all options are
241	! given in the concluding section, headed "Long Description."
242	! A synopsis of the DLSODE Long Description is provided at the
243	! beginning of that section; general topics covered are:
244	! - Elements of the call sequence; optional input and output
245	! - Optional supplemental routines in the DLSODE package
246	! - internal COMMON block
247	!
248	! *Usage:
249	! Communication between the user and the DLSODE package, for normal
250	! situations, is summarized here. This summary describes a subset
251	! of the available options. See "Long Description" for complete
252	! details, including optional communication, nonstandard options,
253	! and instructions for special situations.
254	!
255	! A sample program is given in the "Examples" section.
256	!
257	! Refer to the argument descriptions for the definitions of the
258	! quantities that appear in the following sample declarations.
259	!
260	! For MF = 10,
261	! PARAMETER (LRW = 20 + 16*NEQ, LIW = 20)
262	! For MF = 21 or 22,
263	! PARAMETER (LRW = 22 + 9NEQ + NEQ*2, LIW = 20 + NEQ)
264	! For MF = 24 or 25,
265	! PARAMETER (LRW = 22 + 10NEQ + (2ML+MU)*NEQ,
266	! * LIW = 20 + NEQ)
267	!
268	! EXTERNAL F, JAC
269	! INTEGER NEQ, ITOL, ITASK, ISTATE, IOPT, LRW, IWORK(LIW),
270	! * LIW, MF
271	! KPP_REAL Y(NEQ), T, TOUT, RelTol, AbsTol(ntol), RWORK(LRW)
272	!
273	! CALL DLSODE (F, NEQ, Y, T, TOUT, ITOL, RelTol, AbsTol, ITASK,
274	! * ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF)
275	!
276	! *Arguments:
277	! F :EXT Name of subroutine for right-hand-side vector f.
278	! This name must be declared EXTERNAL in calling
279	! program. The form of F must be:
280	!
281	! SUBROUTINE F (NEQ, T, Y, YDOT)
282	! INTEGER NEQ
283	! KPP_REAL T, Y(), YDOT()
284	!
285	! The inputs are NEQ, T, Y. F is to set
286	!
287	! YDOT(i) = f(i,T,Y(1),Y(2),...,Y(NEQ)),
288	! i = 1, ..., NEQ .
289	!
290	! NEQ :IN Number of first-order ODE's.
291	!
292	! Y :INOUT Array of values of the y(t) vector, of length NEQ.
293	! Input: For the first call, Y should contain the
294	! values of y(t) at t = T. (Y is an input
295	! variable only if ISTATE = 1.)
296	! Output: On return, Y will contain the values at the
297	! new t-value.
298	!
299	! T :INOUT Value of the independent variable. On return it
300	! will be the current value of t (normally TOUT).
301	!
302	! TOUT :IN Next point where output is desired (.NE. T).
303	!
304	! ITOL :IN 1 or 2 according as AbsTol (below) is a scalar or
305	! an array.
306	!
307	! RelTol :IN Relative tolerance parameter (scalar).
308	!
309	! AbsTol :IN Absolute tolerance parameter (scalar or array).
310	! If ITOL = 1, AbsTol need not be dimensioned.
311	! If ITOL = 2, AbsTol must be dimensioned at least NEQ.
312	!
313	! The estimated local error in Y(i) will be controlled
314	! so as to be roughly less (in magnitude) than
315	!
316	! EWT(i) = RelTol*ABS(Y(i)) + AbsTol if ITOL = 1, or
317	! EWT(i) = RelTol*ABS(Y(i)) + AbsTol(i) if ITOL = 2.
318	!
319	! Thus the local error test passes if, in each
320	! component, either the absolute error is less than
321	! AbsTol (or AbsTol(i)), or the relative error is less
322	! than RelTol.
323	!
324	! Use RelTol = 0.0 for pure absolute error control, and
325	! use AbsTol = 0.0 (or AbsTol(i) = 0.0) for pure relative
326	! error control. Caution: Actual (global) errors may
327	! exceed these local tolerances, so choose them
328	! conservatively.
329	!
330	! ITASK :IN Flag indicating the task DLSODE is to perform.
331	! Use ITASK = 1 for normal computation of output
332	! values of y at t = TOUT.
333	!
334	! ISTATE:INOUT Index used for input and output to specify the state
335	! of the calculation.
336	! Input:
337	! 1 This is the first call for a problem.
338	! 2 This is a subsequent call.
339	! Output:
340	! 1 Nothing was done, because TOUT was equal to T.
341	! 2 DLSODE was successful (otherwise, negative).
342	! Note that ISTATE need not be modified after a
343	! successful return.
344	! -1 Excess work done on this call (perhaps wrong
345	! MF).
346	! -2 Excess accuracy requested (tolerances too
347	! small).
348	! -3 Illegal input detected (see printed message).
349	! -4 Repeated error test failures (check all
350	! inputs).
351	! -5 Repeated convergence failures (perhaps bad
352	! Jacobian supplied or wrong choice of MF or
353	! tolerances).
354	! -6 Error weight became zero during problem
355	! (solution component i vanished, and AbsTol or
356	! AbsTol(i) = 0.).
357	!
358	! IOPT :IN Flag indicating whether optional inputs are used:
359	! 0 No.
360	! 1 Yes. (See "Optional inputs" under "Long
361	! Description," Part 1.)
362	!
363	! RWORK :WORK Real work array of length at least:
364	! 20 + 16*NEQ for MF = 10,
365	! 22 + 9NEQ + NEQ*2 for MF = 21 or 22,
366	! 22 + 10NEQ + (2ML + MU)*NEQ for MF = 24 or 25.
367	!
368	! LRW :IN Declared length of RWORK (in user's DIMENSION
369	! statement).
370	!
371	! IWORK :WORK Integer work array of length at least:
372	! 20 for MF = 10,
373	! 20 + NEQ for MF = 21, 22, 24, or 25.
374	!
375	! If MF = 24 or 25, input in IWORK(1),IWORK(2) the
376	! lower and upper Jacobian half-bandwidths ML,MU.
377	!
378	! On return, IWORK contains information that may be
379	! of interest to the user:
380	!
381	! Name Location Meaning
382	! ----- --------- -----------------------------------------
383	! NST IWORK(11) Number of steps taken for the problem so
384	! far.
385	! NFE IWORK(12) Number of f evaluations for the problem
386	! so far.
387	! NJE IWORK(13) Number of Jacobian evaluations (and of
388	! matrix LU decompositions) for the problem
389	! so far.
390	! NQU IWORK(14) Method order last used (successfully).
391	! LENRW IWORK(17) Length of RWORK actually required. This
392	! is defined on normal returns and on an
393	! illegal input return for insufficient
394	! storage.
395	! LENIW IWORK(18) Length of IWORK actually required. This
396	! is defined on normal returns and on an
397	! illegal input return for insufficient
398	! storage.
399	!
400	! LIW :IN Declared length of IWORK (in user's DIMENSION
401	! statement).
402	!
403	! JAC :EXT Name of subroutine for Jacobian matrix (MF =
404	! 21 or 24). If used, this name must be declared
405	! EXTERNAL in calling program. If not used, pass a
406	! dummy name. The form of JAC must be:
407	!
408	! SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD)
409	! INTEGER NEQ, ML, MU, NROWPD
410	! KPP_REAL T, Y(), PD(NROWPD,)
411	!
412	! See item c, under "Description" below for more
413	! information about JAC.
414	!
415	! MF :IN Method flag. Standard values are:
416	! 10 Nonstiff (Adams) method, no Jacobian used.
417	! 21 Stiff (BDF) method, user-supplied full Jacobian.
418	! 22 Stiff method, internally generated full
419	! Jacobian.
420	! 24 Stiff method, user-supplied banded Jacobian.
421	! 25 Stiff method, internally generated banded
422	! Jacobian.
423	!
424	! *Description:
425	! DLSODE solves the initial value problem for stiff or nonstiff
426	! systems of first-order ODE's,
427	!
428	! dy/dt = f(t,y) ,
429	!
430	! or, in component form,
431	!
432	! dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(NEQ))
433	! (i = 1, ..., NEQ) .
434	!
435	! DLSODE is a package based on the GEAR and GEARB packages, and on
436	! the October 23, 1978, version of the tentative ODEPACK user
437	! interface standard, with minor modifications.
438	!
439	! The steps in solving such a problem are as follows.
440	!
441	! a. First write a subroutine of the form
442	!
443	! SUBROUTINE F (NEQ, T, Y, YDOT)
444	! INTEGER NEQ
445	! KPP_REAL T, Y(), YDOT()
446	!
447	! which supplies the vector function f by loading YDOT(i) with
448	! f(i).
449	!
450	! b. Next determine (or guess) whether or not the problem is stiff.
451	! Stiffness occurs when the Jacobian matrix df/dy has an
452	! eigenvalue whose real part is negative and large in magnitude
453	! compared to the reciprocal of the t span of interest. If the
454	! problem is nonstiff, use method flag MF = 10. If it is stiff,
455	! there are four standard choices for MF, and DLSODE requires the
456	! Jacobian matrix in some form. This matrix is regarded either
457	! as full (MF = 21 or 22), or banded (MF = 24 or 25). In the
458	! banded case, DLSODE requires two half-bandwidth parameters ML
459	! and MU. These are, respectively, the widths of the lower and
460	! upper parts of the band, excluding the main diagonal. Thus the
461	! band consists of the locations (i,j) with
462	!
463	! i - ML <= j <= i + MU ,
464	!
465	! and the full bandwidth is ML + MU + 1 .
466	!
467	! c. If the problem is stiff, you are encouraged to supply the
468	! Jacobian directly (MF = 21 or 24), but if this is not feasible,
469	! DLSODE will compute it internally by difference quotients (MF =
470	! 22 or 25). If you are supplying the Jacobian, write a
471	! subroutine of the form
472	!
473	! SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD)
474	! INTEGER NEQ, ML, MU, NRWOPD
475	! KPP_REAL T, Y(), PD(NROWPD,)
476	!
477	! which provides df/dy by loading PD as follows:
478	! - For a full Jacobian (MF = 21), load PD(i,j) with df(i)/dy(j),
479	! the partial derivative of f(i) with respect to y(j). (Ignore
480	! the ML and MU arguments in this case.)
481	! - For a banded Jacobian (MF = 24), load PD(i-j+MU+1,j) with
482	! df(i)/dy(j); i.e., load the diagonal lines of df/dy into the
483	! rows of PD from the top down.
484	! - In either case, only nonzero elements need be loaded.
485	!
486	! d. Write a main program that calls subroutine DLSODE once for each
487	! point at which answers are desired. This should also provide
488	! for possible use of logical unit 6 for output of error messages
489	! by DLSODE.
490	!
491	! Before the first call to DLSODE, set ISTATE = 1, set Y and T to
492	! the initial values, and set TOUT to the first output point. To
493	! continue the integration after a successful return, simply
494	! reset TOUT and call DLSODE again. No other parameters need be
495	! reset.
496	!
497	! *Examples:
498	! The following is a simple example problem, with the coding needed
499	! for its solution by DLSODE. The problem is from chemical kinetics,
500	! and consists of the following three rate equations:
501	!
502	! dy1/dt = -.04y1 + 1.E4y2*y3
503	! dy2/dt = .04y1 - 1.E4y2y3 - 3.E7y2**2
504	! dy3/dt = 3.E7y2*2
505	!
506	! on the interval from t = 0.0 to t = 4.E10, with initial conditions
507	! y1 = 1.0, y2 = y3 = 0. The problem is stiff.
508	!
509	! The following coding solves this problem with DLSODE, using
510	! MF = 21 and printing results at t = .4, 4., ..., 4.E10. It uses
511	! ITOL = 2 and AbsTol much smaller for y2 than for y1 or y3 because y2
512	! has much smaller values. At the end of the run, statistical
513	! quantities of interest are printed.
514	!
515	! EXTERNAL FEX, JEX
516	! INTEGER IOPT, IOUT, ISTATE, ITASK, ITOL, IWORK(23), LIW, LRW,
517	! * MF, NEQ
518	! KPP_REAL AbsTol(3), RelTol, RWORK(58), T, TOUT, Y(3)
519	! NEQ = 3
520	! Y(1) = 1.D0
521	! Y(2) = 0.D0
522	! Y(3) = 0.D0
523	! T = 0.D0
524	! TOUT = .4D0
525	! ITOL = 2
526	! RelTol = 1.D-4
527	! AbsTol(1) = 1.D-6
528	! AbsTol(2) = 1.D-10
529	! AbsTol(3) = 1.D-6
530	! ITASK = 1
531	! ISTATE = 1
532	! IOPT = 0
533	! LRW = 58
534	! LIW = 23
535	! MF = 21
536	! DO 40 IOUT = 1,12
537	! CALL DLSODE (FEX, NEQ, Y, T, TOUT, ITOL, RelTol, AbsTol, ITASK,
538	! * ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JEX, MF)
539	! WRITE(6,20) T, Y(1), Y(2), Y(3)
540	! 20 FORMAT(' At t =',D12.4,' y =',3D14.6)
541	! IF (ISTATE .LT. 0) GO TO 80
542	! 40 TOUT = TOUT*10.D0
543	! WRITE(6,60) IWORK(11), IWORK(12), IWORK(13)
544	! 60 FORMAT(/' No. steps =',i4,', No. f-s =',i4,', No. J-s =',i4)
545	! STOP
546	! 80 WRITE(6,90) ISTATE
547	! 90 FORMAT(///' Error halt.. ISTATE =',I3)
548	! STOP
549	! END
550	!
551	! SUBROUTINE FEX (NEQ, T, Y, YDOT)
552	! INTEGER NEQ
553	! KPP_REAL T, Y(3), YDOT(3)
554	! YDOT(1) = -.04D0Y(1) + 1.D4Y(2)*Y(3)
555	! YDOT(3) = 3.D7Y(2)Y(2)
556	! YDOT(2) = -YDOT(1) - YDOT(3)
557	! RETURN
558	! END
559	!
560	! SUBROUTINE JEX (NEQ, T, Y, ML, MU, PD, NRPD)
561	! INTEGER NEQ, ML, MU, NRPD
562	! KPP_REAL T, Y(3), PD(NRPD,3)
563	! PD(1,1) = -.04D0
564	! PD(1,2) = 1.D4*Y(3)
565	! PD(1,3) = 1.D4*Y(2)
566	! PD(2,1) = .04D0
567	! PD(2,3) = -PD(1,3)
568	! PD(3,2) = 6.D7*Y(2)
569	! PD(2,2) = -PD(1,2) - PD(3,2)
570	! RETURN
571	! END
572	!
573	! The output from this program (on a Cray-1 in single precision)
574	! is as follows.
575	!
576	! At t = 4.0000e-01 y = 9.851726e-01 3.386406e-05 1.479357e-02
577	! At t = 4.0000e+00 y = 9.055142e-01 2.240418e-05 9.446344e-02
578	! At t = 4.0000e+01 y = 7.158050e-01 9.184616e-06 2.841858e-01
579	! At t = 4.0000e+02 y = 4.504846e-01 3.222434e-06 5.495122e-01
580	! At t = 4.0000e+03 y = 1.831701e-01 8.940379e-07 8.168290e-01
581	! At t = 4.0000e+04 y = 3.897016e-02 1.621193e-07 9.610297e-01
582	! At t = 4.0000e+05 y = 4.935213e-03 1.983756e-08 9.950648e-01
583	! At t = 4.0000e+06 y = 5.159269e-04 2.064759e-09 9.994841e-01
584	! At t = 4.0000e+07 y = 5.306413e-05 2.122677e-10 9.999469e-01
585	! At t = 4.0000e+08 y = 5.494530e-06 2.197825e-11 9.999945e-01
586	! At t = 4.0000e+09 y = 5.129458e-07 2.051784e-12 9.999995e-01
587	! At t = 4.0000e+10 y = -7.170603e-08 -2.868241e-13 1.000000e+00
588	!
589	! No. steps = 330, No. f-s = 405, No. J-s = 69
590	!
591	! *Accuracy:
592	! The accuracy of the solution depends on the choice of tolerances
593	! RelTol and AbsTol. Actual (global) errors may exceed these local
594	! tolerances, so choose them conservatively.
595	!
596	! *Cautions:
597	! The work arrays should not be altered between calls to DLSODE for
598	! the same problem, except possibly for the conditional and optional
599	! inputs.
600	!
601	! *Portability:
602	! Since NEQ is dimensioned inside DLSODE, some compilers may object
603	! to a call to DLSODE with NEQ a scalar variable. In this event,
604	! use DIMENSION NEQ. Similar remarks apply to RelTol and AbsTol.
605	!
606	! Note to Cray users:
607	! For maximum efficiency, use the CFT77 compiler. Appropriate
608	! compiler optimization directives have been inserted for CFT77.
609	!
610	! *Reference:
611	! Alan C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE
612	! Solvers," in Scientific Computing, R. S. Stepleman, et al., Eds.
613	! (North-Holland, Amsterdam, 1983), pp. 55-64.
614	!
615	! *Long Description:
616	! The following complete description of the user interface to
617	! DLSODE consists of four parts:
618	!
619	! 1. The call sequence to subroutine DLSODE, which is a driver
620	! routine for the solver. This includes descriptions of both
621	! the call sequence arguments and user-supplied routines.
622	! Following these descriptions is a description of optional
623	! inputs available through the call sequence, and then a
624	! description of optional outputs in the work arrays.
625	!
626	! 2. Descriptions of other routines in the DLSODE package that may
627	! be (optionally) called by the user. These provide the ability
628	! to alter error message handling, save and restore the internal
629	! COMMON, and obtain specified derivatives of the solution y(t).
630	!
631	! 3. Descriptions of COMMON block to be declared in overlay or
632	! similar environments, or to be saved when doing an interrupt
633	! of the problem and continued solution later.
634	!
635	! 4. Description of two routines in the DLSODE package, either of
636	! which the user may replace with his own version, if desired.
637	! These relate to the measurement of errors.
638	!
639	!
640	! Part 1. Call Sequence
641	! ----------------------
642	!
643	! Arguments
644	! ---------
645	! The call sequence parameters used for input only are
646	!
647	! F, NEQ, TOUT, ITOL, RelTol, AbsTol, ITASK, IOPT, LRW, LIW, JAC, MF,
648	!
649	! and those used for both input and output are
650	!
651	! Y, T, ISTATE.
652	!
653	! The work arrays RWORK and IWORK are also used for conditional and
654	! optional inputs and optional outputs. (The term output here
655	! refers to the return from subroutine DLSODE to the user's calling
656	! program.)
657	!
658	! The legality of input parameters will be thoroughly checked on the
659	! initial call for the problem, but not checked thereafter unless a
660	! change in input parameters is flagged by ISTATE = 3 on input.
661	!
662	! The descriptions of the call arguments are as follows.
663	!
664	! F The name of the user-supplied subroutine defining the ODE
665	! system. The system must be put in the first-order form
666	! dy/dt = f(t,y), where f is a vector-valued function of
667	! the scalar t and the vector y. Subroutine F is to compute
668	! the function f. It is to have the form
669	!
670	! SUBROUTINE F (NEQ, T, Y, YDOT)
671	! KPP_REAL T, Y(), YDOT()
672	!
673	! where NEQ, T, and Y are input, and the array YDOT =
674	! f(T,Y) is output. Y and YDOT are arrays of length NEQ.
675	! Subroutine F should not alter Y(1),...,Y(NEQ). F must be
676	! declared EXTERNAL in the calling program.
677	!
678	! Subroutine F may access user-defined quantities in
679	! NEQ(2),... and/or in Y(NEQ+1),..., if NEQ is an array
680	! (dimensioned in F) and/or Y has length exceeding NEQ.
681	! See the descriptions of NEQ and Y below.
682	!
683	! If quantities computed in the F routine are needed
684	! externally to DLSODE, an extra call to F should be made
685	! for this purpose, for consistent and accurate results.
686	! If only the derivative dy/dt is needed, use DINTDY
687	! instead.
688	!
689	! NEQ The size of the ODE system (number of first-order
690	! ordinary differential equations). Used only for input.
691	! NEQ may be decreased, but not increased, during the
692	! problem. If NEQ is decreased (with ISTATE = 3 on input),
693	! the remaining components of Y should be left undisturbed,
694	! if these are to be accessed in F and/or JAC.
695	!
696	! Normally, NEQ is a scalar, and it is generally referred
697	! to as a scalar in this user interface description.
698	! However, NEQ may be an array, with NEQ set to the
699	! system size. (The DLSODE package accesses only NEQ.)
700	! In either case, this parameter is passed as the NEQ
701	! argument in all calls to F and JAC. Hence, if it is an
702	! array, locations NEQ(2),... may be used to store other
703	! integer data and pass it to F and/or JAC. Subroutines
704	! F and/or JAC must include NEQ in a DIMENSION statement
705	! in that case.
706	!
707	! Y A real array for the vector of dependent variables, of
708	! length NEQ or more. Used for both input and output on
709	! the first call (ISTATE = 1), and only for output on
710	! other calls. On the first call, Y must contain the
711	! vector of initial values. On output, Y contains the
712	! computed solution vector, evaluated at T. If desired,
713	! the Y array may be used for other purposes between
714	! calls to the solver.
715	!
716	! This array is passed as the Y argument in all calls to F
717	! and JAC. Hence its length may exceed NEQ, and locations
718	! Y(NEQ+1),... may be used to store other real data and
719	! pass it to F and/or JAC. (The DLSODE package accesses
720	! only Y(1),...,Y(NEQ).)
721	!
722	! T The independent variable. On input, T is used only on
723	! the first call, as the initial point of the integration.
724	! On output, after each call, T is the value at which a
725	! computed solution Y is evaluated (usually the same as
726	! TOUT). On an error return, T is the farthest point
727	! reached.
728	!
729	! TOUT The next value of T at which a computed solution is
730	! desired. Used only for input.
731	!
732	! When starting the problem (ISTATE = 1), TOUT may be equal
733	! to T for one call, then should not equal T for the next
734	! call. For the initial T, an input value of TOUT .NE. T
735	! is used in order to determine the direction of the
736	! integration (i.e., the algebraic sign of the step sizes)
737	! and the rough scale of the problem. Integration in
738	! either direction (forward or backward in T) is permitted.
739	!
740	! If ITASK = 2 or 5 (one-step modes), TOUT is ignored
741	! after the first call (i.e., the first call with
742	! TOUT .NE. T). Otherwise, TOUT is required on every call.
743	!
744	! If ITASK = 1, 3, or 4, the values of TOUT need not be
745	! monotone, but a value of TOUT which backs up is limited
746	! to the current internal T interval, whose endpoints are
747	! TCUR - HU and TCUR. (See "Optional Outputs" below for
748	! TCUR and HU.)
749	!
750	!
751	! ITOL An indicator for the type of error control. See
752	! description below under AbsTol. Used only for input.
753	!
754	! RelTol A relative error tolerance parameter, either a scalar or
755	! an array of length NEQ. See description below under
756	! AbsTol. Input only.
757	!
758	! AbsTol An absolute error tolerance parameter, either a scalar or
759	! an array of length NEQ. Input only.
760	!
761	! The input parameters ITOL, RelTol, and AbsTol determine the
762	! error control performed by the solver. The solver will
763	! control the vector e = (e(i)) of estimated local errors
764	! in Y, according to an inequality of the form
765	!
766	! rms-norm of ( e(i)/EWT(i) ) <= 1,
767	!
768	! where
769	!
770	! EWT(i) = RelTol(i)*ABS(Y(i)) + AbsTol(i),
771	!
772	! and the rms-norm (root-mean-square norm) here is
773	!
774	! rms-norm(v) = SQRT(sum v(i)**2 / NEQ).
775	!
776	! Here EWT = (EWT(i)) is a vector of weights which must
777	! always be positive, and the values of RelTol and AbsTol
778	! should all be nonnegative. The following table gives the
779	! types (scalar/array) of RelTol and AbsTol, and the
780	! corresponding form of EWT(i).
781	!
782	! ITOL RelTol AbsTol EWT(i)
783	! ---- ------ ------ -----------------------------
784	! 1 scalar scalar RelTol*ABS(Y(i)) + AbsTol
785	! 2 scalar array RelTol*ABS(Y(i)) + AbsTol(i)
786	! 3 array scalar RelTol(i)*ABS(Y(i)) + AbsTol
787	! 4 array array RelTol(i)*ABS(Y(i)) + AbsTol(i)
788	!
789	! When either of these parameters is a scalar, it need not
790	! be dimensioned in the user's calling program.
791	!
792	! If none of the above choices (with ITOL, RelTol, and AbsTol
793	! fixed throughout the problem) is suitable, more general
794	! error controls can be obtained by substituting
795	! user-supplied routines for the setting of EWT and/or for
796	! the norm calculation. See Part 4 below.
797	!
798	! If global errors are to be estimated by making a repeated
799	! run on the same problem with smaller tolerances, then all
800	! components of RelTol and AbsTol (i.e., of EWT) should be
801	! scaled down uniformly.
802	!
803	! ITASK An index specifying the task to be performed. Input
804	! only. ITASK has the following values and meanings:
805	! 1 Normal computation of output values of y(t) at
806	! t = TOUT (by overshooting and interpolating).
807	! 2 Take one step only and return.
808	! 3 Stop at the first internal mesh point at or beyond
809	! t = TOUT and return.
810	! 4 Normal computation of output values of y(t) at
811	! t = TOUT but without overshooting t = TCRIT. TCRIT
812	! must be input as RWORK(1). TCRIT may be equal to or
813	! beyond TOUT, but not behind it in the direction of
814	! integration. This option is useful if the problem
815	! has a singularity at or beyond t = TCRIT.
816	! 5 Take one step, without passing TCRIT, and return.
817	! TCRIT must be input as RWORK(1).
818	!
819	! Note: If ITASK = 4 or 5 and the solver reaches TCRIT
820	! (within roundoff), it will return T = TCRIT (exactly) to
821	! indicate this (unless ITASK = 4 and TOUT comes before
822	! TCRIT, in which case answers at T = TOUT are returned
823	! first).
824	!
825	! ISTATE An index used for input and output to specify the state
826	! of the calculation.
827	!
828	! On input, the values of ISTATE are as follows:
829	! 1 This is the first call for the problem
830	! (initializations will be done). See "Note" below.
831	! 2 This is not the first call, and the calculation is to
832	! continue normally, with no change in any input
833	! parameters except possibly TOUT and ITASK. (If ITOL,
834	! RelTol, and/or AbsTol are changed between calls with
835	! ISTATE = 2, the new values will be used but not
836	! tested for legality.)
837	! 3 This is not the first call, and the calculation is to
838	! continue normally, but with a change in input
839	! parameters other than TOUT and ITASK. Changes are
840	! allowed in NEQ, ITOL, RelTol, AbsTol, IOPT, LRW, LIW, MF,
841	! ML, MU, and any of the optional inputs except H0.
842	! (See IWORK description for ML and MU.)
843	!
844	! Note: A preliminary call with TOUT = T is not counted as
845	! a first call here, as no initialization or checking of
846	! input is done. (Such a call is sometimes useful for the
847	! purpose of outputting the initial conditions.) Thus the
848	! first call for which TOUT .NE. T requires ISTATE = 1 on
849	! input.
850	!
851	! On output, ISTATE has the following values and meanings:
852	! 1 Nothing was done, as TOUT was equal to T with
853	! ISTATE = 1 on input.
854	! 2 The integration was performed successfully.
855	! -1 An excessive amount of work (more than MXSTEP steps)
856	! was done on this call, before completing the
857	! requested task, but the integration was otherwise
858	! successful as far as T. (MXSTEP is an optional input
859	! and is normally 500.) To continue, the user may
860	! simply reset ISTATE to a value >1 and call again (the
861	! excess work step counter will be reset to 0). In
862	! addition, the user may increase MXSTEP to avoid this
863	! error return; see "Optional Inputs" below.
864	! -2 Too much accuracy was requested for the precision of
865	! the machine being used. This was detected before
866	! completing the requested task, but the integration
867	! was successful as far as T. To continue, the
868	! tolerance parameters must be reset, and ISTATE must
869	! be set to 3. The optional output TOLSF may be used
870	! for this purpose. (Note: If this condition is
871	! detected before taking any steps, then an illegal
872	! input return (ISTATE = -3) occurs instead.)
873	! -3 Illegal input was detected, before taking any
874	! integration steps. See written message for details.
875	! (Note: If the solver detects an infinite loop of
876	! calls to the solver with illegal input, it will cause
877	! the run to stop.)
878	! -4 There were repeated error-test failures on one
879	! attempted step, before completing the requested task,
880	! but the integration was successful as far as T. The
881	! problem may have a singularity, or the input may be
882	! inappropriate.
883	! -5 There were repeated convergence-test failures on one
884	! attempted step, before completing the requested task,
885	! but the integration was successful as far as T. This
886	! may be caused by an inaccurate Jacobian matrix, if
887	! one is being used.
888	! -6 EWT(i) became zero for some i during the integration.
889	! Pure relative error control (AbsTol(i)=0.0) was
890	! requested on a variable which has now vanished. The
891	! integration was successful as far as T.
892	!
893	! Note: Since the normal output value of ISTATE is 2, it
894	! does not need to be reset for normal continuation. Also,
895	! since a negative input value of ISTATE will be regarded
896	! as illegal, a negative output value requires the user to
897	! change it, and possibly other inputs, before calling the
898	! solver again.
899	!
900	! IOPT An integer flag to specify whether any optional inputs
901	! are being used on this call. Input only. The optional
902	! inputs are listed under a separate heading below.
903	! 0 No optional inputs are being used. Default values
904	! will be used in all cases.
905	! 1 One or more optional inputs are being used.
906	!
907	! RWORK A real working array (double precision). The length of
908	! RWORK must be at least
909	!
910	! 20 + NYH(MAXORD + 1) + 3NEQ + LWM
911	!
912	! where
913	! NYH = the initial value of NEQ,
914	! MAXORD = 12 (if METH = 1) or 5 (if METH = 2) (unless a
915	! smaller value is given as an optional input),
916	! LWM = 0 if MITER = 0,
917	! LWM = NEQ**2 + 2 if MITER = 1 or 2,
918	! LWM = NEQ + 2 if MITER = 3, and
919	! LWM = (2ML + MU + 1)NEQ + 2
920	! if MITER = 4 or 5.
921	! (See the MF description below for METH and MITER.)
922	!
923	! Thus if MAXORD has its default value and NEQ is constant,
924	! this length is:
925	! 20 + 16*NEQ for MF = 10,
926	! 22 + 16NEQ + NEQ*2 for MF = 11 or 12,
927	! 22 + 17*NEQ for MF = 13,
928	! 22 + 17NEQ + (2ML + MU)*NEQ for MF = 14 or 15,
929	! 20 + 9*NEQ for MF = 20,
930	! 22 + 9NEQ + NEQ*2 for MF = 21 or 22,
931	! 22 + 10*NEQ for MF = 23,
932	! 22 + 10NEQ + (2ML + MU)*NEQ for MF = 24 or 25.
933	!
934	! The first 20 words of RWORK are reserved for conditional
935	! and optional inputs and optional outputs.
936	!
937	! The following word in RWORK is a conditional input:
938	! RWORK(1) = TCRIT, the critical value of t which the
939	! solver is not to overshoot. Required if ITASK
940	! is 4 or 5, and ignored otherwise. See ITASK.
941	!
942	! LRW The length of the array RWORK, as declared by the user.
943	! (This will be checked by the solver.)
944	!
945	! IWORK An integer work array. Its length must be at least
946	! 20 if MITER = 0 or 3 (MF = 10, 13, 20, 23), or
947	! 20 + NEQ otherwise (MF = 11, 12, 14, 15, 21, 22, 24, 25).
948	! (See the MF description below for MITER.) The first few
949	! words of IWORK are used for conditional and optional
950	! inputs and optional outputs.
951	!
952	! The following two words in IWORK are conditional inputs:
953	! IWORK(1) = ML These are the lower and upper half-
954	! IWORK(2) = MU bandwidths, respectively, of the banded
955	! Jacobian, excluding the main diagonal.
956	! The band is defined by the matrix locations
957	! (i,j) with i - ML <= j <= i + MU. ML and MU
958	! must satisfy 0 <= ML,MU <= NEQ - 1. These are
959	! required if MITER is 4 or 5, and ignored
960	! otherwise. ML and MU may in fact be the band
961	! parameters for a matrix to which df/dy is only
962	! approximately equal.
963	!
964	! LIW The length of the array IWORK, as declared by the user.
965	! (This will be checked by the solver.)
966	!
967	! Note: The work arrays must not be altered between calls to DLSODE
968	! for the same problem, except possibly for the conditional and
969	! optional inputs, and except for the last 3*NEQ words of RWORK.
970	! The latter space is used for internal scratch space, and so is
971	! available for use by the user outside DLSODE between calls, if
972	! desired (but not for use by F or JAC).
973	!
974	! JAC The name of the user-supplied routine (MITER = 1 or 4) to
975	! compute the Jacobian matrix, df/dy, as a function of the
976	! scalar t and the vector y. (See the MF description below
977	! for MITER.) It is to have the form
978	!
979	! SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD)
980	! KPP_REAL T, Y(), PD(NROWPD,)
981	!
982	! where NEQ, T, Y, ML, MU, and NROWPD are input and the
983	! array PD is to be loaded with partial derivatives
984	! (elements of the Jacobian matrix) on output. PD must be
985	! given a first dimension of NROWPD. T and Y have the same
986	! meaning as in subroutine F.
987	!
988	! In the full matrix case (MITER = 1), ML and MU are
989	! ignored, and the Jacobian is to be loaded into PD in
990	! columnwise manner, with df(i)/dy(j) loaded into PD(i,j).
991	!
992	! In the band matrix case (MITER = 4), the elements within
993	! the band are to be loaded into PD in columnwise manner,
994	! with diagonal lines of df/dy loaded into the rows of PD.
995	! Thus df(i)/dy(j) is to be loaded into PD(i-j+MU+1,j). ML
996	! and MU are the half-bandwidth parameters (see IWORK).
997	! The locations in PD in the two triangular areas which
998	! correspond to nonexistent matrix elements can be ignored
999	! or loaded arbitrarily, as they are overwritten by DLSODE.
1000	!
1001	! JAC need not provide df/dy exactly. A crude approximation
1002	! (possibly with a smaller bandwidth) will do.
1003	!
1004	! In either case, PD is preset to zero by the solver, so
1005	! that only the nonzero elements need be loaded by JAC.
1006	! Each call to JAC is preceded by a call to F with the same
1007	! arguments NEQ, T, and Y. Thus to gain some efficiency,
1008	! intermediate quantities shared by both calculations may
1009	! be saved in a user COMMON block by F and not recomputed
1010	! by JAC, if desired. Also, JAC may alter the Y array, if
1011	! desired. JAC must be declared EXTERNAL in the calling
1012	! program.
1013	!
1014	! Subroutine JAC may access user-defined quantities in
1015	! NEQ(2),... and/or in Y(NEQ+1),... if NEQ is an array
1016	! (dimensioned in JAC) and/or Y has length exceeding
1017	! NEQ. See the descriptions of NEQ and Y above.
1018	!
1019	! MF The method flag. Used only for input. The legal values
1020	! of MF are 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24,
1021	! and 25. MF has decimal digits METH and MITER:
1022	! MF = 10*METH + MITER .
1023	!
1024	! METH indicates the basic linear multistep method:
1025	! 1 Implicit Adams method.
1026	! 2 Method based on backward differentiation formulas
1027	! (BDF's).
1028	!
1029	! MITER indicates the corrector iteration method:
1030	! 0 Functional iteration (no Jacobian matrix is
1031	! involved).
1032	! 1 Chord iteration with a user-supplied full (NEQ by
1033	! NEQ) Jacobian.
1034	! 2 Chord iteration with an internally generated
1035	! (difference quotient) full Jacobian (using NEQ
1036	! extra calls to F per df/dy value).
1037	! 3 Chord iteration with an internally generated
1038	! diagonal Jacobian approximation (using one extra call
1039	! to F per df/dy evaluation).
1040	! 4 Chord iteration with a user-supplied banded Jacobian.
1041	! 5 Chord iteration with an internally generated banded
1042	! Jacobian (using ML + MU + 1 extra calls to F per
1043	! df/dy evaluation).
1044	!
1045	! If MITER = 1 or 4, the user must supply a subroutine JAC
1046	! (the name is arbitrary) as described above under JAC.
1047	! For other values of MITER, a dummy argument can be used.
1048	!
1049	! Optional Inputs
1050	! ---------------
1051	! The following is a list of the optional inputs provided for in the
1052	! call sequence. (See also Part 2.) For each such input variable,
1053	! this table lists its name as used in this documentation, its
1054	! location in the call sequence, its meaning, and the default value.
1055	! The use of any of these inputs requires IOPT = 1, and in that case
1056	! all of these inputs are examined. A value of zero for any of
1057	! these optional inputs will cause the default value to be used.
1058	! Thus to use a subset of the optional inputs, simply preload
1059	! locations 5 to 10 in RWORK and IWORK to 0.0 and 0 respectively,
1060	! and then set those of interest to nonzero values.
1061	!
1062	! Name Location Meaning and default value
1063	! ------ --------- -----------------------------------------------
1064	! H0 RWORK(5) Step size to be attempted on the first step.
1065	! The default value is determined by the solver.
1066	! HMAX RWORK(6) Maximum absolute step size allowed. The
1067	! default value is infinite.
1068	! HMIN RWORK(7) Minimum absolute step size allowed. The
1069	! default value is 0. (This lower bound is not
1070	! enforced on the final step before reaching
1071	! TCRIT when ITASK = 4 or 5.)
1072	! MAXORD IWORK(5) Maximum order to be allowed. The default value
1073	! is 12 if METH = 1, and 5 if METH = 2. (See the
1074	! MF description above for METH.) If MAXORD
1075	! exceeds the default value, it will be reduced
1076	! to the default value. If MAXORD is changed
1077	! during the problem, it may cause the current
1078	! order to be reduced.
1079	! MXSTEP IWORK(6) Maximum number of (internally defined) steps
1080	! allowed during one call to the solver. The
1081	! default value is 500.
1082	! MXHNIL IWORK(7) Maximum number of messages printed (per
1083	! problem) warning that T + H = T on a step
1084	! (H = step size). This must be positive to
1085	! result in a nondefault value. The default
1086	! value is 10.
1087	!
1088	! Optional Outputs
1089	! ----------------
1090	! As optional additional output from DLSODE, the variables listed
1091	! below are quantities related to the performance of DLSODE which
1092	! are available to the user. These are communicated by way of the
1093	! work arrays, but also have internal mnemonic names as shown.
1094	! Except where stated otherwise, all of these outputs are defined on
1095	! any successful return from DLSODE, and on any return with ISTATE =
1096	! -1, -2, -4, -5, or -6. On an illegal input return (ISTATE = -3),
1097	! they will be unchanged from their existing values (if any), except
1098	! possibly for TOLSF, LENRW, and LENIW. On any error return,
1099	! outputs relevant to the error will be defined, as noted below.
1100	!
1101	! Name Location Meaning
1102	! ----- --------- ------------------------------------------------
1103	! HU RWORK(11) Step size in t last used (successfully).
1104	! HCUR RWORK(12) Step size to be attempted on the next step.
1105	! TCUR RWORK(13) Current value of the independent variable which
1106	! the solver has actually reached, i.e., the
1107	! current internal mesh point in t. On output,
1108	! TCUR will always be at least as far as the
1109	! argument T, but may be farther (if interpolation
1110	! was done).
1111	! TOLSF RWORK(14) Tolerance scale factor, greater than 1.0,
1112	! computed when a request for too much accuracy
1113	! was detected (ISTATE = -3 if detected at the
1114	! start of the problem, ISTATE = -2 otherwise).
1115	! If ITOL is left unaltered but RelTol and AbsTol are
1116	! uniformly scaled up by a factor of TOLSF for the
1117	! next call, then the solver is deemed likely to
1118	! succeed. (The user may also ignore TOLSF and
1119	! alter the tolerance parameters in any other way
1120	! appropriate.)
1121	! NST IWORK(11) Number of steps taken for the problem so far.
1122	! NFE IWORK(12) Number of F evaluations for the problem so far.
1123	! NJE IWORK(13) Number of Jacobian evaluations (and of matrix LU
1124	! decompositions) for the problem so far.
1125	! NQU IWORK(14) Method order last used (successfully).
1126	! NQCUR IWORK(15) Order to be attempted on the next step.
1127	! IMXER IWORK(16) Index of the component of largest magnitude in
1128	! the weighted local error vector ( e(i)/EWT(i) ),
1129	! on an error return with ISTATE = -4 or -5.
1130	! LENRW IWORK(17) Length of RWORK actually required. This is
1131	! defined on normal returns and on an illegal
1132	! input return for insufficient storage.
1133	! LENIW IWORK(18) Length of IWORK actually required. This is
1134	! defined on normal returns and on an illegal
1135	! input return for insufficient storage.
1136	!
1137	! The following two arrays are segments of the RWORK array which may
1138	! also be of interest to the user as optional outputs. For each
1139	! array, the table below gives its internal name, its base address
1140	! in RWORK, and its description.
1141	!
1142	! Name Base address Description
1143	! ---- ------------ ----------------------------------------------
1144	! YH 21 The Nordsieck history array, of size NYH by
1145	! (NQCUR + 1), where NYH is the initial value of
1146	! NEQ. For j = 0,1,...,NQCUR, column j + 1 of
1147	! YH contains HCUR**j/factorial(j) times the jth
1148	! derivative of the interpolating polynomial
1149	! currently representing the solution, evaluated
1150	! at t = TCUR.
1151	! ACOR LENRW-NEQ+1 Array of size NEQ used for the accumulated
1152	! corrections on each step, scaled on output to
1153	! represent the estimated local error in Y on
1154	! the last step. This is the vector e in the
1155	! description of the error control. It is
1156	! defined only on successful return from DLSODE.
1157	!
1158	!
1159	! Part 2. Other Callable Routines
1160	! --------------------------------
1161	!
1162	! The following are optional calls which the user may make to gain
1163	! additional capabilities in conjunction with DLSODE.
1164	!
1165	! Form of call Function
1166	! ------------------------ ----------------------------------------
1167	! CALL XSETUN(LUN) Set the logical unit number, LUN, for
1168	! output of messages from DLSODE, if the
1169	! default is not desired. The default
1170	! value of LUN is 6. This call may be made
1171	! at any time and will take effect
1172	! immediately.
1173	! CALL XSETF(MFLAG) Set a flag to control the printing of
1174	! messages by DLSODE. MFLAG = 0 means do
1175	! not print. (Danger: this risks losing
1176	! valuable information.) MFLAG = 1 means
1177	! print (the default). This call may be
1178	! made at any time and will take effect
1179	! immediately.
1180	! CALL DSRCOM(RSAV,ISAV,JOB) Saves and restores the contents of the
1181	! internal COMMON blocks used by DLSODE
1182	! (see Part 3 below). RSAV must be a
1183	! real array of length 218 or more, and
1184	! ISAV must be an integer array of length
1185	! 37 or more. JOB = 1 means save COMMON
1186	! into RSAV/ISAV. JOB = 2 means restore
1187	! COMMON from same. DSRCOM is useful if
1188	! one is interrupting a run and restarting
1189	! later, or alternating between two or
1190	! more problems solved with DLSODE.
1191	! CALL DINTDY(,,,,,) Provide derivatives of y, of various
1192	! (see below) orders, at a specified point t, if
1193	! desired. It may be called only after a
1194	! successful return from DLSODE. Detailed
1195	! instructions follow.
1196	!
1197	! Detailed instructions for using DINTDY
1198	! --------------------------------------
1199	! The form of the CALL is:
1200	!
1201	! CALL DINTDY (T, K, RWORK(21), NYH, DKY, IFLAG)
1202	!
1203	! The input parameters are:
1204	!
1205	! T Value of independent variable where answers are
1206	! desired (normally the same as the T last returned by
1207	! DLSODE). For valid results, T must lie between
1208	! TCUR - HU and TCUR. (See "Optional Outputs" above
1209	! for TCUR and HU.)
1210	! K Integer order of the derivative desired. K must
1211	! satisfy 0 <= K <= NQCUR, where NQCUR is the current
1212	! order (see "Optional Outputs"). The capability
1213	! corresponding to K = 0, i.e., computing y(t), is
1214	! already provided by DLSODE directly. Since
1215	! NQCUR >= 1, the first derivative dy/dt is always
1216	! available with DINTDY.
1217	! RWORK(21) The base address of the history array YH.
1218	! NYH Column length of YH, equal to the initial value of NEQ.
1219	!
1220	! The output parameters are:
1221	!
1222	! DKY Real array of length NEQ containing the computed value
1223	! of the Kth derivative of y(t).
1224	! IFLAG Integer flag, returned as 0 if K and T were legal,
1225	! -1 if K was illegal, and -2 if T was illegal.
1226	! On an error return, a message is also written.
1227	!
1228	!
1229	! Part 3. Common Blocks
1230	! ----------------------
1231	!
1232	! If DLSODE is to be used in an overlay situation, the user must
1233	! declare, in the primary overlay, the variables in:
1234	! (1) the call sequence to DLSODE,
1235	! (2) the internal COMMON block /DLS001/, of length 255
1236	! (218 double precision words followed by 37 integer words).
1237	!
1238	! If DLSODE is used on a system in which the contents of internal
1239	! COMMON blocks are not preserved between calls, the user should
1240	! declare the above COMMON block in his main program to insure that
1241	! its contents are preserved.
1242	!
1243	! If the solution of a given problem by DLSODE is to be interrupted
1244	! and then later continued, as when restarting an interrupted run or
1245	! alternating between two or more problems, the user should save,
1246	! following the return from the last DLSODE call prior to the
1247	! interruption, the contents of the call sequence variables and the
1248	! internal COMMON block, and later restore these values before the
1249	! next DLSODE call for that problem. In addition, if XSETUN and/or
1250	! XSETF was called for non-default handling of error messages, then
1251	! these calls must be repeated. To save and restore the COMMON
1252	! block, use subroutine DSRCOM (see Part 2 above).
1253	!
1254	!
1255	! Part 4. Optionally Replaceable Solver Routines
1256	! -----------------------------------------------
1257	!
1258	! Below are descriptions of two routines in the DLSODE package which
1259	! relate to the measurement of errors. Either routine can be
1260	! replaced by a user-supplied version, if desired. However, since
1261	! such a replacement may have a major impact on performance, it
1262	! should be done only when absolutely necessary, and only with great
1263	! caution. (Note: The means by which the package version of a
1264	! routine is superseded by the user's version may be system-
1265	! dependent.)
1266	!
1267	! DEWSET
1268	! ------
1269	! The following subroutine is called just before each internal
1270	! integration step, and sets the array of error weights, EWT, as
1271	! described under ITOL/RelTol/AbsTol above:
1272	!
1273	! SUBROUTINE DEWSET (NEQ, ITOL, RelTol, AbsTol, YCUR, EWT)
1274	!
1275	! where NEQ, ITOL, RelTol, and AbsTol are as in the DLSODE call
1276	! sequence, YCUR contains the current dependent variable vector,
1277	! and EWT is the array of weights set by DEWSET.
1278	!
1279	! If the user supplies this subroutine, it must return in EWT(i)
1280	! (i = 1,...,NEQ) a positive quantity suitable for comparing errors
1281	! in Y(i) to. The EWT array returned by DEWSET is passed to the
1282	! DVNORM routine (see below), and also used by DLSODE in the
1283	! computation of the optional output IMXER, the diagonal Jacobian
1284	! approximation, and the increments for difference quotient
1285	! Jacobians.
1286	!
1287	! In the user-supplied version of DEWSET, it may be desirable to use
1288	! the current values of derivatives of y. Derivatives up to order NQ
1289	! are available from the history array YH, described above under
1290	! optional outputs. In DEWSET, YH is identical to the YCUR array,
1291	! extended to NQ + 1 columns with a column length of NYH and scale
1292	! factors of H**j/factorial(j). On the first call for the problem,
1293	! given by NST = 0, NQ is 1 and H is temporarily set to 1.0.
1294	! NYH is the initial value of NEQ. The quantities NQ, H, and NST
1295	! can be obtained by including in SEWSET the statements:
1296	! KPP_REAL RLS
1297	! COMMON /DLS001/ RLS(218),ILS(37)
1298	! NQ = ILS(33)
1299	! NST = ILS(34)
1300	! H = RLS(212)
1301	! Thus, for example, the current value of dy/dt can be obtained as
1302	! YCUR(NYH+i)/H (i=1,...,NEQ) (and the division by H is unnecessary
1303	! when NST = 0).
1304	!
1305	! DVNORM
1306	! ------
1307	! DVNORM is a real function routine which computes the weighted
1308	! root-mean-square norm of a vector v:
1309	!
1310	! d = DVNORM (n, v, w)
1311	!
1312	! where:
1313	! n = the length of the vector,
1314	! v = real array of length n containing the vector,
1315	! w = real array of length n containing weights,
1316	! d = SQRT( (1/n) * sum(v(i)w(i))*2 ).
1317	!
1318	! DVNORM is called with n = NEQ and with w(i) = 1.0/EWT(i), where
1319	! EWT is as set by subroutine DEWSET.
1320	!
1321	! If the user supplies this function, it should return a nonnegative
1322	! value of DVNORM suitable for use in the error control in DLSODE.
1323	! None of the arguments should be altered by DVNORM. For example, a
1324	! user-supplied DVNORM routine might:
1325	! - Substitute a max-norm of (v(i)*w(i)) for the rms-norm, or
1326	! - Ignore some components of v in the norm, with the effect of
1327	! suppressing the error control on those components of Y.
1328	! ---------------------------------------------------------------------
1329	!***ROUTINES CALLED DEWSET, DINTDY, DUMACH, DSTODE, DVNORM, XERRWD
1330	!***COMMON BLOCKS DLS001
1331	!***REVISION HISTORY (YYYYMMDD)
1332	! 19791129 DATE WRITTEN
1333	! 19791213 Minor changes to declarations; DELP init. in STODE.
1334	! 19800118 Treat NEQ as array; integer declarations added throughout;
1335	! minor changes to prologue.
1336	! 19800306 Corrected TESCO(1,NQP1) setting in CFODE.
1337	! 19800519 Corrected access of YH on forced order reduction;
1338	! numerous corrections to prologues and other comments.
1339	! 19800617 In main driver, added loading of SQRT(UROUND) in RWORK;
1340	! minor corrections to main prologue.
1341	! 19800923 Added zero initialization of HU and NQU.
1342	! 19801218 Revised XERRWD routine; minor corrections to main prologue.
1343	! 19810401 Minor changes to comments and an error message.
1344	! 19810814 Numerous revisions: replaced EWT by 1/EWT; used flags
1345	! JCUR, ICF, IERPJ, IERSL between STODE and subordinates;
1346	! added tuning parameters CCMAX, MAXCOR, MSBP, MXNCF;
1347	! reorganized returns from STODE; reorganized type decls.;
1348	! fixed message length in XERRWD; changed default LUNIT to 6;
1349	! changed Common lengths; changed comments throughout.
1350	! 19870330 Major update by ACH: corrected comments throughout;
1351	! removed TRET from Common; rewrote EWSET with 4 loops;
1352	! fixed t test in INTDY; added Cray directives in STODE;
1353	! in STODE, fixed DELP init. and logic around PJAC call;
1354	! combined routines to save/restore Common;
1355	! passed LEVEL = 0 in error message calls (except run abort).
1356	! 19890426 Modified prologue to SLATEC/LDOC format. (FNF)
1357	! 19890501 Many improvements to prologue. (FNF)
1358	! 19890503 A few final corrections to prologue. (FNF)
1359	! 19890504 Minor cosmetic changes. (FNF)
1360	! 19890510 Corrected description of Y in Arguments section. (FNF)
1361	! 19890517 Minor corrections to prologue. (FNF)
1362	! 19920514 Updated with prologue edited 891025 by G. Shaw for manual.
1363	! 19920515 Converted source lines to upper case. (FNF)
1364	! 19920603 Revised XERRWD calls using mixed upper-lower case. (ACH)
1365	! 19920616 Revised prologue comment regarding CFT. (ACH)
1366	! 19921116 Revised prologue comments regarding Common. (ACH).
1367	! 19930326 Added comment about non-reentrancy. (FNF)
1368	! 19930723 Changed D1MACH to DUMACH. (FNF)
1369	! 19930801 Removed ILLIN and NTREP from Common (affects driver logic);
1370	! minor changes to prologue and internal comments;
1371	! changed Hollerith strings to quoted strings;
1372	! changed internal comments to mixed case;
1373	! replaced XERRWD with new version using character type;
1374	! changed dummy dimensions from 1 to *. (ACH)
1375	! 19930809 Changed to generic intrinsic names; changed names of
1376	! subprograms and Common blocks to DLSODE etc. (ACH)
1377	! 19930929 Eliminated use of REAL intrinsic; other minor changes. (ACH)
1378	! 20010412 Removed all 'own' variables from Common block /DLS001/
1379	! (affects declarations in 6 routines). (ACH)
1380	! 20010509 Minor corrections to prologue. (ACH)
1381	! 20031105 Restored 'own' variables to Common block /DLS001/, to
1382	! enable interrupt/restart feature. (ACH)
1383	! 20031112 Added SAVE statements for data-loaded constants.
1384	!
1385	!***END PROLOGUE DLSODE
1386	!
1387	!*Internal Notes:
1388	!
1389	! Other Routines in the DLSODE Package.
1390	!
1391	! In addition to Subroutine DLSODE, the DLSODE package includes the
1392	! following subroutines and function routines:
1393	! DINTDY computes an interpolated value of the y vector at t = TOUT.
1394	! DSTODE is the core integrator, which does one step of the
1395	! integration and the associated error control.
1396	! DCFODE sets all method coefficients and test constants.
1397	! DPREPJ computes and preprocesses the Jacobian matrix J = df/dy
1398	! and the Newton iteration matrix P = I - hl0J.
1399	! DSOLSY manages solution of linear system in chord iteration.
1400	! DEWSET sets the error weight vector EWT before each step.
1401	! DVNORM computes the weighted R.M.S. norm of a vector.
1402	! DSRCOM is a user-callable routine to save and restore
1403	! the contents of the internal Common block.
1404	! DGEFA and DGESL are routines from LINPACK for solving full
1405	! systems of linear algebraic equations.
1406	! DGBFA and DGBSL are routines from LINPACK for solving banded
1407	! linear systems.
1408	! DUMACH computes the unit roundoff in a machine-independent manner.
1409	! XERRWD, XSETUN, XSETF, IXSAV, IUMACH handle the printing of all
1410	! error messages and warnings. XERRWD is machine-dependent.
1411	! Note: DVNORM, DUMACH, IXSAV, and IUMACH are function routines.
1412	! All the others are subroutines.
1413	!
1414	!**End
1415	!
1416	! Declare externals.
1417	! Note: they are now internal
1418	!EXTERNAL DPREPJ, DSOLSY
1419	!KPP_REAL DUMACH, DVNORM
1420	!
1421	! Declare all other variables.
1422	INTEGER INIT, MXSTEP, MXHNIL, NHNIL, NSLAST, NYH, IOWNS, &
1423	ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, &
1424	LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, &
1425	MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
1426	INTEGER I, I1, I2, IFLAG, IMXER, KGO, LF0, &
1427	LENIW, LENRW, LENWM, ML, MORD(2), MU, MXHNL0, MXSTP0
1428	KPP_REAL ROWNS, &
1429	CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
1430	KPP_REAL AbsTolI, AYI, BIG, EWTI, H0, HMAX, HMX, RH, RelTolI, &
1431	TCRIT, TDIST, TNEXT, TOL, TOLSF, TP, SIZE, SUM, W0
1432
1433	LOGICAL IHIT
1434	CHARACTER*80 MSG
1435	SAVE MORD, MXSTP0, MXHNL0
1436	!-----------------------------------------------------------------------
1437	! The following internal Common block contains
1438	! (a) variables which are local to any subroutine but whose values must
1439	! be preserved between calls to the routine ("own" variables), and
1440	! (b) variables which are communicated between subroutines.
1441	! The block DLS001 is declared in subroutines DLSODE, DINTDY, DSTODE,
1442	! DPREPJ, and DSOLSY.
1443	! Groups of variables are replaced by dummy arrays in the Common
1444	! declarations in routines where those variables are not used.
1445	!-----------------------------------------------------------------------
1446	COMMON /DLS001/ ROWNS(209), &
1447	CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, &
1448	INIT, MXSTEP, MXHNIL, NHNIL, NSLAST, NYH, IOWNS(6), &
1449	ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, &
1450	LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, &
1451	MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
1452	!
1453	DATA MORD(1),MORD(2)/12,5/, MXSTP0/5000/, MXHNL0/10/
1454	!-----------------------------------------------------------------------
1455	! Block A.
1456	! This code block is executed on every call.
1457	! It tests ISTATE and ITASK for legality and branches appropriately.
1458	! If ISTATE .GT. 1 but the flag INIT shows that initialization has
1459	! not yet been done, an error return occurs.
1460	! If ISTATE = 1 and TOUT = T, return immediately.
1461	!-----------------------------------------------------------------------
1462	!
1463	!***FIRST EXECUTABLE STATEMENT DLSODE
1464	IF (ISTATE .LT. 1 .OR. ISTATE .GT. 3) GO TO 601
1465	IF (ITASK .LT. 1 .OR. ITASK .GT. 5) GO TO 602
1466	IF (ISTATE .EQ. 1) GO TO 10
1467	IF (INIT .EQ. 0) GO TO 603
1468	IF (ISTATE .EQ. 2) GO TO 200
1469	GO TO 20
1470	10 INIT = 0
1471	IF (TOUT .EQ. T) RETURN
1472	!-----------------------------------------------------------------------
1473	! Block B.
1474	! The next code block is executed for the initial call (ISTATE = 1),
1475	! or for a continuation call with parameter changes (ISTATE = 3).
1476	! It contains checking of all inputs and various initializations.
1477	!
1478	! First check legality of the non-optional inputs NEQ, ITOL, IOPT,
1479	! MF, ML, and MU.
1480	!-----------------------------------------------------------------------
1481	20 IF (NEQ .LE. 0) GO TO 604
1482	IF (ISTATE .EQ. 1) GO TO 25
1483	IF (NEQ .GT. N) GO TO 605
1484	25 N = NEQ
1485	IF (ITOL .LT. 1 .OR. ITOL .GT. 4) GO TO 606
1486	IF (IOPT .LT. 0 .OR. IOPT .GT. 1) GO TO 607
1487	METH = MF/10
1488	MITER = MF - 10*METH
1489	IF (METH .LT. 1 .OR. METH .GT. 2) GO TO 608
1490	IF (MITER .LT. 0 .OR. MITER .GT. 5) GO TO 608
1491	IF (MITER .LE. 3) GO TO 30
1492	ML = IWORK(1)
1493	MU = IWORK(2)
1494	IF (ML .LT. 0 .OR. ML .GE. N) GO TO 609
1495	IF (MU .LT. 0 .OR. MU .GE. N) GO TO 610
1496	30 CONTINUE
1497	! Next process and check the optional inputs. --------------------------
1498	IF (IOPT .EQ. 1) GO TO 40
1499	MAXORD = MORD(METH)
1500	MXSTEP = MXSTP0
1501	MXHNIL = MXHNL0
1502	IF (ISTATE .EQ. 1) H0 = 0.0D0
1503	HMXI = 0.0D0
1504	HMIN = 0.0D0
1505	GO TO 60
1506	40 MAXORD = IWORK(5)
1507	IF (MAXORD .LT. 0) GO TO 611
1508	IF (MAXORD .EQ. 0) MAXORD = 100
1509	MAXORD = MIN(MAXORD,MORD(METH))
1510	MXSTEP = IWORK(6)
1511	IF (MXSTEP .LT. 0) GO TO 612
1512	IF (MXSTEP .EQ. 0) MXSTEP = MXSTP0
1513	MXHNIL = IWORK(7)
1514	IF (MXHNIL .LT. 0) GO TO 613
1515	IF (MXHNIL .EQ. 0) MXHNIL = MXHNL0
1516	IF (ISTATE .NE. 1) GO TO 50
1517	H0 = RWORK(5)
1518	IF ((TOUT - T)*H0 .LT. 0.0D0) GO TO 614
1519	50 HMAX = RWORK(6)
1520	IF (HMAX .LT. 0.0D0) GO TO 615
1521	HMXI = 0.0D0
1522	IF (HMAX .GT. 0.0D0) HMXI = 1.0D0/HMAX
1523	HMIN = RWORK(7)
1524	IF (HMIN .LT. 0.0D0) GO TO 616
1525	!-----------------------------------------------------------------------
1526	! Set work array pointers and check lengths LRW and LIW.
1527	! Pointers to segments of RWORK and IWORK are named by prefixing L to
1528	! the name of the segment. E.g., the segment YH starts at RWORK(LYH).
1529	! Segments of RWORK (in order) are denoted YH, WM, EWT, SAVF, ACOR.
1530	!-----------------------------------------------------------------------
1531	60 LYH = 21
1532	IF (ISTATE .EQ. 1) NYH = N
1533	LWM = LYH + (MAXORD + 1)*NYH
1534	IF (MITER .EQ. 0) LENWM = 0
1535	IF (MITER .EQ. 1 .OR. MITER .EQ. 2) LENWM = N*N + 2
1536	IF (MITER .EQ. 3) LENWM = N + 2
1537	IF (MITER .GE. 4) LENWM = (2ML + MU + 1)N + 2
1538	LEWT = LWM + LENWM
1539	LSAVF = LEWT + N
1540	LACOR = LSAVF + N
1541	LENRW = LACOR + N - 1
1542	IWORK(17) = LENRW
1543	LIWM = 1
1544	LENIW = 20 + N
1545	IF (MITER .EQ. 0 .OR. MITER .EQ. 3) LENIW = 20
1546	IWORK(18) = LENIW
1547	IF (LENRW .GT. LRW) GO TO 617
1548	IF (LENIW .GT. LIW) GO TO 618
1549	! Check RelTol and AbsTol for legality. ------------------------------------
1550	RelTolI = RelTol(1)
1551	AbsTolI = AbsTol(1)
1552	DO 70 I = 1,N
1553	IF (ITOL .GE. 3) RelTolI = RelTol(I)
1554	IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) AbsTolI = AbsTol(I)
1555	IF (RelTolI .LT. 0.0D0) GO TO 619
1556	IF (AbsTolI .LT. 0.0D0) GO TO 620
1557	70 CONTINUE
1558	IF (ISTATE .EQ. 1) GO TO 100
1559	! If ISTATE = 3, set flag to signal parameter changes to DSTODE. -------
1560	JSTART = -1
1561	IF (NQ .LE. MAXORD) GO TO 90
1562	! MAXORD was reduced below NQ. Copy YH(*,MAXORD+2) into SAVF. ---------
1563	DO 80 I = 1,N
1564	80 RWORK(I+LSAVF-1) = RWORK(I+LWM-1)
1565	! Reload WM(1) = RWORK(LWM), since LWM may have changed. ---------------
1566	90 IF (MITER .GT. 0) RWORK(LWM) = SQRT(UROUND)
1567	IF (N .EQ. NYH) GO TO 200
1568	! NEQ was reduced. Zero part of YH to avoid undefined references. -----
1569	I1 = LYH + L*NYH
1570	I2 = LYH + (MAXORD + 1)*NYH - 1
1571	IF (I1 .GT. I2) GO TO 200
1572	DO 95 I = I1,I2
1573	95 RWORK(I) = 0.0D0
1574	GO TO 200
1575	!-----------------------------------------------------------------------
1576	! Block C.
1577	! The next block is for the initial call only (ISTATE = 1).
1578	! It contains all remaining initializations, the initial call to F,
1579	! and the calculation of the initial step size.
1580	! The error weights in EWT are inverted after being loaded.
1581	!-----------------------------------------------------------------------
1582	100 UROUND = DUMACH()
1583	TN = T
1584	IF (ITASK .NE. 4 .AND. ITASK .NE. 5) GO TO 110
1585	TCRIT = RWORK(1)
1586	IF ((TCRIT - TOUT)*(TOUT - T) .LT. 0.0D0) GO TO 625
1587	IF (H0 .NE. 0.0D0 .AND. (T + H0 - TCRIT)*H0 .GT. 0.0D0) &
1588	H0 = TCRIT - T
1589	110 JSTART = 0
1590	IF (MITER .GT. 0) RWORK(LWM) = SQRT(UROUND)
1591	NHNIL = 0
1592	NST = 0
1593	NJE = 0
1594	NSLAST = 0
1595	HU = 0.0D0
1596	NQU = 0
1597	CCMAX = 0.3D0
1598	MAXCOR = 3
1599	MSBP = 20
1600	MXNCF = 10
1601	! Initial call to F. (LF0 points to YH(*,2).) -------------------------
1602	LF0 = LYH + NYH
1603	CALL F (NEQ, T, Y, RWORK(LF0))
1604	NFE = 1
1605	! Load the initial value vector in YH. ---------------------------------
1606	DO 115 I = 1,N
1607	115 RWORK(I+LYH-1) = Y(I)
1608	! Load and invert the EWT array. (H is temporarily set to 1.0.) -------
1609	NQ = 1
1610	H = 1.0D0
1611	CALL DEWSET (N, ITOL, RelTol, AbsTol, RWORK(LYH), RWORK(LEWT))
1612	DO 120 I = 1,N
1613	IF (RWORK(I+LEWT-1) .LE. 0.0D0) GO TO 621
1614	120 RWORK(I+LEWT-1) = 1.0D0/RWORK(I+LEWT-1)
1615	!-----------------------------------------------------------------------
1616	! The coding below computes the step size, H0, to be attempted on the
1617	! first step, unless the user has supplied a value for this.
1618	! First check that TOUT - T differs significantly from zero.
1619	! A scalar tolerance quantity TOL is computed, as MAX(RelTol(I))
1620	! if this is positive, or MAX(AbsTol(I)/ABS(Y(I))) otherwise, adjusted
1621	! so as to be between 100*UROUND and 1.0E-3.
1622	! Then the computed value H0 is given by..
1623	! NEQ
1624	! H02 = TOL / ( w0-2 + (1/NEQ) * SUM ( f(i)/ywt(i) )**2 )
1625	! 1
1626	! where w0 = MAX ( ABS(T), ABS(TOUT) ),
1627	! f(i) = i-th component of initial value of f,
1628	! ywt(i) = EWT(i)/TOL (a weight for y(i)).
1629	! The sign of H0 is inferred from the initial values of TOUT and T.
1630	!-----------------------------------------------------------------------
1631	IF (H0 .NE. 0.0D0) GO TO 180
1632	TDIST = ABS(TOUT - T)
1633	W0 = MAX(ABS(T),ABS(TOUT))
1634	IF (TDIST .LT. 2.0D0UROUNDW0) GO TO 622
1635	TOL = RelTol(1)
1636	IF (ITOL .LE. 2) GO TO 140
1637	DO 130 I = 1,N
1638	130 TOL = MAX(TOL,RelTol(I))
1639	140 IF (TOL .GT. 0.0D0) GO TO 160
1640	AbsTolI = AbsTol(1)
1641	DO 150 I = 1,N
1642	IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) AbsTolI = AbsTol(I)
1643	AYI = ABS(Y(I))
1644	IF (AYI .NE. 0.0D0) TOL = MAX(TOL,AbsTolI/AYI)
1645	150 CONTINUE
1646	160 TOL = MAX(TOL,100.0D0*UROUND)
1647	TOL = MIN(TOL,0.001D0)
1648	SUM = DVNORM (N, RWORK(LF0), RWORK(LEWT))
1649	SUM = 1.0D0/(TOLW0W0) + TOLSUM*2
1650	H0 = 1.0D0/SQRT(SUM)
1651	H0 = MIN(H0,TDIST)
1652	H0 = SIGN(H0,TOUT-T)
1653	! Adjust H0 if necessary to meet HMAX bound. ---------------------------
1654	180 RH = ABS(H0)*HMXI
1655	IF (RH .GT. 1.0D0) H0 = H0/RH
1656	! Load H with H0 and scale YH(*,2) by H0. ------------------------------
1657	H = H0
1658	DO 190 I = 1,N
1659	190 RWORK(I+LF0-1) = H0*RWORK(I+LF0-1)
1660	GO TO 270
1661	!-----------------------------------------------------------------------
1662	! Block D.
1663	! The next code block is for continuation calls only (ISTATE = 2 or 3)
1664	! and is to check stop conditions before taking a step.
1665	!-----------------------------------------------------------------------
1666	200 NSLAST = NST
1667	GO TO (210, 250, 220, 230, 240), ITASK
1668	210 IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 250
1669	CALL DINTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG)
1670	IF (IFLAG .NE. 0) GO TO 627
1671	T = TOUT
1672	GO TO 420
1673	220 TP = TN - HU(1.0D0 + 100.0D0UROUND)
1674	IF ((TP - TOUT)*H .GT. 0.0D0) GO TO 623
1675	IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 250
1676	GO TO 400
1677	230 TCRIT = RWORK(1)
1678	IF ((TN - TCRIT)*H .GT. 0.0D0) GO TO 624
1679	IF ((TCRIT - TOUT)*H .LT. 0.0D0) GO TO 625
1680	IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 245
1681	CALL DINTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG)
1682	IF (IFLAG .NE. 0) GO TO 627
1683	T = TOUT
1684	GO TO 420
1685	240 TCRIT = RWORK(1)
1686	IF ((TN - TCRIT)*H .GT. 0.0D0) GO TO 624
1687	245 HMX = ABS(TN) + ABS(H)
1688	IHIT = ABS(TN - TCRIT) .LE. 100.0D0UROUNDHMX
1689	IF (IHIT) GO TO 400
1690	TNEXT = TN + H(1.0D0 + 4.0D0UROUND)
1691	IF ((TNEXT - TCRIT)*H .LE. 0.0D0) GO TO 250
1692	H = (TCRIT - TN)(1.0D0 - 4.0D0UROUND)
1693	IF (ISTATE .EQ. 2) JSTART = -2
1694	!-----------------------------------------------------------------------
1695	! Block E.
1696	! The next block is normally executed for all calls and contains
1697	! the call to the one-step core integrator DSTODE.
1698	!
1699	! This is a looping point for the integration steps.
1700	!
1701	! First check for too many steps being taken, update EWT (if not at
1702	! start of problem), check for too much accuracy being requested, and
1703	! check for H below the roundoff level in T.
1704	!-----------------------------------------------------------------------
1705	250 CONTINUE
1706	IF ((NST-NSLAST) .GE. MXSTEP) GO TO 500
1707	CALL DEWSET (N, ITOL, RelTol, AbsTol, RWORK(LYH), RWORK(LEWT))
1708	DO 260 I = 1,N
1709	IF (RWORK(I+LEWT-1) .LE. 0.0D0) GO TO 510
1710	260 RWORK(I+LEWT-1) = 1.0D0/RWORK(I+LEWT-1)
1711	270 TOLSF = UROUND*DVNORM (N, RWORK(LYH), RWORK(LEWT))
1712	IF (TOLSF .LE. 1.0D0) GO TO 280
1713	TOLSF = TOLSF*2.0D0
1714	IF (NST .EQ. 0) GO TO 626
1715	GO TO 520
1716	280 IF ((TN + H) .NE. TN) GO TO 290
1717	NHNIL = NHNIL + 1
1718	IF (NHNIL .GT. MXHNIL) GO TO 290
1719	MSG = 'DLSODE- Warning..internal T (=R1) and H (=R2) are'
1720	CALL XERRWD (MSG, 50, 101, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
1721	MSG=' such that in the machine, T + H = T on the next step '
1722	CALL XERRWD (MSG, 60, 101, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
1723	MSG = ' (H = step size). Solver will continue anyway'
1724	CALL XERRWD (MSG, 50, 101, 0, 0, 0, 0, 2, TN, H)
1725	IF (NHNIL .LT. MXHNIL) GO TO 290
1726	MSG = 'DLSODE- Above warning has been issued I1 times. '
1727	CALL XERRWD (MSG, 50, 102, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
1728	MSG = ' It will not be issued again for this problem'
1729	CALL XERRWD (MSG, 50, 102, 0, 1, MXHNIL, 0, 0, 0.0D0, 0.0D0)
1730	290 CONTINUE
1731	!-----------------------------------------------------------------------
1732	! CALL DSTODE(NEQ,Y,YH,NYH,YH,EWT,SAVF,ACOR,WM,IWM,F,JAC,DPREPJ,DSOLSY)
1733	!-----------------------------------------------------------------------
1734	CALL DSTODE (NEQ, Y, RWORK(LYH), NYH, RWORK(LYH), RWORK(LEWT), &
1735	RWORK(LSAVF), RWORK(LACOR), RWORK(LWM), IWORK(LIWM), &
1736	F, JAC)
1737	!F, JAC, DPREPJ, DSOLSY)
1738	KGO = 1 - KFLAG
1739	GO TO (300, 530, 540), KGO
1740	!-----------------------------------------------------------------------
1741	! Block F.
1742	! The following block handles the case of a successful return from the
1743	! core integrator (KFLAG = 0). Test for stop conditions.
1744	!-----------------------------------------------------------------------
1745	300 INIT = 1
1746	GO TO (310, 400, 330, 340, 350), ITASK
1747	! ITASK = 1. If TOUT has been reached, interpolate. -------------------
1748	310 IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 250
1749	CALL DINTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG)
1750	T = TOUT
1751	GO TO 420
1752	! ITASK = 3. Jump to exit if TOUT was reached. ------------------------
1753	330 IF ((TN - TOUT)*H .GE. 0.0D0) GO TO 400
1754	GO TO 250
1755	! ITASK = 4. See if TOUT or TCRIT was reached. Adjust H if necessary.
1756	340 IF ((TN - TOUT)*H .LT. 0.0D0) GO TO 345
1757	CALL DINTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG)
1758	T = TOUT
1759	GO TO 420
1760	345 HMX = ABS(TN) + ABS(H)
1761	IHIT = ABS(TN - TCRIT) .LE. 100.0D0UROUNDHMX
1762	IF (IHIT) GO TO 400
1763	TNEXT = TN + H(1.0D0 + 4.0D0UROUND)
1764	IF ((TNEXT - TCRIT)*H .LE. 0.0D0) GO TO 250
1765	H = (TCRIT - TN)(1.0D0 - 4.0D0UROUND)
1766	JSTART = -2
1767	GO TO 250
1768	! ITASK = 5. See if TCRIT was reached and jump to exit. ---------------
1769	350 HMX = ABS(TN) + ABS(H)
1770	IHIT = ABS(TN - TCRIT) .LE. 100.0D0UROUNDHMX
1771	!-----------------------------------------------------------------------
1772	! Block G.
1773	! The following block handles all successful returns from DLSODE.
1774	! If ITASK .NE. 1, Y is loaded from YH and T is set accordingly.
1775	! ISTATE is set to 2, and the optional outputs are loaded into the
1776	! work arrays before returning.
1777	!-----------------------------------------------------------------------
1778	400 DO 410 I = 1,N
1779	410 Y(I) = RWORK(I+LYH-1)
1780	T = TN
1781	IF (ITASK .NE. 4 .AND. ITASK .NE. 5) GO TO 420
1782	IF (IHIT) T = TCRIT
1783	420 ISTATE = 2
1784	RWORK(11) = HU
1785	RWORK(12) = H
1786	RWORK(13) = TN
1787	IWORK(11) = NST
1788	IWORK(12) = NFE
1789	IWORK(13) = NJE
1790	IWORK(14) = NQU
1791	IWORK(15) = NQ
1792	RETURN
1793	!-----------------------------------------------------------------------
1794	! Block H.
1795	! The following block handles all unsuccessful returns other than
1796	! those for illegal input. First the error message routine is called.
1797	! If there was an error test or convergence test failure, IMXER is set.
1798	! Then Y is loaded from YH and T is set to TN. The optional outputs
1799	! are loaded into the work arrays before returning.
1800	!-----------------------------------------------------------------------
1801	! The maximum number of steps was taken before reaching TOUT. ----------
1802	500 MSG = 'DLSODE- At current T (=R1), MXSTEP (=I1) steps '
1803	CALL XERRWD (MSG, 50, 201, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
1804	MSG = ' taken on this call before reaching TOUT '
1805	CALL XERRWD (MSG, 50, 201, 0, 1, MXSTEP, 0, 1, TN, 0.0D0)
1806	ISTATE = -1
1807	GO TO 580
1808	! EWT(I) .LE. 0.0 for some I (not at start of problem). ----------------
1809	510 EWTI = RWORK(LEWT+I-1)
1810	MSG = 'DLSODE- At T (=R1), EWT(I1) has become R2 .LE. 0.'
1811	CALL XERRWD (MSG, 50, 202, 0, 1, I, 0, 2, TN, EWTI)
1812	ISTATE = -6
1813	GO TO 580
1814	! Too much accuracy requested for machine precision. -------------------
1815	520 MSG = 'DLSODE- At T (=R1), too much accuracy requested '
1816	CALL XERRWD (MSG, 50, 203, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
1817	MSG = ' for precision of machine.. see TOLSF (=R2) '
1818	CALL XERRWD (MSG, 50, 203, 0, 0, 0, 0, 2, TN, TOLSF)
1819	RWORK(14) = TOLSF
1820	ISTATE = -2
1821	GO TO 580
1822	! KFLAG = -1. Error test failed repeatedly or with ABS(H) = HMIN. -----
1823	530 MSG = 'DLSODE- At T(=R1) and step size H(=R2), the error'
1824	CALL XERRWD (MSG, 50, 204, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
1825	MSG = ' test failed repeatedly or with ABS(H) = HMIN'
1826	CALL XERRWD (MSG, 50, 204, 0, 0, 0, 0, 2, TN, H)
1827	ISTATE = -4
1828	GO TO 560
1829	! KFLAG = -2. Convergence failed repeatedly or with ABS(H) = HMIN. ----
1830	540 MSG = 'DLSODE- At T (=R1) and step size H (=R2), the '
1831	CALL XERRWD (MSG, 50, 205, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
1832	MSG = ' corrector convergence failed repeatedly '
1833	CALL XERRWD (MSG, 50, 205, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
1834	MSG = ' or with ABS(H) = HMIN '
1835	CALL XERRWD (MSG, 30, 205, 0, 0, 0, 0, 2, TN, H)
1836	ISTATE = -5
1837	! Compute IMXER if relevant. -------------------------------------------
1838	560 BIG = 0.0D0
1839	IMXER = 1
1840	DO 570 I = 1,N
1841	SIZE = ABS(RWORK(I+LACOR-1)*RWORK(I+LEWT-1))
1842	IF (BIG .GE. SIZE) GO TO 570
1843	BIG = SIZE
1844	IMXER = I
1845	570 CONTINUE
1846	IWORK(16) = IMXER
1847	! Set Y vector, T, and optional outputs. -------------------------------
1848	580 DO 590 I = 1,N
1849	590 Y(I) = RWORK(I+LYH-1)
1850	T = TN
1851	RWORK(11) = HU
1852	RWORK(12) = H
1853	RWORK(13) = TN
1854	IWORK(11) = NST
1855	IWORK(12) = NFE
1856	IWORK(13) = NJE
1857	IWORK(14) = NQU
1858	IWORK(15) = NQ
1859	RETURN
1860	!-----------------------------------------------------------------------
1861	! Block I.
1862	! The following block handles all error returns due to illegal input
1863	! (ISTATE = -3), as detected before calling the core integrator.
1864	! First the error message routine is called. If the illegal input
1865	! is a negative ISTATE, the run is aborted (apparent infinite loop).
1866	!-----------------------------------------------------------------------
1867	601 MSG = 'DLSODE- ISTATE (=I1) illegal '
1868	CALL XERRWD (MSG, 30, 1, 0, 1, ISTATE, 0, 0, 0.0D0, 0.0D0)
1869	IF (ISTATE .LT. 0) GO TO 800
1870	GO TO 700
1871	602 MSG = 'DLSODE- ITASK (=I1) illegal '
1872	CALL XERRWD (MSG, 30, 2, 0, 1, ITASK, 0, 0, 0.0D0, 0.0D0)
1873	GO TO 700
1874	603 MSG = 'DLSODE- ISTATE .GT. 1 but DLSODE not initialized '
1875	CALL XERRWD (MSG, 50, 3, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
1876	GO TO 700
1877	604 MSG = 'DLSODE- NEQ (=I1) .LT. 1 '
1878	CALL XERRWD (MSG, 30, 4, 0, 1, NEQ, 0, 0, 0.0D0, 0.0D0)
1879	GO TO 700
1880	605 MSG = 'DLSODE- ISTATE = 3 and NEQ increased (I1 to I2) '
1881	CALL XERRWD (MSG, 50, 5, 0, 2, N, NEQ, 0, 0.0D0, 0.0D0)
1882	GO TO 700
1883	606 MSG = 'DLSODE- ITOL (=I1) illegal '
1884	CALL XERRWD (MSG, 30, 6, 0, 1, ITOL, 0, 0, 0.0D0, 0.0D0)
1885	GO TO 700
1886	607 MSG = 'DLSODE- IOPT (=I1) illegal '
1887	CALL XERRWD (MSG, 30, 7, 0, 1, IOPT, 0, 0, 0.0D0, 0.0D0)
1888	GO TO 700
1889	608 MSG = 'DLSODE- MF (=I1) illegal '
1890	CALL XERRWD (MSG, 30, 8, 0, 1, MF, 0, 0, 0.0D0, 0.0D0)
1891	GO TO 700
1892	609 MSG = 'DLSODE- ML (=I1) illegal.. .LT.0 or .GE.NEQ (=I2)'
1893	CALL XERRWD (MSG, 50, 9, 0, 2, ML, NEQ, 0, 0.0D0, 0.0D0)
1894	GO TO 700
1895	610 MSG = 'DLSODE- MU (=I1) illegal.. .LT.0 or .GE.NEQ (=I2)'
1896	CALL XERRWD (MSG, 50, 10, 0, 2, MU, NEQ, 0, 0.0D0, 0.0D0)
1897	GO TO 700
1898	611 MSG = 'DLSODE- MAXORD (=I1) .LT. 0 '
1899	CALL XERRWD (MSG, 30, 11, 0, 1, MAXORD, 0, 0, 0.0D0, 0.0D0)
1900	GO TO 700
1901	612 MSG = 'DLSODE- MXSTEP (=I1) .LT. 0 '
1902	CALL XERRWD (MSG, 30, 12, 0, 1, MXSTEP, 0, 0, 0.0D0, 0.0D0)
1903	GO TO 700
1904	613 MSG = 'DLSODE- MXHNIL (=I1) .LT. 0 '
1905	CALL XERRWD (MSG, 30, 13, 0, 1, MXHNIL, 0, 0, 0.0D0, 0.0D0)
1906	GO TO 700
1907	614 MSG = 'DLSODE- TOUT (=R1) behind T (=R2) '
1908	CALL XERRWD (MSG, 40, 14, 0, 0, 0, 0, 2, TOUT, T)
1909	MSG = ' Integration direction is given by H0 (=R1) '
1910	CALL XERRWD (MSG, 50, 14, 0, 0, 0, 0, 1, H0, 0.0D0)
1911	GO TO 700
1912	615 MSG = 'DLSODE- HMAX (=R1) .LT. 0.0 '
1913	CALL XERRWD (MSG, 30, 15, 0, 0, 0, 0, 1, HMAX, 0.0D0)
1914	GO TO 700
1915	616 MSG = 'DLSODE- HMIN (=R1) .LT. 0.0 '
1916	CALL XERRWD (MSG, 30, 16, 0, 0, 0, 0, 1, HMIN, 0.0D0)
1917	GO TO 700
1918	617 CONTINUE
1919	MSG='DLSODE- RWORK length needed, LENRW (=I1), exceeds LRW (=I2)'
1920	CALL XERRWD (MSG, 60, 17, 0, 2, LENRW, LRW, 0, 0.0D0, 0.0D0)
1921	GO TO 700
1922	618 CONTINUE
1923	MSG='DLSODE- IWORK length needed, LENIW (=I1), exceeds LIW (=I2)'
1924	CALL XERRWD (MSG, 60, 18, 0, 2, LENIW, LIW, 0, 0.0D0, 0.0D0)
1925	GO TO 700
1926	619 MSG = 'DLSODE- RelTol(I1) is R1 .LT. 0.0 '
1927	CALL XERRWD (MSG, 40, 19, 0, 1, I, 0, 1, RelTolI, 0.0D0)
1928	GO TO 700
1929	620 MSG = 'DLSODE- AbsTol(I1) is R1 .LT. 0.0 '
1930	CALL XERRWD (MSG, 40, 20, 0, 1, I, 0, 1, AbsTolI, 0.0D0)
1931	GO TO 700
1932	621 EWTI = RWORK(LEWT+I-1)
1933	MSG = 'DLSODE- EWT(I1) is R1 .LE. 0.0 '
1934	CALL XERRWD (MSG, 40, 21, 0, 1, I, 0, 1, EWTI, 0.0D0)
1935	GO TO 700
1936	622 CONTINUE
1937	MSG='DLSODE- TOUT (=R1) too close to T(=R2) to start integration'
1938	CALL XERRWD (MSG, 60, 22, 0, 0, 0, 0, 2, TOUT, T)
1939	GO TO 700
1940	623 CONTINUE
1941	MSG='DLSODE- ITASK = I1 and TOUT (=R1) behind TCUR - HU (= R2) '
1942	CALL XERRWD (MSG, 60, 23, 0, 1, ITASK, 0, 2, TOUT, TP)
1943	GO TO 700
1944	624 CONTINUE
1945	MSG='DLSODE- ITASK = 4 OR 5 and TCRIT (=R1) behind TCUR (=R2) '
1946	CALL XERRWD (MSG, 60, 24, 0, 0, 0, 0, 2, TCRIT, TN)
1947	GO TO 700
1948	625 CONTINUE
1949	MSG='DLSODE- ITASK = 4 or 5 and TCRIT (=R1) behind TOUT (=R2) '
1950	CALL XERRWD (MSG, 60, 25, 0, 0, 0, 0, 2, TCRIT, TOUT)
1951	GO TO 700
1952	626 MSG = 'DLSODE- At start of problem, too much accuracy '
1953	CALL XERRWD (MSG, 50, 26, 0, 0, 0, 0, 0, 0.0D0, 0.0D0)
1954	MSG=' requested for precision of machine.. See TOLSF (=R1) '
1955	CALL XERRWD (MSG, 60, 26, 0, 0, 0, 0, 1, TOLSF, 0.0D0)
1956	RWORK(14) = TOLSF
1957	GO TO 700
1958	627 MSG = 'DLSODE- Trouble in DINTDY. ITASK = I1, TOUT = R1'
1959	CALL XERRWD (MSG, 50, 27, 0, 1, ITASK, 0, 1, TOUT, 0.0D0)
1960	!
1961	700 ISTATE = -3
1962	RETURN
1963	!
1964	800 MSG = 'DLSODE- Run aborted.. apparent infinite loop '
1965	CALL XERRWD (MSG, 50, 303, 2, 0, 0, 0, 0, 0.0D0, 0.0D0)
1966	RETURN
1967	!----------------------- END OF SUBROUTINE DLSODE ----------------------
1968	!END SUBROUTINE DLSODE
1969	CONTAINS
1970
1971
1972	!DECK DUMACH
1973	KPP_REAL FUNCTION DUMACH ()
1974	!***BEGIN PROLOGUE DUMACH
1975	!***PURPOSE Compute the unit roundoff of the machine.
1976	!***CATEGORY R1
1977	!***TYPE KPP_REAL (RUMACH-S, DUMACH-D)
1978	!***KEYWORDS MACHINE CONSTANTS
1979	!***AUTHOR Hindmarsh, Alan C., (LLNL)
1980	!***DESCRIPTION
1981	! *Usage:
1982	! KPP_REAL A, DUMACH
1983	! A = DUMACH()
1984	!
1985	! *Function Return Values:
1986	! A : the unit roundoff of the machine.
1987	!
1988	! *Description:
1989	! The unit roundoff is defined as the smallest positive machine
1990	! number u such that 1.0 + u .ne. 1.0. This is computed by DUMACH
1991	! in a machine-independent manner.
1992	!
1993	!***REFERENCES (NONE)
1994	!***ROUTINES CALLED DUMSUM
1995	!***REVISION HISTORY (YYYYMMDD)
1996	! 19930216 DATE WRITTEN
1997	! 19930818 Added SLATEC-format prologue. (FNF)
1998	! 20030707 Added DUMSUM to force normal storage of COMP. (ACH)
1999	!***END PROLOGUE DUMACH
2000	!
2001	KPP_REAL U, COMP
2002	!***FIRST EXECUTABLE STATEMENT DUMACH
2003	U = 1.0D0
2004	10 U = U*0.5D0
2005	CALL DUMSUM(1.0D0, U, COMP)
2006	IF (COMP .NE. 1.0D0) GO TO 10
2007	DUMACH = U*2.0D0
2008	RETURN
2009	!----------------------- End of Function DUMACH ------------------------
2010	END FUNCTION DUMACH
2011
2012	SUBROUTINE DUMSUM(A,B,C)
2013	! Routine to force normal storing of A + B, for DUMACH.
2014	KPP_REAL A, B, C
2015	C = A + B
2016	RETURN
2017	END SUBROUTINE DUMSUM
2018	!DECK DCFODE
2019	SUBROUTINE DCFODE (METH, ELCO, TESCO)
2020	!***BEGIN PROLOGUE DCFODE
2021	!***SUBSIDIARY
2022	!***PURPOSE Set ODE integrator coefficients.
2023	!***TYPE KPP_REAL (SCFODE-S, DCFODE-D)
2024	!***AUTHOR Hindmarsh, Alan C., (LLNL)
2025	!***DESCRIPTION
2026	!
2027	! DCFODE is called by the integrator routine to set coefficients
2028	! needed there. The coefficients for the current method, as
2029	! given by the value of METH, are set for all orders and saved.
2030	! The maximum order assumed here is 12 if METH = 1 and 5 if METH = 2.
2031	! (A smaller value of the maximum order is also allowed.)
2032	! DCFODE is called once at the beginning of the problem,
2033	! and is not called again unless and until METH is changed.
2034	!
2035	! The ELCO array contains the basic method coefficients.
2036	! The coefficients el(i), 1 .le. i .le. nq+1, for the method of
2037	! order nq are stored in ELCO(i,nq). They are given by a genetrating
2038	! polynomial, i.e.,
2039	! l(x) = el(1) + el(2)x + ... + el(nq+1)x**nq.
2040	! For the implicit Adams methods, l(x) is given by
2041	! dl/dx = (x+1)(x+2)...*(x+nq-1)/factorial(nq-1), l(-1) = 0.
2042	! For the BDF methods, l(x) is given by
2043	! l(x) = (x+1)(x+2) ... *(x+nq)/K,
2044	! where K = factorial(nq)*(1 + 1/2 + ... + 1/nq).
2045	!
2046	! The TESCO array contains test constants used for the
2047	! local error test and the selection of step size and/or order.
2048	! At order nq, TESCO(k,nq) is used for the selection of step
2049	! size at order nq - 1 if k = 1, at order nq if k = 2, and at order
2050	! nq + 1 if k = 3.
2051	!
2052	!***SEE ALSO DLSODE
2053	!***ROUTINES CALLED (NONE)
2054	!***REVISION HISTORY (YYMMDD)
2055	! 791129 DATE WRITTEN
2056	! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
2057	! 890503 Minor cosmetic changes. (FNF)
2058	! 930809 Renamed to allow single/double precision versions. (ACH)
2059	!***END PROLOGUE DCFODE
2060	!**End
2061	INTEGER METH
2062	INTEGER I, IB, NQ, NQM1, NQP1
2063	KPP_REAL ELCO(13,12), TESCO(3,12), PC(12)
2064	KPP_REAL AGAMQ, FNQ, FNQM1, PINT, RAGQ, RQFAC, RQ1FAC, TSIGN, XPIN
2065	!
2066	!***FIRST EXECUTABLE STATEMENT DCFODE
2067	GO TO (100, 200), METH
2068	!
2069	100 ELCO(1,1) = 1.0D0
2070	ELCO(2,1) = 1.0D0
2071	TESCO(1,1) = 0.0D0
2072	TESCO(2,1) = 2.0D0
2073	TESCO(1,2) = 1.0D0
2074	TESCO(3,12) = 0.0D0
2075	PC(1) = 1.0D0
2076	RQFAC = 1.0D0
2077	DO 140 NQ = 2,12
2078	!-----------------------------------------------------------------------
2079	! The PC array will contain the coefficients of the polynomial
2080	! p(x) = (x+1)(x+2)...*(x+nq-1).
2081	! Initially, p(x) = 1.
2082	!-----------------------------------------------------------------------
2083	RQ1FAC = RQFAC
2084	RQFAC = RQFAC/NQ
2085	NQM1 = NQ - 1
2086	FNQM1 = NQM1
2087	NQP1 = NQ + 1
2088	! Form coefficients of p(x)*(x+nq-1). ----------------------------------
2089	PC(NQ) = 0.0D0
2090	DO 110 IB = 1,NQM1
2091	I = NQP1 - IB
2092	110 PC(I) = PC(I-1) + FNQM1*PC(I)
2093	PC(1) = FNQM1*PC(1)
2094	! Compute integral, -1 to 0, of p(x) and x*p(x). -----------------------
2095	PINT = PC(1)
2096	XPIN = PC(1)/2.0D0
2097	TSIGN = 1.0D0
2098	DO 120 I = 2,NQ
2099	TSIGN = -TSIGN
2100	PINT = PINT + TSIGN*PC(I)/I
2101	120 XPIN = XPIN + TSIGN*PC(I)/(I+1)
2102	! Store coefficients in ELCO and TESCO. --------------------------------
2103	ELCO(1,NQ) = PINT*RQ1FAC
2104	ELCO(2,NQ) = 1.0D0
2105	DO 130 I = 2,NQ
2106	130 ELCO(I+1,NQ) = RQ1FAC*PC(I)/I
2107	AGAMQ = RQFAC*XPIN
2108	RAGQ = 1.0D0/AGAMQ
2109	TESCO(2,NQ) = RAGQ
2110	IF (NQ .LT. 12) TESCO(1,NQP1) = RAGQ*RQFAC/NQP1
2111	TESCO(3,NQM1) = RAGQ
2112	140 CONTINUE
2113	RETURN
2114	!
2115	200 PC(1) = 1.0D0
2116	RQ1FAC = 1.0D0
2117	DO 230 NQ = 1,5
2118	!-----------------------------------------------------------------------
2119	! The PC array will contain the coefficients of the polynomial
2120	! p(x) = (x+1)(x+2)...*(x+nq).
2121	! Initially, p(x) = 1.
2122	!-----------------------------------------------------------------------
2123	FNQ = NQ
2124	NQP1 = NQ + 1
2125	! Form coefficients of p(x)*(x+nq). ------------------------------------
2126	PC(NQP1) = 0.0D0
2127	DO 210 IB = 1,NQ
2128	I = NQ + 2 - IB
2129	210 PC(I) = PC(I-1) + FNQ*PC(I)
2130	PC(1) = FNQ*PC(1)
2131	! Store coefficients in ELCO and TESCO. --------------------------------
2132	DO 220 I = 1,NQP1
2133	220 ELCO(I,NQ) = PC(I)/PC(2)
2134	ELCO(2,NQ) = 1.0D0
2135	TESCO(1,NQ) = RQ1FAC
2136	TESCO(2,NQ) = NQP1/ELCO(1,NQ)
2137	TESCO(3,NQ) = (NQ+2)/ELCO(1,NQ)
2138	RQ1FAC = RQ1FAC/FNQ
2139	230 CONTINUE
2140	RETURN
2141	!----------------------- END OF SUBROUTINE DCFODE ----------------------
2142	END SUBROUTINE DCFODE
2143	!DECK DINTDY
2144	SUBROUTINE DINTDY (T, K, YH, NYH, DKY, IFLAG)
2145	!***BEGIN PROLOGUE DINTDY
2146	!***SUBSIDIARY
2147	!***PURPOSE Interpolate solution derivatives.
2148	!***TYPE KPP_REAL (SINTDY-S, DINTDY-D)
2149	!***AUTHOR Hindmarsh, Alan C., (LLNL)
2150	!***DESCRIPTION
2151	!
2152	! DINTDY computes interpolated values of the K-th derivative of the
2153	! dependent variable vector y, and stores it in DKY. This routine
2154	! is called within the package with K = 0 and T = TOUT, but may
2155	! also be called by the user for any K up to the current order.
2156	! (See detailed instructions in the usage documentation.)
2157	!
2158	! The computed values in DKY are gotten by interpolation using the
2159	! Nordsieck history array YH. This array corresponds uniquely to a
2160	! vector-valued polynomial of degree NQCUR or less, and DKY is set
2161	! to the K-th derivative of this polynomial at T.
2162	! The formula for DKY is:
2163	! q
2164	! DKY(i) = sum c(j,K) * (T - tn)*(j-K) h*(-j) YH(i,j+1)
2165	! j=K
2166	! where c(j,K) = j(j-1)...*(j-K+1), q = NQCUR, tn = TCUR, h = HCUR.
2167	! The quantities nq = NQCUR, l = nq+1, N = NEQ, tn, and h are
2168	! communicated by COMMON. The above sum is done in reverse order.
2169	! IFLAG is returned negative if either K or T is out of bounds.
2170	!
2171	!***SEE ALSO DLSODE
2172	!***ROUTINES CALLED XERRWD
2173	!***COMMON BLOCKS DLS001
2174	!***REVISION HISTORY (YYMMDD)
2175	! 791129 DATE WRITTEN
2176	! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
2177	! 890503 Minor cosmetic changes. (FNF)
2178	! 930809 Renamed to allow single/double precision versions. (ACH)
2179	! 010418 Reduced size of Common block /DLS001/. (ACH)
2180	! 031105 Restored 'own' variables to Common block /DLS001/, to
2181	! enable interrupt/restart feature. (ACH)
2182	! 050427 Corrected roundoff decrement in TP. (ACH)
2183	!***END PROLOGUE DINTDY
2184	!**End
2185	INTEGER K, NYH, IFLAG
2186	KPP_REAL T, YH(NYH,), DKY()
2187	INTEGER IOWND, IOWNS, &
2188	ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, &
2189	LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, &
2190	MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
2191	KPP_REAL ROWNS, &
2192	CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
2193	COMMON /DLS001/ ROWNS(209), &
2194	CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, &
2195	IOWND(6), IOWNS(6), &
2196	ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, &
2197	LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, &
2198	MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
2199	INTEGER I, IC, J, JB, JB2, JJ, JJ1, JP1
2200	KPP_REAL C, R, S, TP
2201	CHARACTER*80 MSG
2202	!
2203	!***FIRST EXECUTABLE STATEMENT DINTDY
2204	IFLAG = 0
2205	IF (K .LT. 0 .OR. K .GT. NQ) GO TO 80
2206	TP = TN - HU - 100.0D0UROUNDSIGN(ABS(TN) + ABS(HU), HU)
2207	IF ((T-TP)*(T-TN) .GT. 0.0D0) GO TO 90
2208	!
2209	S = (T - TN)/H
2210	IC = 1
2211	IF (K .EQ. 0) GO TO 15
2212	JJ1 = L - K
2213	DO 10 JJ = JJ1,NQ
2214	10 IC = IC*JJ
2215	15 C = IC
2216	DO 20 I = 1,N
2217	20 DKY(I) = C*YH(I,L)
2218	IF (K .EQ. NQ) GO TO 55
2219	JB2 = NQ - K
2220	DO 50 JB = 1,JB2
2221	J = NQ - JB
2222	JP1 = J + 1
2223	IC = 1
2224	IF (K .EQ. 0) GO TO 35
2225	JJ1 = JP1 - K
2226	DO 30 JJ = JJ1,J
2227	30 IC = IC*JJ
2228	35 C = IC
2229	DO 40 I = 1,N
2230	40 DKY(I) = CYH(I,JP1) + SDKY(I)
2231	50 CONTINUE
2232	IF (K .EQ. 0) RETURN
2233	55 R = H**(-K)
2234	DO 60 I = 1,N
2235	60 DKY(I) = R*DKY(I)
2236	RETURN
2237	!
2238	80 MSG = 'DINTDY- K (=I1) illegal '
2239	CALL XERRWD (MSG, 30, 51, 0, 1, K, 0, 0, 0.0D0, 0.0D0)
2240	IFLAG = -1
2241	RETURN
2242	90 MSG = 'DINTDY- T (=R1) illegal '
2243	CALL XERRWD (MSG, 30, 52, 0, 0, 0, 0, 1, T, 0.0D0)
2244	MSG=' T not in interval TCUR - HU (= R1) to TCUR (=R2) '
2245	CALL XERRWD (MSG, 60, 52, 0, 0, 0, 0, 2, TP, TN)
2246	IFLAG = -2
2247	RETURN
2248	!----------------------- END OF SUBROUTINE DINTDY ----------------------
2249	END SUBROUTINE DINTDY
2250	!DECK DPREPJ
2251	SUBROUTINE DPREPJ (NEQ, Y, YH, NYH, EWT, FTEM, SAVF, WM, IWM, F, JAC)
2252	!***BEGIN PROLOGUE DPREPJ
2253	!***SUBSIDIARY
2254	!***PURPOSE Compute and process Newton iteration matrix.
2255	!***TYPE KPP_REAL (SPREPJ-S, DPREPJ-D)
2256	!***AUTHOR Hindmarsh, Alan C., (LLNL)
2257	!***DESCRIPTION
2258	!
2259	! DPREPJ is called by DSTODE to compute and process the matrix
2260	! P = I - hel(1)J , where J is an approximation to the Jacobian.
2261	! Here J is computed by the user-supplied routine JAC if
2262	! MITER = 1 or 4, or by finite differencing if MITER = 2, 3, or 5.
2263	! If MITER = 3, a diagonal approximation to J is used.
2264	! J is stored in WM and replaced by P. If MITER .ne. 3, P is then
2265	! subjected to LU decomposition in preparation for later solution
2266	! of linear systems with P as coefficient matrix. This is done
2267	! by DGEFA if MITER = 1 or 2, and by DGBFA if MITER = 4 or 5.
2268	!
2269	! In addition to variables described in DSTODE and DLSODE prologues,
2270	! communication with DPREPJ uses the following:
2271	! Y = array containing predicted values on entry.
2272	! FTEM = work array of length N (ACOR in DSTODE).
2273	! SAVF = array containing f evaluated at predicted y.
2274	! WM = real work space for matrices. On output it contains the
2275	! inverse diagonal matrix if MITER = 3 and the LU decomposition
2276	! of P if MITER is 1, 2 , 4, or 5.
2277	! Storage of matrix elements starts at WM(3).
2278	! WM also contains the following matrix-related data:
2279	! WM(1) = SQRT(UROUND), used in numerical Jacobian increments.
2280	! WM(2) = H*EL0, saved for later use if MITER = 3.
2281	! IWM = integer work space containing pivot information, starting at
2282	! IWM(21), if MITER is 1, 2, 4, or 5. IWM also contains band
2283	! parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5.
2284	! EL0 = EL(1) (input).
2285	! IERPJ = output error flag, = 0 if no trouble, .gt. 0 if
2286	! P matrix found to be singular.
2287	! JCUR = output flag = 1 to indicate that the Jacobian matrix
2288	! (or approximation) is now current.
2289	! This routine also uses the COMMON variables EL0, H, TN, UROUND,
2290	! MITER, N, NFE, and NJE.
2291	!
2292	!***SEE ALSO DLSODE
2293	!***ROUTINES CALLED DGBFA, DGEFA, DVNORM
2294	!***COMMON BLOCKS DLS001
2295	!***REVISION HISTORY (YYMMDD)
2296	! 791129 DATE WRITTEN
2297	! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
2298	! 890504 Minor cosmetic changes. (FNF)
2299	! 930809 Renamed to allow single/double precision versions. (ACH)
2300	! 010418 Reduced size of Common block /DLS001/. (ACH)
2301	! 031105 Restored 'own' variables to Common block /DLS001/, to
2302	! enable interrupt/restart feature. (ACH)
2303	!***END PROLOGUE DPREPJ
2304	!**End
2305	EXTERNAL F, JAC
2306	INTEGER NEQ, NYH, IWM(*)
2307	KPP_REAL Y(), YH(NYH,), EWT(), FTEM(), SAVF(), WM()
2308	INTEGER IOWND, IOWNS, IER, &
2309	ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, &
2310	LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, &
2311	MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
2312	KPP_REAL ROWNS, &
2313	CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
2314	COMMON /DLS001/ ROWNS(209), &
2315	CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, &
2316	IOWND(6), IOWNS(6), &
2317	ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, &
2318	LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, &
2319	MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
2320	INTEGER I, I1, I2, ISING, II, J, J1, JJ, LENP, &
2321	MBA, MBAND, MEB1, MEBAND, ML, ML3, MU, NP1
2322	KPP_REAL CON, DI, FAC, HL0, R, R0, SRUR, YI, YJ, YJJ
2323	!KPP_REAL DVNORM
2324	!
2325	!***FIRST EXECUTABLE STATEMENT DPREPJ
2326	NJE = NJE + 1
2327	IERPJ = 0
2328	JCUR = 1
2329	HL0 = H*EL0
2330	CON = -HL0
2331
2332	#ifdef FULL_ALGEBRA
2333	LENP = N*N
2334	DO i = 1,LENP
2335	WM(i+2) = 0.0D0
2336	END DO
2337	CALL JAC_CHEM (NEQ, TN, Y, WM(3))
2338	DO I = 1,LENP
2339	WM(I+2) = WM(I+2)*CON
2340	END DO
2341	! Add identity matrix
2342	J = 3
2343	NP1 = N + 1
2344	DO I = 1,N
2345	WM(J) = WM(J) + 1.0D0
2346	J = J + NP1
2347	END DO
2348	! Do LU decomposition on P
2349	CALL DGETRF(N,N,WM(3),N,IWM(21),ISING)
2350	#else
2351	CALL JAC_CHEM (NEQ, TN, Y, WM(3))
2352	DO i = 1,LU_NONZERO
2353	WM(i+2) = WM(i+2)*CON
2354	END DO
2355	! Add identity matrix
2356	DO i = 1,N
2357	j = 2+LU_DIAG(i)
2358	WM(j) = WM(j) + 1.0D0
2359	END DO
2360	! Do LU decomposition on P
2361	CALL KppDecomp(WM(3),IER)
2362	#endif
2363	IF (IER .NE. 0) IERPJ = 1
2364	RETURN
2365	!----------------------- END OF SUBROUTINE DPREPJ ----------------------
2366	END SUBROUTINE DPREPJ
2367	!DECK DSOLSY
2368	SUBROUTINE DSOLSY (WM, IWM, X, TEM)
2369	!***BEGIN PROLOGUE DSOLSY
2370	!***SUBSIDIARY
2371	!***PURPOSE ODEPACK linear system solver.
2372	!***TYPE KPP_REAL (SSOLSY-S, DSOLSY-D)
2373	!***AUTHOR Hindmarsh, Alan C., (LLNL)
2374	!***DESCRIPTION
2375	!
2376	! This routine manages the solution of the linear system arising from
2377	! a chord iteration. It is called if MITER .ne. 0.
2378	! If MITER is 1 or 2, it calls DGESL to accomplish this.
2379	! If MITER = 3 it updates the coefficient h*EL0 in the diagonal
2380	! matrix, and then computes the solution.
2381	! If MITER is 4 or 5, it calls DGBSL.
2382	! Communication with DSOLSY uses the following variables:
2383	! WM = real work space containing the inverse diagonal matrix if
2384	! MITER = 3 and the LU decomposition of the matrix otherwise.
2385	! Storage of matrix elements starts at WM(3).
2386	! WM also contains the following matrix-related data:
2387	! WM(1) = SQRT(UROUND) (not used here),
2388	! WM(2) = HL0, the previous value of h*EL0, used if MITER = 3.
2389	! IWM = integer work space containing pivot information, starting at
2390	! IWM(21), if MITER is 1, 2, 4, or 5. IWM also contains band
2391	! parameters ML = IWM(1) and MU = IWM(2) if MITER is 4 or 5.
2392	! X = the right-hand side vector on input, and the solution vector
2393	! on output, of length N.
2394	! TEM = vector of work space of length N, not used in this version.
2395	! IERSL = output flag (in COMMON). IERSL = 0 if no trouble occurred.
2396	! IERSL = 1 if a singular matrix arose with MITER = 3.
2397	! This routine also uses the COMMON variables EL0, H, MITER, and N.
2398	!
2399	!***SEE ALSO DLSODE
2400	!***ROUTINES CALLED DGBSL, DGESL
2401	!***COMMON BLOCKS DLS001
2402	!***REVISION HISTORY (YYMMDD)
2403	! 791129 DATE WRITTEN
2404	! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
2405	! 890503 Minor cosmetic changes. (FNF)
2406	! 930809 Renamed to allow single/double precision versions. (ACH)
2407	! 010418 Reduced size of Common block /DLS001/. (ACH)
2408	! 031105 Restored 'own' variables to Common block /DLS001/, to
2409	! enable interrupt/restart feature. (ACH)
2410	!***END PROLOGUE DSOLSY
2411	!**End
2412	INTEGER IWM(*)
2413	KPP_REAL WM(), X(), TEM(*)
2414	INTEGER IOWND, IOWNS, &
2415	ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, &
2416	LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, &
2417	MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
2418	KPP_REAL ROWNS, &
2419	CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
2420	COMMON /DLS001/ ROWNS(209), &
2421	CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, &
2422	IOWND(6), IOWNS(6), &
2423	ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, &
2424	LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, &
2425	MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
2426	INTEGER I, MEBAND, ML, MU
2427	KPP_REAL DI, HL0, PHL0, R
2428	#ifdef FULL_ALGEBRA
2429	INTEGER ISING
2430	#endif
2431	!
2432	!***FIRST EXECUTABLE STATEMENT DSOLSY
2433	IERSL = 0
2434	#ifdef FULL_ALGEBRA
2435	CALL DGETRS ('N',N,1,WM(3),N,IWM(21),X,N,ISING)
2436	#else
2437	CALL KppSolve(WM(3),X)
2438	#endif
2439	RETURN
2440	!----------------------- END OF SUBROUTINE DSOLSY ----------------------
2441	END SUBROUTINE DSOLSY
2442	!DECK DSRCOM
2443	SUBROUTINE DSRCOM (RSAV, ISAV, JOB)
2444	!***BEGIN PROLOGUE DSRCOM
2445	!***SUBSIDIARY
2446	!***PURPOSE Save/restore ODEPACK COMMON blocks.
2447	!***TYPE KPP_REAL (SSRCOM-S, DSRCOM-D)
2448	!***AUTHOR Hindmarsh, Alan C., (LLNL)
2449	!***DESCRIPTION
2450	!
2451	! This routine saves or restores (depending on JOB) the contents of
2452	! the COMMON block DLS001, which is used internally
2453	! by one or more ODEPACK solvers.
2454	!
2455	! RSAV = real array of length 218 or more.
2456	! ISAV = integer array of length 37 or more.
2457	! JOB = flag indicating to save or restore the COMMON blocks:
2458	! JOB = 1 if COMMON is to be saved (written to RSAV/ISAV)
2459	! JOB = 2 if COMMON is to be restored (read from RSAV/ISAV)
2460	! A call with JOB = 2 presumes a prior call with JOB = 1.
2461	!
2462	!***SEE ALSO DLSODE
2463	!***ROUTINES CALLED (NONE)
2464	!***COMMON BLOCKS DLS001
2465	!***REVISION HISTORY (YYMMDD)
2466	! 791129 DATE WRITTEN
2467	! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
2468	! 890503 Minor cosmetic changes. (FNF)
2469	! 921116 Deleted treatment of block /EH0001/. (ACH)
2470	! 930801 Reduced Common block length by 2. (ACH)
2471	! 930809 Renamed to allow single/double precision versions. (ACH)
2472	! 010418 Reduced Common block length by 209+12. (ACH)
2473	! 031105 Restored 'own' variables to Common block /DLS001/, to
2474	! enable interrupt/restart feature. (ACH)
2475	! 031112 Added SAVE statement for data-loaded constants.
2476	!***END PROLOGUE DSRCOM
2477	!**End
2478	INTEGER ISAV(*), JOB
2479	INTEGER ILS
2480	INTEGER I, LENILS, LENRLS
2481	KPP_REAL RSAV(*), RLS
2482	SAVE LENRLS, LENILS
2483	COMMON /DLS001/ RLS(218), ILS(37)
2484	DATA LENRLS/218/, LENILS/37/
2485	!
2486	!***FIRST EXECUTABLE STATEMENT DSRCOM
2487	IF (JOB .EQ. 2) GO TO 100
2488	!
2489	DO 10 I = 1,LENRLS
2490	10 RSAV(I) = RLS(I)
2491	DO 20 I = 1,LENILS
2492	20 ISAV(I) = ILS(I)
2493	RETURN
2494	!
2495	100 CONTINUE
2496	DO 110 I = 1,LENRLS
2497	110 RLS(I) = RSAV(I)
2498	DO 120 I = 1,LENILS
2499	120 ILS(I) = ISAV(I)
2500	RETURN
2501	!----------------------- END OF SUBROUTINE DSRCOM ----------------------
2502	END SUBROUTINE DSRCOM
2503	!DECK DSTODE
2504	SUBROUTINE DSTODE (NEQ, Y, YH, NYH, YH1, EWT, SAVF, ACOR, &
2505	WM, IWM, F, JAC)
2506	!WM, IWM, F, JAC, PJAC, SLVS)
2507	!***BEGIN PROLOGUE DSTODE
2508	!***SUBSIDIARY
2509	!***PURPOSE Performs one step of an ODEPACK integration.
2510	!***TYPE KPP_REAL (SSTODE-S, DSTODE-D)
2511	!***AUTHOR Hindmarsh, Alan C., (LLNL)
2512	!***DESCRIPTION
2513	!
2514	! DSTODE performs one step of the integration of an initial value
2515	! problem for a system of ordinary differential equations.
2516	! Note: DSTODE is independent of the value of the iteration method
2517	! indicator MITER, when this is .ne. 0, and hence is independent
2518	! of the type of chord method used, or the Jacobian structure.
2519	! Communication with DSTODE is done with the following variables:
2520	!
2521	! NEQ = integer array containing problem size in NEQ, and
2522	! passed as the NEQ argument in all calls to F and JAC.
2523	! Y = an array of length .ge. N used as the Y argument in
2524	! all calls to F and JAC.
2525	! YH = an NYH by LMAX array containing the dependent variables
2526	! and their approximate scaled derivatives, where
2527	! LMAX = MAXORD + 1. YH(i,j+1) contains the approximate
2528	! j-th derivative of y(i), scaled by h**j/factorial(j)
2529	! (j = 0,1,...,NQ). on entry for the first step, the first
2530	! two columns of YH must be set from the initial values.
2531	! NYH = a constant integer .ge. N, the first dimension of YH.
2532	! YH1 = a one-dimensional array occupying the same space as YH.
2533	! EWT = an array of length N containing multiplicative weights
2534	! for local error measurements. Local errors in Y(i) are
2535	! compared to 1.0/EWT(i) in various error tests.
2536	! SAVF = an array of working storage, of length N.
2537	! Also used for input of YH(*,MAXORD+2) when JSTART = -1
2538	! and MAXORD .lt. the current order NQ.
2539	! ACOR = a work array of length N, used for the accumulated
2540	! corrections. On a successful return, ACOR(i) contains
2541	! the estimated one-step local error in Y(i).
2542	! WM,IWM = real and integer work arrays associated with matrix
2543	! operations in chord iteration (MITER .ne. 0).
2544	! PJAC = name of routine to evaluate and preprocess Jacobian matrix
2545	! and P = I - hel0JAC, if a chord method is being used.
2546	! SLVS = name of routine to solve linear system in chord iteration.
2547	! CCMAX = maximum relative change in h*el0 before PJAC is called.
2548	! H = the step size to be attempted on the next step.
2549	! H is altered by the error control algorithm during the
2550	! problem. H can be either positive or negative, but its
2551	! sign must remain constant throughout the problem.
2552	! HMIN = the minimum absolute value of the step size h to be used.
2553	! HMXI = inverse of the maximum absolute value of h to be used.
2554	! HMXI = 0.0 is allowed and corresponds to an infinite hmax.
2555	! HMIN and HMXI may be changed at any time, but will not
2556	! take effect until the next change of h is considered.
2557	! TN = the independent variable. TN is updated on each step taken.
2558	! JSTART = an integer used for input only, with the following
2559	! values and meanings:
2560	! 0 perform the first step.
2561	! .gt.0 take a new step continuing from the last.
2562	! -1 take the next step with a new value of H, MAXORD,
2563	! N, METH, MITER, and/or matrix parameters.
2564	! -2 take the next step with a new value of H,
2565	! but with other inputs unchanged.
2566	! On return, JSTART is set to 1 to facilitate continuation.
2567	! KFLAG = a completion code with the following meanings:
2568	! 0 the step was succesful.
2569	! -1 the requested error could not be achieved.
2570	! -2 corrector convergence could not be achieved.
2571	! -3 fatal error in PJAC or SLVS.
2572	! A return with KFLAG = -1 or -2 means either
2573	! abs(H) = HMIN or 10 consecutive failures occurred.
2574	! On a return with KFLAG negative, the values of TN and
2575	! the YH array are as of the beginning of the last
2576	! step, and H is the last step size attempted.
2577	! MAXORD = the maximum order of integration method to be allowed.
2578	! MAXCOR = the maximum number of corrector iterations allowed.
2579	! MSBP = maximum number of steps between PJAC calls (MITER .gt. 0).
2580	! MXNCF = maximum number of convergence failures allowed.
2581	! METH/MITER = the method flags. See description in driver.
2582	! N = the number of first-order differential equations.
2583	! The values of CCMAX, H, HMIN, HMXI, TN, JSTART, KFLAG, MAXORD,
2584	! MAXCOR, MSBP, MXNCF, METH, MITER, and N are communicated via COMMON.
2585	!
2586	!***SEE ALSO DLSODE
2587	!***ROUTINES CALLED DCFODE, DVNORM
2588	!***COMMON BLOCKS DLS001
2589	!***REVISION HISTORY (YYMMDD)
2590	! 791129 DATE WRITTEN
2591	! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
2592	! 890503 Minor cosmetic changes. (FNF)
2593	! 930809 Renamed to allow single/double precision versions. (ACH)
2594	! 010418 Reduced size of Common block /DLS001/. (ACH)
2595	! 031105 Restored 'own' variables to Common block /DLS001/, to
2596	! enable interrupt/restart feature. (ACH)
2597	!***END PROLOGUE DSTODE
2598	!**End
2599	EXTERNAL F, JAC !, PJAC, SLVS
2600	INTEGER NEQ, NYH, IWM(*)
2601	KPP_REAL Y(), YH(NYH,), YH1(), EWT(), SAVF(*), &
2602	ACOR(), WM()
2603	INTEGER IOWND, IALTH, IPUP, LMAX, MEO, NQNYH, NSLP, &
2604	ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, &
2605	LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, &
2606	MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
2607	INTEGER I, I1, IREDO, IRET, J, JB, M, NCF, NEWQ
2608	KPP_REAL CONIT, CRATE, EL, ELCO, HOLD, RMAX, TESCO, &
2609	CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND
2610	KPP_REAL DCON, DDN, DEL, DELP, DSM, DUP, EXDN, EXSM, &
2611	EXUP,R, RH, RHDN, RHSM, RHUP, TOLD
2612	!KPP_REAL DVNORM
2613	COMMON /DLS001/ CONIT, CRATE, EL(13), ELCO(13,12), &
2614	HOLD, RMAX, TESCO(3,12), &
2615	CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND, &
2616	IOWND(6), IALTH, IPUP, LMAX, MEO, NQNYH, NSLP, &
2617	ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, &
2618	LYH, LEWT, LACOR, LSAVF, LWM, LIWM, METH, MITER, &
2619	MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
2620	!
2621	!***FIRST EXECUTABLE STATEMENT DSTODE
2622	KFLAG = 0
2623	TOLD = TN
2624	NCF = 0
2625	IERPJ = 0
2626	IERSL = 0
2627	JCUR = 0
2628	ICF = 0
2629	DELP = 0.0D0
2630	IF (JSTART .GT. 0) GO TO 200
2631	IF (JSTART .EQ. -1) GO TO 100
2632	IF (JSTART .EQ. -2) GO TO 160
2633	!-----------------------------------------------------------------------
2634	! On the first call, the order is set to 1, and other variables are
2635	! initialized. RMAX is the maximum ratio by which H can be increased
2636	! in a single step. It is initially 1.E4 to compensate for the small
2637	! initial H, but then is normally equal to 10. If a failure
2638	! occurs (in corrector convergence or error test), RMAX is set to 2
2639	! for the next increase.
2640	!-----------------------------------------------------------------------
2641	LMAX = MAXORD + 1
2642	NQ = 1
2643	L = 2
2644	IALTH = 2
2645	RMAX = 10000.0D0
2646	RC = 0.0D0
2647	EL0 = 1.0D0
2648	CRATE = 0.7D0
2649	HOLD = H
2650	MEO = METH
2651	NSLP = 0
2652	IPUP = MITER
2653	IRET = 3
2654	GO TO 140
2655	!-----------------------------------------------------------------------
2656	! The following block handles preliminaries needed when JSTART = -1.
2657	! IPUP is set to MITER to force a matrix update.
2658	! If an order increase is about to be considered (IALTH = 1),
2659	! IALTH is reset to 2 to postpone consideration one more step.
2660	! If the caller has changed METH, DCFODE is called to reset
2661	! the coefficients of the method.
2662	! If the caller has changed MAXORD to a value less than the current
2663	! order NQ, NQ is reduced to MAXORD, and a new H chosen accordingly.
2664	! If H is to be changed, YH must be rescaled.
2665	! If H or METH is being changed, IALTH is reset to L = NQ + 1
2666	! to prevent further changes in H for that many steps.
2667	!-----------------------------------------------------------------------
2668	100 IPUP = MITER
2669	LMAX = MAXORD + 1
2670	IF (IALTH .EQ. 1) IALTH = 2
2671	IF (METH .EQ. MEO) GO TO 110
2672	CALL DCFODE (METH, ELCO, TESCO)
2673	MEO = METH
2674	IF (NQ .GT. MAXORD) GO TO 120
2675	IALTH = L
2676	IRET = 1
2677	GO TO 150
2678	110 IF (NQ .LE. MAXORD) GO TO 160
2679	120 NQ = MAXORD
2680	L = LMAX
2681	DO 125 I = 1,L
2682	125 EL(I) = ELCO(I,NQ)
2683	NQNYH = NQ*NYH
2684	RC = RC*EL(1)/EL0
2685	EL0 = EL(1)
2686	CONIT = 0.5D0/(NQ+2)
2687	DDN = DVNORM (N, SAVF, EWT)/TESCO(1,L)
2688	EXDN = 1.0D0/L
2689	RHDN = 1.0D0/(1.3D0DDN*EXDN + 0.0000013D0)
2690	RH = MIN(RHDN,1.0D0)
2691	IREDO = 3
2692	IF (H .EQ. HOLD) GO TO 170
2693	RH = MIN(RH,ABS(H/HOLD))
2694	H = HOLD
2695	GO TO 175
2696	!-----------------------------------------------------------------------
2697	! DCFODE is called to get all the integration coefficients for the
2698	! current METH. Then the EL vector and related constants are reset
2699	! whenever the order NQ is changed, or at the start of the problem.
2700	!-----------------------------------------------------------------------
2701	140 CALL DCFODE (METH, ELCO, TESCO)
2702	150 DO 155 I = 1,L
2703	155 EL(I) = ELCO(I,NQ)
2704	NQNYH = NQ*NYH
2705	RC = RC*EL(1)/EL0
2706	EL0 = EL(1)
2707	CONIT = 0.5D0/(NQ+2)
2708	GO TO (160, 170, 200), IRET
2709	!-----------------------------------------------------------------------
2710	! If H is being changed, the H ratio RH is checked against
2711	! RMAX, HMIN, and HMXI, and the YH array rescaled. IALTH is set to
2712	! L = NQ + 1 to prevent a change of H for that many steps, unless
2713	! forced by a convergence or error test failure.
2714	!-----------------------------------------------------------------------
2715	160 IF (H .EQ. HOLD) GO TO 200
2716	RH = H/HOLD
2717	H = HOLD
2718	IREDO = 3
2719	GO TO 175
2720	170 RH = MAX(RH,HMIN/ABS(H))
2721	175 RH = MIN(RH,RMAX)
2722	RH = RH/MAX(1.0D0,ABS(H)HMXIRH)
2723	R = 1.0D0
2724	DO 180 J = 2,L
2725	R = R*RH
2726	DO 180 I = 1,N
2727	180 YH(I,J) = YH(I,J)*R
2728	H = H*RH
2729	RC = RC*RH
2730	IALTH = L
2731	IF (IREDO .EQ. 0) GO TO 690
2732	!-----------------------------------------------------------------------
2733	! This section computes the predicted values by effectively
2734	! multiplying the YH array by the Pascal Triangle matrix.
2735	! RC is the ratio of new to old values of the coefficient H*EL(1).
2736	! When RC differs from 1 by more than CCMAX, IPUP is set to MITER
2737	! to force PJAC to be called, if a Jacobian is involved.
2738	! In any case, PJAC is called at least every MSBP steps.
2739	!-----------------------------------------------------------------------
2740	200 IF (ABS(RC-1.0D0) .GT. CCMAX) IPUP = MITER
2741	IF (NST .GE. NSLP+MSBP) IPUP = MITER
2742	TN = TN + H
2743	I1 = NQNYH + 1
2744	DO 215 JB = 1,NQ
2745	I1 = I1 - NYH
2746	!dir$ ivdep
2747	DO 210 I = I1,NQNYH
2748	210 YH1(I) = YH1(I) + YH1(I+NYH)
2749	215 CONTINUE
2750	!-----------------------------------------------------------------------
2751	! Up to MAXCOR corrector iterations are taken. A convergence test is
2752	! made on the R.M.S. norm of each correction, weighted by the error
2753	! weight vector EWT. The sum of the corrections is accumulated in the
2754	! vector ACOR(i). The YH array is not altered in the corrector loop.
2755	!-----------------------------------------------------------------------
2756	220 M = 0
2757	DO 230 I = 1,N
2758	230 Y(I) = YH(I,1)
2759	CALL F (NEQ, TN, Y, SAVF)
2760	NFE = NFE + 1
2761	IF (IPUP .LE. 0) GO TO 250
2762	!-----------------------------------------------------------------------
2763	! If indicated, the matrix P = I - hel(1)J is reevaluated and
2764	! preprocessed before starting the corrector iteration. IPUP is set
2765	! to 0 as an indicator that this has been done.
2766	!-----------------------------------------------------------------------
2767	!CALL PJAC (NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC)
2768	CALL DPREPJ(NEQ, Y, YH, NYH, EWT, ACOR, SAVF, WM, IWM, F, JAC)
2769	IPUP = 0
2770	RC = 1.0D0
2771	NSLP = NST
2772	CRATE = 0.7D0
2773	IF (IERPJ .NE. 0) GO TO 430
2774	250 DO 260 I = 1,N
2775	260 ACOR(I) = 0.0D0
2776	270 IF (MITER .NE. 0) GO TO 350
2777	!-----------------------------------------------------------------------
2778	! In the case of functional iteration, update Y directly from
2779	! the result of the last function evaluation.
2780	!-----------------------------------------------------------------------
2781	DO 290 I = 1,N
2782	SAVF(I) = H*SAVF(I) - YH(I,2)
2783	290 Y(I) = SAVF(I) - ACOR(I)
2784	DEL = DVNORM (N, Y, EWT)
2785	DO 300 I = 1,N
2786	Y(I) = YH(I,1) + EL(1)*SAVF(I)
2787	300 ACOR(I) = SAVF(I)
2788	GO TO 400
2789	!-----------------------------------------------------------------------
2790	! In the case of the chord method, compute the corrector error,
2791	! and solve the linear system with that as right-hand side and
2792	! P as coefficient matrix.
2793	!-----------------------------------------------------------------------
2794	350 DO 360 I = 1,N
2795	360 Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I))
2796	!CALL SLVS (WM, IWM, Y, SAVF)
2797	CALL DSOLSY(WM, IWM, Y, SAVF)
2798	IF (IERSL .LT. 0) GO TO 430
2799	IF (IERSL .GT. 0) GO TO 410
2800	DEL = DVNORM (N, Y, EWT)
2801	DO 380 I = 1,N
2802	ACOR(I) = ACOR(I) + Y(I)
2803	380 Y(I) = YH(I,1) + EL(1)*ACOR(I)
2804	!-----------------------------------------------------------------------
2805	! Test for convergence. If M.gt.0, an estimate of the convergence
2806	! rate constant is stored in CRATE, and this is used in the test.
2807	!-----------------------------------------------------------------------
2808	400 IF (M .NE. 0) CRATE = MAX(0.2D0*CRATE,DEL/DELP)
2809	DCON = DELMIN(1.0D0,1.5D0CRATE)/(TESCO(2,NQ)*CONIT)
2810	IF (DCON .LE. 1.0D0) GO TO 450
2811	M = M + 1
2812	IF (M .EQ. MAXCOR) GO TO 410
2813	IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP) GO TO 410
2814	DELP = DEL
2815	CALL F (NEQ, TN, Y, SAVF)
2816	NFE = NFE + 1
2817	GO TO 270
2818	!-----------------------------------------------------------------------
2819	! The corrector iteration failed to converge.
2820	! If MITER .ne. 0 and the Jacobian is out of date, PJAC is called for
2821	! the next try. Otherwise the YH array is retracted to its values
2822	! before prediction, and H is reduced, if possible. If H cannot be
2823	! reduced or MXNCF failures have occurred, exit with KFLAG = -2.
2824	!-----------------------------------------------------------------------
2825	410 IF (MITER .EQ. 0 .OR. JCUR .EQ. 1) GO TO 430
2826	ICF = 1
2827	IPUP = MITER
2828	GO TO 220
2829	430 ICF = 2
2830	NCF = NCF + 1
2831	RMAX = 2.0D0
2832	TN = TOLD
2833	I1 = NQNYH + 1
2834	DO 445 JB = 1,NQ
2835	I1 = I1 - NYH
2836	!dir$ ivdep
2837	DO 440 I = I1,NQNYH
2838	440 YH1(I) = YH1(I) - YH1(I+NYH)
2839	445 CONTINUE
2840	IF (IERPJ .LT. 0 .OR. IERSL .LT. 0) GO TO 680
2841	IF (ABS(H) .LE. HMIN*1.00001D0) GO TO 670
2842	IF (NCF .EQ. MXNCF) GO TO 670
2843	RH = 0.25D0
2844	IPUP = MITER
2845	IREDO = 1
2846	GO TO 170
2847	!-----------------------------------------------------------------------
2848	! The corrector has converged. JCUR is set to 0
2849	! to signal that the Jacobian involved may need updating later.
2850	! The local error test is made and control passes to statement 500
2851	! if it fails.
2852	!-----------------------------------------------------------------------
2853	450 JCUR = 0
2854	IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ)
2855	IF (M .GT. 0) DSM = DVNORM (N, ACOR, EWT)/TESCO(2,NQ)
2856	IF (DSM .GT. 1.0D0) GO TO 500
2857	!-----------------------------------------------------------------------
2858	! After a successful step, update the YH array.
2859	! Consider changing H if IALTH = 1. Otherwise decrease IALTH by 1.
2860	! If IALTH is then 1 and NQ .lt. MAXORD, then ACOR is saved for
2861	! use in a possible order increase on the next step.
2862	! If a change in H is considered, an increase or decrease in order
2863	! by one is considered also. A change in H is made only if it is by a
2864	! factor of at least 1.1. If not, IALTH is set to 3 to prevent
2865	! testing for that many steps.
2866	!-----------------------------------------------------------------------
2867	KFLAG = 0
2868	IREDO = 0
2869	NST = NST + 1
2870	HU = H
2871	NQU = NQ
2872	DO 470 J = 1,L
2873	DO 470 I = 1,N
2874	470 YH(I,J) = YH(I,J) + EL(J)*ACOR(I)
2875	IALTH = IALTH - 1
2876	IF (IALTH .EQ. 0) GO TO 520
2877	IF (IALTH .GT. 1) GO TO 700
2878	IF (L .EQ. LMAX) GO TO 700
2879	DO 490 I = 1,N
2880	490 YH(I,LMAX) = ACOR(I)
2881	GO TO 700
2882	!-----------------------------------------------------------------------
2883	! The error test failed. KFLAG keeps track of multiple failures.
2884	! Restore TN and the YH array to their previous values, and prepare
2885	! to try the step again. Compute the optimum step size for this or
2886	! one lower order. After 2 or more failures, H is forced to decrease
2887	! by a factor of 0.2 or less.
2888	!-----------------------------------------------------------------------
2889	500 KFLAG = KFLAG - 1
2890	TN = TOLD
2891	I1 = NQNYH + 1
2892	DO 515 JB = 1,NQ
2893	I1 = I1 - NYH
2894	!dir$ ivdep
2895	DO 510 I = I1,NQNYH
2896	510 YH1(I) = YH1(I) - YH1(I+NYH)
2897	515 CONTINUE
2898	RMAX = 2.0D0
2899	IF (ABS(H) .LE. HMIN*1.00001D0) GO TO 660
2900	IF (KFLAG .LE. -3) GO TO 640
2901	IREDO = 2
2902	RHUP = 0.0D0
2903	GO TO 540
2904	!-----------------------------------------------------------------------
2905	! Regardless of the success or failure of the step, factors
2906	! RHDN, RHSM, and RHUP are computed, by which H could be multiplied
2907	! at order NQ - 1, order NQ, or order NQ + 1, respectively.
2908	! In the case of failure, RHUP = 0.0 to avoid an order increase.
2909	! The largest of these is determined and the new order chosen
2910	! accordingly. If the order is to be increased, we compute one
2911	! additional scaled derivative.
2912	!-----------------------------------------------------------------------
2913	520 RHUP = 0.0D0
2914	IF (L .EQ. LMAX) GO TO 540
2915	DO 530 I = 1,N
2916	530 SAVF(I) = ACOR(I) - YH(I,LMAX)
2917	DUP = DVNORM (N, SAVF, EWT)/TESCO(3,NQ)
2918	EXUP = 1.0D0/(L+1)
2919	RHUP = 1.0D0/(1.4D0DUP*EXUP + 0.0000014D0)
2920	540 EXSM = 1.0D0/L
2921	RHSM = 1.0D0/(1.2D0DSM*EXSM + 0.0000012D0)
2922	RHDN = 0.0D0
2923	IF (NQ .EQ. 1) GO TO 560
2924	DDN = DVNORM (N, YH(1,L), EWT)/TESCO(1,NQ)
2925	EXDN = 1.0D0/NQ
2926	RHDN = 1.0D0/(1.3D0DDN*EXDN + 0.0000013D0)
2927	560 IF (RHSM .GE. RHUP) GO TO 570
2928	IF (RHUP .GT. RHDN) GO TO 590
2929	GO TO 580
2930	570 IF (RHSM .LT. RHDN) GO TO 580
2931	NEWQ = NQ
2932	RH = RHSM
2933	GO TO 620
2934	580 NEWQ = NQ - 1
2935	RH = RHDN
2936	IF (KFLAG .LT. 0 .AND. RH .GT. 1.0D0) RH = 1.0D0
2937	GO TO 620
2938	590 NEWQ = L
2939	RH = RHUP
2940	IF (RH .LT. 1.1D0) GO TO 610
2941	R = EL(L)/L
2942	DO 600 I = 1,N
2943	600 YH(I,NEWQ+1) = ACOR(I)*R
2944	GO TO 630
2945	610 IALTH = 3
2946	GO TO 700
2947	620 IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1D0)) GO TO 610
2948	IF (KFLAG .LE. -2) RH = MIN(RH,0.2D0)
2949	!-----------------------------------------------------------------------
2950	! If there is a change of order, reset NQ, l, and the coefficients.
2951	! In any case H is reset according to RH and the YH array is rescaled.
2952	! Then exit from 690 if the step was OK, or redo the step otherwise.
2953	!-----------------------------------------------------------------------
2954	IF (NEWQ .EQ. NQ) GO TO 170
2955	630 NQ = NEWQ
2956	L = NQ + 1
2957	IRET = 2
2958	GO TO 150
2959	!-----------------------------------------------------------------------
2960	! Control reaches this section if 3 or more failures have occured.
2961	! If 10 failures have occurred, exit with KFLAG = -1.
2962	! It is assumed that the derivatives that have accumulated in the
2963	! YH array have errors of the wrong order. Hence the first
2964	! derivative is recomputed, and the order is set to 1. Then
2965	! H is reduced by a factor of 10, and the step is retried,
2966	! until it succeeds or H reaches HMIN.
2967	!-----------------------------------------------------------------------
2968	640 IF (KFLAG .EQ. -10) GO TO 660
2969	RH = 0.1D0
2970	RH = MAX(HMIN/ABS(H),RH)
2971	H = H*RH
2972	DO 645 I = 1,N
2973	645 Y(I) = YH(I,1)
2974	CALL F (NEQ, TN, Y, SAVF)
2975	NFE = NFE + 1
2976	DO 650 I = 1,N
2977	650 YH(I,2) = H*SAVF(I)
2978	IPUP = MITER
2979	IALTH = 5
2980	IF (NQ .EQ. 1) GO TO 200
2981	NQ = 1
2982	L = 2
2983	IRET = 3
2984	GO TO 150
2985	!-----------------------------------------------------------------------
2986	! All returns are made through this section. H is saved in HOLD
2987	! to allow the caller to change H on the next step.
2988	!-----------------------------------------------------------------------
2989	660 KFLAG = -1
2990	GO TO 720
2991	670 KFLAG = -2
2992	GO TO 720
2993	680 KFLAG = -3
2994	GO TO 720
2995	690 RMAX = 10.0D0
2996	700 R = 1.0D0/TESCO(2,NQU)
2997	DO 710 I = 1,N
2998	710 ACOR(I) = ACOR(I)*R
2999	720 HOLD = H
3000	JSTART = 1
3001	RETURN
3002	!----------------------- END OF SUBROUTINE DSTODE ----------------------
3003	END SUBROUTINE DSTODE
3004	!DECK DEWSET
3005	SUBROUTINE DEWSET (N, ITOL, RelTol, AbsTol, YCUR, EWT)
3006	!***BEGIN PROLOGUE DEWSET
3007	!***SUBSIDIARY
3008	!***PURPOSE Set error weight vector.
3009	!***TYPE KPP_REAL (SEWSET-S, DEWSET-D)
3010	!***AUTHOR Hindmarsh, Alan C., (LLNL)
3011	!***DESCRIPTION
3012	!
3013	! This subroutine sets the error weight vector EWT according to
3014	! EWT(i) = RelTol(i)*ABS(YCUR(i)) + AbsTol(i), i = 1,...,N,
3015	! with the subscript on RelTol and/or AbsTol possibly replaced by 1 above,
3016	! depending on the value of ITOL.
3017	!
3018	!***SEE ALSO DLSODE
3019	!***ROUTINES CALLED (NONE)
3020	!***REVISION HISTORY (YYMMDD)
3021	! 791129 DATE WRITTEN
3022	! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
3023	! 890503 Minor cosmetic changes. (FNF)
3024	! 930809 Renamed to allow single/double precision versions. (ACH)
3025	!***END PROLOGUE DEWSET
3026	!**End
3027	INTEGER N, ITOL
3028	INTEGER I
3029	KPP_REAL RelTol(), AbsTol(), YCUR(N), EWT(N)
3030	!
3031	!***FIRST EXECUTABLE STATEMENT DEWSET
3032	GO TO (10, 20, 30, 40), ITOL
3033	10 CONTINUE
3034	DO 15 I = 1,N
3035	15 EWT(I) = RelTol(1)*ABS(YCUR(I)) + AbsTol(1)
3036	RETURN
3037	20 CONTINUE
3038	DO 25 I = 1,N
3039	25 EWT(I) = RelTol(1)*ABS(YCUR(I)) + AbsTol(I)
3040	RETURN
3041	30 CONTINUE
3042	DO 35 I = 1,N
3043	35 EWT(I) = RelTol(I)*ABS(YCUR(I)) + AbsTol(1)
3044	RETURN
3045	40 CONTINUE
3046	DO 45 I = 1,N
3047	45 EWT(I) = RelTol(I)*ABS(YCUR(I)) + AbsTol(I)
3048	RETURN
3049	!----------------------- END OF SUBROUTINE DEWSET ----------------------
3050	END SUBROUTINE DEWSET
3051	!DECK DVNORM
3052	KPP_REAL FUNCTION DVNORM (N, V, W)
3053	!***BEGIN PROLOGUE DVNORM
3054	!***SUBSIDIARY
3055	!***PURPOSE Weighted root-mean-square vector norm.
3056	!***TYPE KPP_REAL (SVNORM-S, DVNORM-D)
3057	!***AUTHOR Hindmarsh, Alan C., (LLNL)
3058	!***DESCRIPTION
3059	!
3060	! This function routine computes the weighted root-mean-square norm
3061	! of the vector of length N contained in the array V, with weights
3062	! contained in the array W of length N:
3063	! DVNORM = SQRT( (1/N) * SUM( V(i)W(i) )*2 )
3064	!
3065	!***SEE ALSO DLSODE
3066	!***ROUTINES CALLED (NONE)
3067	!***REVISION HISTORY (YYMMDD)
3068	! 791129 DATE WRITTEN
3069	! 890501 Modified prologue to SLATEC/LDOC format. (FNF)
3070	! 890503 Minor cosmetic changes. (FNF)
3071	! 930809 Renamed to allow single/double precision versions. (ACH)
3072	!***END PROLOGUE DVNORM
3073	!**End
3074	INTEGER N, I
3075	KPP_REAL V(N), W(N), SUM
3076	!
3077	!***FIRST EXECUTABLE STATEMENT DVNORM
3078	SUM = 0.0D0
3079	DO 10 I = 1,N
3080	10 SUM = SUM + (V(I)W(I))*2
3081	DVNORM = SQRT(SUM/N)
3082	RETURN
3083	!----------------------- END OF FUNCTION DVNORM ------------------------
3084	END FUNCTION DVNORM
3085	!DECK XERRWD
3086	SUBROUTINE XERRWD (MSG, NMES, NERR, LEVEL, NI, I1, I2, NR, R1, R2)
3087	!***BEGIN PROLOGUE XERRWD
3088	!***SUBSIDIARY
3089	!***PURPOSE Write error message with values.
3090	!***CATEGORY R3C
3091	!***TYPE KPP_REAL (XERRWV-S, XERRWD-D)
3092	!***AUTHOR Hindmarsh, Alan C., (LLNL)
3093	!***DESCRIPTION
3094	!
3095	! Subroutines XERRWD, XSETF, XSETUN, and the function routine IXSAV,
3096	! as given here, constitute a simplified version of the SLATEC error
3097	! handling package.
3098	!
3099	! All arguments are input arguments.
3100	!
3101	! MSG = The message (character array).
3102	! NMES = The length of MSG (number of characters).
3103	! NERR = The error number (not used).
3104	! LEVEL = The error level..
3105	! 0 or 1 means recoverable (control returns to caller).
3106	! 2 means fatal (run is aborted--see note below).
3107	! NI = Number of integers (0, 1, or 2) to be printed with message.
3108	! I1,I2 = Integers to be printed, depending on NI.
3109	! NR = Number of reals (0, 1, or 2) to be printed with message.
3110	! R1,R2 = Reals to be printed, depending on NR.
3111	!
3112	! Note.. this routine is machine-dependent and specialized for use
3113	! in limited context, in the following ways..
3114	! 1. The argument MSG is assumed to be of type CHARACTER, and
3115	! the message is printed with a format of (1X,A).
3116	! 2. The message is assumed to take only one line.
3117	! Multi-line messages are generated by repeated calls.
3118	! 3. If LEVEL = 2, control passes to the statement STOP
3119	! to abort the run. This statement may be machine-dependent.
3120	! 4. R1 and R2 are assumed to be in double precision and are printed
3121	! in D21.13 format.
3122	!
3123	!***ROUTINES CALLED IXSAV
3124	!***REVISION HISTORY (YYMMDD)
3125	! 920831 DATE WRITTEN
3126	! 921118 Replaced MFLGSV/LUNSAV by IXSAV. (ACH)
3127	! 930329 Modified prologue to SLATEC format. (FNF)
3128	! 930407 Changed MSG from CHARACTER*1 array to variable. (FNF)
3129	! 930922 Minor cosmetic change. (FNF)
3130	!***END PROLOGUE XERRWD
3131	!
3132	!*Internal Notes:
3133	!
3134	! For a different default logical unit number, IXSAV (or a subsidiary
3135	! routine that it calls) will need to be modified.
3136	! For a different run-abort command, change the statement following
3137	! statement 100 at the end.
3138	!-----------------------------------------------------------------------
3139	! Subroutines called by XERRWD.. None
3140	! Function routine called by XERRWD.. IXSAV
3141	!-----------------------------------------------------------------------
3142	!**End
3143	!
3144	! Declare arguments.
3145	!
3146	KPP_REAL R1, R2
3147	INTEGER NMES, NERR, LEVEL, NI, I1, I2, NR
3148	CHARACTER() MSG
3149	!
3150	! Declare local variables.
3151	!
3152	INTEGER LUNIT, MESFLG !, IXSAV
3153	!
3154	! Get logical unit number and message print flag.
3155	!
3156	!***FIRST EXECUTABLE STATEMENT XERRWD
3157	LUNIT = IXSAV (1, 0, .FALSE.)
3158	MESFLG = IXSAV (2, 0, .FALSE.)
3159	IF (MESFLG .EQ. 0) GO TO 100
3160	!
3161	! Write the message.
3162	!
3163	WRITE (LUNIT,10) MSG
3164	10 FORMAT(1X,A)
3165	IF (NI .EQ. 1) WRITE (LUNIT, 20) I1
3166	20 FORMAT(6X,'In above message, I1 =',I10)
3167	IF (NI .EQ. 2) WRITE (LUNIT, 30) I1,I2
3168	30 FORMAT(6X,'In above message, I1 =',I10,3X,'I2 =',I10)
3169	IF (NR .EQ. 1) WRITE (LUNIT, 40) R1
3170	40 FORMAT(6X,'In above message, R1 =',D21.13)
3171	IF (NR .EQ. 2) WRITE (LUNIT, 50) R1,R2
3172	50 FORMAT(6X,'In above, R1 =',D21.13,3X,'R2 =',D21.13)
3173	!
3174	! Abort the run if LEVEL = 2.
3175	!
3176	100 IF (LEVEL .NE. 2) RETURN
3177	STOP
3178	!----------------------- End of Subroutine XERRWD ----------------------
3179	END SUBROUTINE XERRWD
3180	!DECK XSETF
3181	SUBROUTINE XSETF (MFLAG)
3182	!***BEGIN PROLOGUE XSETF
3183	!***PURPOSE Reset the error print control flag.
3184	!***CATEGORY R3A
3185	!***TYPE ALL (XSETF-A)
3186	!***KEYWORDS ERROR CONTROL
3187	!***AUTHOR Hindmarsh, Alan C., (LLNL)
3188	!***DESCRIPTION
3189	!
3190	! XSETF sets the error print control flag to MFLAG:
3191	! MFLAG=1 means print all messages (the default).
3192	! MFLAG=0 means no printing.
3193	!
3194	!***SEE ALSO XERRWD, XERRWV
3195	!***REFERENCES (NONE)
3196	!***ROUTINES CALLED IXSAV
3197	!***REVISION HISTORY (YYMMDD)
3198	! 921118 DATE WRITTEN
3199	! 930329 Added SLATEC format prologue. (FNF)
3200	! 930407 Corrected SEE ALSO section. (FNF)
3201	! 930922 Made user-callable, and other cosmetic changes. (FNF)
3202	!***END PROLOGUE XSETF
3203	!
3204	! Subroutines called by XSETF.. None
3205	! Function routine called by XSETF.. IXSAV
3206	!-----------------------------------------------------------------------
3207	!**End
3208	INTEGER MFLAG, JUNK !, IXSAV
3209	!
3210	!***FIRST EXECUTABLE STATEMENT XSETF
3211	IF (MFLAG .EQ. 0 .OR. MFLAG .EQ. 1) JUNK = IXSAV (2,MFLAG,.TRUE.)
3212	RETURN
3213	!----------------------- End of Subroutine XSETF -----------------------
3214	END SUBROUTINE XSETF
3215	!DECK XSETUN
3216	SUBROUTINE XSETUN (LUN)
3217	!***BEGIN PROLOGUE XSETUN
3218	!***PURPOSE Reset the logical unit number for error messages.
3219	!***CATEGORY R3B
3220	!***TYPE ALL (XSETUN-A)
3221	!***KEYWORDS ERROR CONTROL
3222	!***DESCRIPTION
3223	!
3224	! XSETUN sets the logical unit number for error messages to LUN.
3225	!
3226	!***AUTHOR Hindmarsh, Alan C., (LLNL)
3227	!***SEE ALSO XERRWD, XERRWV
3228	!***REFERENCES (NONE)
3229	!***ROUTINES CALLED IXSAV
3230	!***REVISION HISTORY (YYMMDD)
3231	! 921118 DATE WRITTEN
3232	! 930329 Added SLATEC format prologue. (FNF)
3233	! 930407 Corrected SEE ALSO section. (FNF)
3234	! 930922 Made user-callable, and other cosmetic changes. (FNF)
3235	!***END PROLOGUE XSETUN
3236	!
3237	! Subroutines called by XSETUN.. None
3238	! Function routine called by XSETUN.. IXSAV
3239	!-----------------------------------------------------------------------
3240	!**End
3241	INTEGER LUN, JUNK !, IXSAV
3242	!
3243	!***FIRST EXECUTABLE STATEMENT XSETUN
3244	IF (LUN .GT. 0) JUNK = IXSAV (1,LUN,.TRUE.)
3245	RETURN
3246	!----------------------- End of Subroutine XSETUN ----------------------
3247	END SUBROUTINE XSETUN
3248	!DECK IXSAV
3249	INTEGER FUNCTION IXSAV (IPAR, IVALUE, ISET)
3250	!***BEGIN PROLOGUE IXSAV
3251	!***SUBSIDIARY
3252	!***PURPOSE Save and recall error message control parameters.
3253	!***CATEGORY R3C
3254	!***TYPE ALL (IXSAV-A)
3255	!***AUTHOR Hindmarsh, Alan C., (LLNL)
3256	!***DESCRIPTION
3257	!
3258	! IXSAV saves and recalls one of two error message parameters:
3259	! LUNIT, the logical unit number to which messages are printed, and
3260	! MESFLG, the message print flag.
3261	! This is a modification of the SLATEC library routine J4SAVE.
3262	!
3263	! Saved local variables..
3264	! LUNIT = Logical unit number for messages. The default is obtained
3265	! by a call to IUMACH (may be machine-dependent).
3266	! MESFLG = Print control flag..
3267	! 1 means print all messages (the default).
3268	! 0 means no printing.
3269	!
3270	! On input..
3271	! IPAR = Parameter indicator (1 for LUNIT, 2 for MESFLG).
3272	! IVALUE = The value to be set for the parameter, if ISET = .TRUE.
3273	! ISET = Logical flag to indicate whether to read or write.
3274	! If ISET = .TRUE., the parameter will be given
3275	! the value IVALUE. If ISET = .FALSE., the parameter
3276	! will be unchanged, and IVALUE is a dummy argument.
3277	!
3278	! On return..
3279	! IXSAV = The (old) value of the parameter.
3280	!
3281	!***SEE ALSO XERRWD, XERRWV
3282	!***ROUTINES CALLED IUMACH
3283	!***REVISION HISTORY (YYMMDD)
3284	! 921118 DATE WRITTEN
3285	! 930329 Modified prologue to SLATEC format. (FNF)
3286	! 930915 Added IUMACH call to get default output unit. (ACH)
3287	! 930922 Minor cosmetic changes. (FNF)
3288	! 010425 Type declaration for IUMACH added. (ACH)
3289	!***END PROLOGUE IXSAV
3290	!
3291	! Subroutines called by IXSAV.. None
3292	! Function routine called by IXSAV.. IUMACH
3293	!-----------------------------------------------------------------------
3294	!**End
3295	LOGICAL ISET
3296	INTEGER IPAR, IVALUE
3297	!-----------------------------------------------------------------------
3298	INTEGER LUNIT, MESFLG!, IUMACH
3299	!-----------------------------------------------------------------------
3300	! The following Fortran-77 declaration is to cause the values of the
3301	! listed (local) variables to be saved between calls to this routine.
3302	!-----------------------------------------------------------------------
3303	SAVE LUNIT, MESFLG
3304	DATA LUNIT/-1/, MESFLG/1/
3305	!
3306	!***FIRST EXECUTABLE STATEMENT IXSAV
3307	IF (IPAR .EQ. 1) THEN
3308	IF (LUNIT .EQ. -1) LUNIT = IUMACH()
3309	IXSAV = LUNIT
3310	IF (ISET) LUNIT = IVALUE
3311	ENDIF
3312	!
3313	IF (IPAR .EQ. 2) THEN
3314	IXSAV = MESFLG
3315	IF (ISET) MESFLG = IVALUE
3316	ENDIF
3317	!
3318	RETURN
3319	!----------------------- End of Function IXSAV -------------------------
3320	END FUNCTION IXSAV
3321	!DECK IUMACH
3322	INTEGER FUNCTION IUMACH()
3323	!***BEGIN PROLOGUE IUMACH
3324	!***PURPOSE Provide standard output unit number.
3325	!***CATEGORY R1
3326	!***TYPE INTEGER (IUMACH-I)
3327	!***KEYWORDS MACHINE CONSTANTS
3328	!***AUTHOR Hindmarsh, Alan C., (LLNL)
3329	!***DESCRIPTION
3330	! *Usage:
3331	! INTEGER LOUT, IUMACH
3332	! LOUT = IUMACH()
3333	!
3334	! *Function Return Values:
3335	! LOUT : the standard logical unit for Fortran output.
3336	!
3337	!***REFERENCES (NONE)
3338	!***ROUTINES CALLED (NONE)
3339	!***REVISION HISTORY (YYMMDD)
3340	! 930915 DATE WRITTEN
3341	! 930922 Made user-callable, and other cosmetic changes. (FNF)
3342	!***END PROLOGUE IUMACH
3343	!
3344	!*Internal Notes:
3345	! The built-in value of 6 is standard on a wide range of Fortran
3346	! systems. This may be machine-dependent.
3347	!**End
3348	!***FIRST EXECUTABLE STATEMENT IUMACH
3349	IUMACH = 6
3350	!
3351	RETURN
3352	!----------------------- End of Function IUMACH ------------------------
3353	END FUNCTION IUMACH
3354
3355	!---- END OF SUBROUTINE DLSODE AND ITS INTERNAL PROCEDURES
3356	END SUBROUTINE DLSODE
3357	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3358	SUBROUTINE FUN_CHEM(N, T, V, FCT)
3359	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3360
3361	USE KPP_ROOT_Parameters
3362	USE KPP_ROOT_Global
3363	USE KPP_ROOT_Function, ONLY: Fun
3364	USE KPP_ROOT_Rates
3365
3366	IMPLICIT NONE
3367
3368	INTEGER :: N
3369	KPP_REAL :: V(NVAR), FCT(NVAR), T
3370
3371	! TOLD = TIME
3372	! TIME = T
3373	! CALL Update_SUN()
3374	! CALL Update_RCONST()
3375	! CALL Update_PHOTO()
3376	! TIME = TOLD
3377
3378	CALL Fun(V, FIX, RCONST, FCT)
3379
3380	!Nfun=Nfun+1
3381
3382	END SUBROUTINE FUN_CHEM
3383
3384
3385	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3386	SUBROUTINE JAC_CHEM (N, T, V, JF)
3387	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3388
3389	USE KPP_ROOT_Parameters
3390	USE KPP_ROOT_Global
3391	USE KPP_ROOT_JacobianSP
3392	USE KPP_ROOT_Jacobian, ONLY: Jac_SP
3393	USE KPP_ROOT_Rates
3394
3395	IMPLICIT NONE
3396
3397	KPP_REAL :: V(NVAR), T
3398	INTEGER :: N
3399	#ifdef FULL_ALGEBRA
3400	INTEGER :: I, J
3401	KPP_REAL :: JV(LU_NONZERO), JF(NVAR,NVAR)
3402	#else
3403	KPP_REAL :: JF(LU_NONZERO)
3404	#endif
3405
3406	! TOLD = TIME
3407	! TIME = T
3408	! CALL Update_SUN()
3409	! CALL Update_RCONST()
3410	! CALL Update_PHOTO()
3411	! TIME = TOLD
3412
3413	#ifdef FULL_ALGEBRA
3414	CALL Jac_SP(V, FIX, RCONST, JV)
3415	DO j=1,NVAR
3416	DO i=1,NVAR
3417	JF(i,j) = 0.0d0
3418	END DO
3419	END DO
3420	DO i=1,LU_NONZERO
3421	JF(LU_IROW(i),LU_ICOL(i)) = JV(i)
3422	END DO
3423	#else
3424	CALL Jac_SP(V, FIX, RCONST, JF)
3425	#endif
3426	!Njac=Njac+1
3427
3428	END SUBROUTINE JAC_CHEM
3429
3430
3431	END MODULE KPP_ROOT_Integrator

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/kpp_lsode.f90 @ 4416

Download in other formats: