Home

Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

kpp_odessa_ddm.f @ 3534

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

Line
1	SUBROUTINE INTEGRATE( NSENSIT, Y, TIN, TOUT )
2
3	INCLUDE 'KPP_ROOT_Parameters.h'
4	INCLUDE 'KPP_ROOT_Global.h'
5
6	C TIN - Start Time
7	REAL*8 TIN
8	C TOUT - End Time
9	REAL*8 TOUT
10	C Concentrations and Sensitivities
11	REAL*8 Y(NVAR,NSENSIT+1), PARAMS(NSENSIT)
12	C --- Note: Y contains: (1:NVAR) concentrations, followed by
13	C --- (1:NVAR) sensitivities w.r.t. first parameter, followed by
14	C --- etc., followed by
15	C --- (1:NVAR) sensitivities w.r.t. NSENSIT's parameter
16	INTEGER i
17
18	INTEGER LIW, LRW
19	C PARAMETER (LRW = 22 + 8(NSENSIT+1)NVAR + NVAR**2 + NVAR)
20	C PARAMETER (LIW = 21 + NVAR + NSENSIT)
21	C REAL*8 RWORK(LRW)
22	C INTEGER IWORK(LIW)
23	C Note: the following dynamic allocation is not standard F77 and may not work on
24	C some systems. Declare LRW, LIW parameters as above with some upper bound used for NSENSIT
25	REAL8 RWORK(22 + 8(NSENSIT+1)NVAR + NVAR*2 + NVAR)
26	INTEGER IWORK(21 + NVAR + NSENSIT)
27
28	INTEGER IOPT(3), NEQ(2)
29
30	EXTERNAL FUNC_CHEM,JAC,DFUNC_CHEMDPAR
31
32	MF = 21 ! --- BDF plus user-supplied Jacobian
33
34	LRW = 22 + 8(NSENSIT+1)NVAR + NVAR**2 + NVAR
35	LIW = 21 + NVAR + NSENSIT
36
37	NEQ(1) = NVAR ! --- No. of Variables
38	NEQ(2) = NSENSIT ! --- No of parameters
39
40	ITOL=1 ! --- 1=Scalar Tolerances; 4 = VECTOR TOLERANCES
41	ITASK=1 ! --- Normal Output
42	ISTATE=1
43	IOPT(1)=1 ! --- 0= No optional parameters, 1=Optional parameters
44	IOPT(2)=1 ! --- 1=Perform sensitivity analysis; 0 if not
45	IOPT(3)=1 ! --- 1 if DFUNC_CHEMDPAR supplied by the user;
46	! --- 0 if finite differences are to be used
47	C --- Set optional parameters
48	DO 10 i=1,LRW
49	RWORK(i) = 0.0D0
50	10 CONTINUE
51	DO 20 i=1,LIW
52	IWORK(i) = 0
53	20 CONTINUE
54
55	RWORK(5) = STEPMIN ! THE STEP SIZE TO BE ATTEMPTED ON THE FIRST STEP.
56	RWORK(6) = STEPMAX ! THE MAXIMUM ABSOLUTE STEP SIZE ALLOWED.
57	RWORK(7) = 0.0D0 ! THE MINIMUM ABSOLUTE STEP SIZE ALLOWED.
58	IWORK(6) = 5000 ! MAXIMUM NUMBER OF (INTERNALLY DEFINED) STEPS
59
60	CALL KPP_ODESSA( FUNC_CHEM,DFUNC_CHEMDPAR,NEQ,Y,PARAMS,TIN,TOUT,
61	& ITOL,RTOL,ATOL,
62	1 ITASK,ISTATE,IOPT,RWORK,LRW,IWORK,LIW,
63	& JAC,MF)
64
65	IF (ISTATE.LT.0) THEN
66	print *,'KPP_ODESSA: Unsucessfull exit at T=',
67	& TIN,' (ISTATE=',ISTATE,')'
68	ENDIF
69
70	RETURN
71	END
72
73
74
75
76	SUBROUTINE FUNC_CHEM (N, T, V, PARAMS, FCT)
77	INCLUDE 'KPP_ROOT_Parameters.h'
78	INCLUDE 'KPP_ROOT_Global.h'
79
80	DIMENSION V(NVAR), PARAMS(*), FCT(NVAR)
81	TOLD = TIME
82	TIME = T
83	CALL Update_SUN()
84	CALL Update_RCONST()
85	TIME = TOLD
86	CALL Fun(V, FIX, RCONST, FCT)
87	RETURN
88	END
89
90	SUBROUTINE DFUNC_CHEMDPAR (N, T, V, PARAMS, DFCT, JPAR)
91	INCLUDE 'KPP_ROOT_Parameters.h'
92	INCLUDE 'KPP_ROOT_Global.h'
93	C --- NCOEFF = number of rate coefficients w.r.t. which we differentiate
94	C (note that in some applications NCOEFF may be different than NSENSIT)
95	C JCOEFF(1:NCOEFF) are the indices of rate coefficients w.r.t. which we differentiate
96	INTEGER NCOEFF, JCOEFF(NREACT)
97	COMMON /DDMRCOEFF/ NCOEFF, JCOEFF
98
99	DIMENSION V(NVAR), PARAMS(*), DFCT(NVAR)
100	INTEGER JPAR, i, JC(1)
101	IF (DDMTYPE .EQ. 0) THEN
102	C This setting is required for sensitivities w.r.t. initial conditions
103	DO i=1,NVAR
104	DFCT(i) = 0.d0
105	END DO
106	ELSE
107	C This setting is required for sensitivities w.r.t. rate coefficients
108	C ... and should be changed by the user for other applications
109	JC(1) = JCOEFF(JPAR)
110	CALL dFun_dRcoeff(V, FIX, 1, JC, DFCT )
111	END IF
112	RETURN
113	END
114
115	SUBROUTINE JAC (N, T, V, PARAMS, ML, MU, JS, NROWPD)
116	INCLUDE 'KPP_ROOT_Parameters.h'
117	INCLUDE 'KPP_ROOT_Sparse.h'
118
119	REAL8 V(NVAR), PARAMS(), JS(LU_NONZERO)
120	INTEGER ML, MU, NROWPD
121	TOLD = TIME
122	TIME = T
123	CALL Update_SUN()
124	CALL Update_RCONST()
125	TIME = TOLD
126	DO i=1,LU_NONZERO
127	JS(i) = 0.0D0
128	END DO
129	CALL Jac_SP(V, FIX, RCONST, JS)
130	RETURN
131	END
132
133
134	C ALGORITHM 658, COLLECTED ALGORITHMS FROM ACM.
135	C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
136	C VOL. 14, NO. 1, P.61.
137	C-----------------------------------------------------------------------
138	C THIS IS THE SEPTEMBER 1, 1986 VERSION OF ODESSA..
139	C AN ORDINARY DIFFERENTIAL EQUATION KppSolveR WITH EXPLICIT SIMULTANEOUS
140	C SENSITIVITY ANALYSIS.
141	C
142	C THIS PACKAGE IS A MODIFICATION OF THE AUGUST 13, 1981 VERSION OF
143	C LSODE.. LIVERMORE KppSolveR FOR ORDINARY DIFFERENTIAL EQUATIONS.
144	C THIS VERSION IS IN DOUBLE PRECISION.
145	C
146	C ODESSA KppSolveS FOR THE FIRST-ORDER SENSITIVITY COEFFICIENTS..
147	C DY(I)/DP, FOR A SINGLE PARAMETER, OR,
148	C DY(I)/DP(J), FOR MULTIPLE PARAMETERS,
149	C ASSOCIATED WITH A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS..
150	C DY/DT = F(Y,T;P).
151	C-----------------------------------------------------------------------
152	C REFERENCES...
153	C
154	C 1. JORGE R. LEIS AND MARK A. KRAMER, THE SIMULTANEOUS SOLUTION AND
155	C EXPLICIT SENSITIVITY ANALYSIS OF SYSTEMS DESCRIBED BY ORDINARY
156	C DIFFERENTIAL EQUATIONS. SUBMITTED TO ACM TRANS. MATH. SOFTWARE,
157	C (1985).
158	C
159	C 2. JORGE R. LEIS AND MARK A. KRAMER, ODESSA - AN ORDINARY DIFFERENTIA
160	C EQUATION KppSolveR WITH EXPLICIT SIMULTANEOUS SENSITIVITY ANALYSIS.
161	C SUBMITTED TO ACM TRANS. MATH. SOFTWARE, (1985).
162	C
163	C 3. ALAN C. HINDMARSH, LSODE AND LSODI, TWO NEW INITIAL VALUE
164	C ORDINARY DIFFERENTIAL EQUATION KppSolveRS, ACM-SIGNUM NEWSLETTER,
165	C VOL. 15, NO. 4 (1980), PP. 10-11.
166	C-----------------------------------------------------------------------
167	C PROBLEM STATEMENT..
168	C
169	C THE ODESSA MODIFICATION OF THE LSODE PACKAGE PROVIDES THE OPTION TO
170	C CALCULATE FIRST-ORDER SENSITIVITY COEFFICIENTS FOR A SYSTEM OF STIFF
171	C OR NON-STIFF EXPLICIT ORDINARY DIFFERENTIAL EQUATIONS OF THE GENERAL
172	C FORM :
173	C
174	C DY/DT = F(Y,T;P) (1)
175	C
176	C WHERE Y IS AN N-DIMENSIONAL DEPENDENT VARIABLE VECTOR, T IS THE
177	C INDEPENDENT INTEGRATION VARIABLE, AND P IS AN NPAR-DIMENSIONAL
178	C CONSTANT VECTOR. THE GOVERNING EQUATIONS FOR THE FIRST-ORDER
179	C SENSITIVITY COEFFICIENTS ARE GIVEN BY :
180	C
181	C S'(T) = J(T)*S(T) + DF/DP (2)
182	C
183	C WHERE
184	C
185	C S(T) = DY(T)/DP (= SENSITIVITY FUNCTIONS)
186	C S'(T) = D(DY(T)/DP)/DT
187	C J(T) = DF(Y,T;P)/DY(T) (= JACOBIAN MATRIX)
188	C AND DF/DP = DF(Y,T;P)/DP (= INHOMOGENEITY MATRIX)
189	C
190	C SOLUTION OF EQUATIONS (1) AND (2) PROCEEDS SIMULTANEOUSLY VIA AN
191	C EXTENSION OF THE LSODE PACKAGE AS DESCRIBED IN [1].
192	C----------------------------------------------------------------------
193	C ACKNOWLEDGEMENT : THE FOLLOWING ODESSA PACKAGE DOCUMENTATION IS A
194	C MODIFICATION OF THE LSODE DOCUMENTATION WHICH
195	C ACCOMPANIES THE LSODE PACKAGE CODE.
196	C----------------------------------------------------------------------
197	C SUMMARY OF USAGE.
198	C
199	C COMMUNICATION BETWEEN THE USER AND THE ODESSA PACKAGE, FOR NORMAL
200	C SITUATIONS, IS SUMMARIZED HERE. THIS SUMMARY DESCRIBES ONLY A SUBSET
201	C OF THE FULL SET OF OPTIONS AVAILABLE. SEE THE FULL DESCRIPTION FOR
202	C DETAILS, INCLUDING OPTIONAL COMMUNICATION, NONSTANDARD OPTIONS,
203	C AND INSTRUCTIONS FOR SPECIAL SITUATIONS. SEE ALSO THE EXAMPLE
204	C PROBLEM (WITH PROGRAM AND OUTPUT) FOLLOWING THIS SUMMARY.
205	C
206	C A. FIRST PROVIDE A SUBROUTINE OF THE FORM..
207	C SUBROUTINE F (N, T, Y, PAR, YDOT)
208	C DOUBLE PRECISION T, Y, PAR, YDOT
209	C DIMENSION Y(N), YDOT(N), PAR(NPAR)
210	C WHICH SUPPLIES THE VECTOR FUNCTION F BY LOADING YDOT(I) WITH F(I).
211	C N IS THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS IN THE
212	C ABOVE MODEL. NPAR IS THE NUMBER OF MODEL PARAMETERS FOR WHICH
213	C VECTOR SENSITIVITY FUNCTIONS ARE DESIRED. YOU ARE ALSO ENCOURAGED
214	C TO PROVIDE A SUBROUTINE OF THE FORM..
215	C SUBROUTINE DF (N, T, Y, PAR, DFDP, JPAR)
216	C DOUBLE PRECISION T, Y, PAR, DFDP
217	C DIMENSION Y(N), PAR(NPAR), DFDP(N)
218	C GO TO (1,...,NPAR) JPAR
219	C 1 DFDP(1) = DF(1)/DP(1)
220	C .
221	C DFDP(I) = DF(I)/DP(1)
222	C .
223	C DFDP(N) = DF(N)/DP(1)
224	C RETURN
225	C 2 DFDP(1) = DF(1)/DP(2)
226	C .
227	C DFDP(I) = DF(I)/DP(2)
228	C .
229	C DFDP(N) = DF(N)/DP(2)
230	C RETURN
231	C . .
232	C . .
233	C RETURN
234	C NPAR DFDP(1) = DF(1)/DP(NPAR)
235	C .
236	C DFDP(I) = DF(I)/DP(NPAR)
237	C .
238	C DFDP(N) = DF(N)/DP(NPAR)
239	C RETURN
240	C END
241	C ONLY NONZERO ELEMENTS NEED BE LOADED. IF THIS IS NOT FEASIBLE,
242	C ODESSA WILL GENERATE THIS MATRIX INTERNALLY BY DIFFERENCE QUOTIENTS.
243	C
244	C B. NEXT DETERMINE (OR GUESS) WHETHER OR NOT THE PROBLEM IS STIFF.
245	C STIFFNESS OCCURS WHEN THE JACOBIAN MATRIX DF/DY HAS AN EIGENVALUE
246	C WHOSE REAL PART IS NEGATIVE AND LARGE IN MAGNITUDE, COMPARED TO THE
247	C RECIPROCAL OF THE T SPAN OF INTEREST. IF THE PROBLEM IS NONSTIFF,
248	C USE METH = 10. IF IT IS STIFF, USE METH = 20. THE USER IS REQUIRED
249	C TO INPUT THE METHOD FLAG MF = 10*METH + MITER. THERE ARE FOUR
250	C STANDARD CHOICES FOR MITER WHEN A SENSITIVITY ANALYSIS IS DESIRED,
251	C AND ODESSA REQUIRES THE JACOBIAN MATRIX IN SOME FORM.
252	C THIS MATRIX IS REGARDED EITHER AS FULL (MITER = 1 OR 2),
253	C OR BANDED (MITER = 4 OR 5). IN THE BANDED CASE, ODESSA REQUIRES TWO
254	C HALF-BANDWIDTH PARAMETERS ML AND MU. THESE ARE, RESPECTIVELY, THE
255	C WIDTHS OF THE LOWER AND UPPER PARTS OF THE BAND, EXCLUDING THE MAIN
256	C DIAGONAL. THUS THE BAND CONSISTS OF THE LOCATIONS (I,J) WITH
257	C I-ML .LE. J .LE. I+MU, AND THE FULL BANDWIDTH IS ML+MU+1.
258	C
259	C C. YOU ARE ENCOURAGED TO SUPPLY THE JACOBIAN DIRECTLY (MF = 11, 14,
260	C 21, OR 24), BUT IF THIS IS NOT FEASIBLE, ODESSA WILL COMPUTE IT
261	C INTERNALLY BY DIFFERENCE QUOTIENTS (MF = 12, 15, 22, OR 25). IF YOU
262	C ARE SUPPLYING THE JACOBIAN, PROVIDE A SUBROUTINE OF THE FORM..
263	C SUBROUTINE JAC (NEQ, T, Y, PAR, ML, MU, PD, NROWPD)
264	C DOUBLE PRECISION T, Y, PAR, PD
265	C DIMENSION Y(N), PD(NROWPD,N), PAR(NPAR)
266	C WHICH SUPPLIES DF/DY BY LOADING PD AS FOLLOWS..
267	C FOR A FULL JACOBIAN (MF = 11, OR 21), LOAD PD(I,J) WITH DF(I)/DY(J),
268	C THE PARTIAL DERIVATIVE OF F(I) WITH RESPECT TO Y(J). (IGNORE THE
269	C ML AND MU ARGUMENTS IN THIS CASE.)
270	C FOR A BANDED JACOBIAN (MF = 14, OR 24), LOAD PD(I-J+MU+1,J) WITH
271	C DF(I)/DY(J), I.E. LOAD THE DIAGONAL LINES OF DF/DY INTO THE ROWS OF
272	C PD FROM THE TOP DOWN.
273	C IN EITHER CASE, ONLY NONZERO ELEMENTS NEED BE LOADED.
274	C
275	C D. WRITE A MAIN PROGRAM WHICH CALLS SUBROUTINE ODESSA ONCE FOR
276	C EACH POINT AT WHICH ANSWERS ARE DESIRED. THIS SHOULD ALSO PROVIDE
277	C FOR POSSIBLE USE OF LOGICAL UNIT 6 FOR OUTPUT OF ERROR MESSAGES BY
278	C ODESSA. ON THE FIRST CALL TO ODESSA, SUPPLY ARGUMENTS AS FOLLOWS..
279	C F = NAME OF SUBROUTINE FOR RIGHT-HAND SIDE VECTOR F (MODEL).
280	C THIS NAME MUST BE DECLARED EXTERNAL IN CALLING PROGRAM.
281	C DF = NAME OF SUBROUTINE FOR INHOMOGENEITY MATRIX DF/DP.
282	C IF USED (IDF = 1), THIS NAME MUST BE DECLARED EXTERNAL IN
283	C CALLING PROGRAM. IF NOT USED (IDF = 0), PASS A DUMMY NAME.
284	C N = NUMBER OF FIRST ORDER ODE-S IN MODEL; LOAD INTO NEQ(1).
285	C NPAR = NUMBER OF MODEL PARAMETERS OF INTEREST; LOAD INTO NEQ(2).
286	C Y = AN (N) BY (NPAR+1) REAL ARRAY OF INITIAL VALUES..
287	C Y(I,1) , I = 1,N , CONTAIN THE STATE, OR MODEL, DEPENDENT
288	C VARIABLES,
289	C Y(I,J) , J = 2,NPAR , CONTAIN THE DEPENDENT SENSITIVITY
290	C COEFFICIENTS.
291	C PAR = A REAL ARRAY OF LENGTH NPAR CONTAINING MODEL PARAMETERS
292	C OF INTEREST.
293	C T = THE INITIAL VALUE OF THE INDEPENDENT VARIABLE.
294	C TOUT = FIRST POINT WHERE OUTPUT IS DESIRED (.NE. T).
295	C ITOL = 1, 2, 3, OR 4 ACCORDING AS RTOL, ATOL (BELOW) ARE SCALARS
296	C OR ARRAYS.
297	C RTOL = RELATIVE TOLERANCE PARAMETER (SCALAR OR (N) BY (NPAR+1)
298	C ARRAY).
299	C ATOL = ABSOLUTE TOLERANCE PARAMETER (SCALAR OR (N) BY (NPAR+1)
300	C ARRAY).
301	C THE ESTIMATED LOCAL ERROR IN Y(I,J) WILL BE CONTROLLED SO AS
302	C TO BE ROUGHLY LESS (IN MAGNITUDE) THAN
303	C EWT(I,J) = RTOL*ABS(Y(I,J)) + ATOL IF ITOL = 1,
304	C EWT(I,J) = RTOL*ABS(Y(I,J)) + ATOL(I,J) IF ITOL = 2,
305	C EWT(I,J) = RTOL(I,J)*ABS(Y(I,J) + ATOL IF ITOL = 3, OR
306	C EWT(I,J) = RTOL(I,J)*ABS(Y(I,J) + ATOL(I,J) IF ITOL = 4.
307	C THUS THE LOCAL ERROR TEST PASSES IF, IN EACH COMPONENT,
308	C EITHER THE ABSOLUTE ERROR IS LESS THAN ATOL (OR ATOL(I,J)),
309	C OR THE RELATIVE ERROR IS LESS THAN RTOL (OR RTOL(I,J)).
310	C USE RTOL = 0.0 FOR PURE ABSOLUTE ERROR CONTROL, AND
311	C USE ATOL = 0.0 FOR PURE RELATIVE ERROR CONTROL.
312	C CAUTION.. ACTUAL (GLOBAL) ERRORS MAY EXCEED THESE LOCAL
313	C TOLERANCES, SO CHOOSE THEM CONSERVATIVELY.
314	C ITASK = 1 FOR NORMAL COMPUTATION OF OUTPUT VALUES OF Y AT T = TOUT.
315	C ISTATE = INTEGER FLAG (INPUT AND OUTPUT). SET ISTATE = 1.
316	C IOPT = 0, TO INDICATE NO OPTIONAL INPUTS FOR INTEGRATION;
317	C LOAD INTO IOPT(1).
318	C ISOPT = 1, TO INDICATE SENSITIVITY ANALYSIS, = 0, TO INDICATE
319	C NO SENSITIVITY ANALYSIS; LOAD INTO IOPT(2).
320	C IDF = 1, IF SUBROUTINE DF (ABOVE) IS SUPPLIED BY THE USER,
321	C = 0, OTHERWISE; LOAD INTO IOPT(3).
322	C RWORK = REAL WORK ARRAY OF LENGTH AT LEAST..
323	C 22 + 16N + N*2 FOR MF = 11 OR 12,
324	C 22 + 17N + (2ML + MU)*N FOR MF = 14 OR 15,
325	C 22 + 9N + N*2 FOR MF = 21 OR 22,
326	C 22 + 10N + (2ML + MU)*N FOR MF = 24 OR 25,
327	C IF ISOPT = 0, OR..
328	C 22 + 15(NPAR+1)N + N**2 + N FOR MF = 11 OR 12,
329	C 24 + 15(NPAR+1)N + (2ML+MU+2)N + N FOR MF = 14 OR 15,
330	C 22 + 8(NPAR+1)N + N**2 + N FOR MF = 21 OR 22,
331	C 24 + 8(NPAR+1)N + (2ML+MU+2)N + N FOR MF = 24 OR 25,
332	C IF ISOPT = 1.
333	C LRW = DECLARED LENGTH OF RWORK (IN USER-S DIMENSION STATEMENT).
334	C IWORK = INTEGER WORK ARRAY OF LENGTH AT LEAST..
335	C 20 + N IF ISOPT = 0,
336	C 21 + N + NPAR IF ISOPT = 1.
337	C IF MITER = 4 OR 5, INPUT IN IWORK(1),IWORK(2) THE LOWER
338	C AND UPPER HALF-BANDWIDTHS ML,MU (EXCLUDING MAIN DIAGONAL).
339	C LIW = DECLARED LENGTH OF IWORK (IN USER-S DIMENSION STATEMENT).
340	C JAC = NAME OF SUBROUTINE FOR JACOBIAN MATRIX.
341	C IF USED, THIS NAME MUST BE DECLARED EXTERNAL IN CALLING
342	C PROGRAM. IF NOT USED, PASS A DUMMY NAME.
343	C MF = METHOD FLAG. STANDARD VALUES FOR ISOPT = 0 ARE..
344	C 10 FOR NONSTIFF (ADAMS) METHOD, NO JACOBIAN USED.
345	C 21 FOR STIFF (BDF) METHOD, USER-SUPPLIED FULL JACOBIAN.
346	C 22 FOR STIFF METHOD, INTERNALLY GENERATED FULL JACOBIAN.
347	C 24 FOR STIFF METHOD, USER-SUPPLIED BANDED JACOBIAN.
348	C 25 FOR STIFF METHOD, INTERNALLY GENERATED BANDED JACOBIAN.
349	C IF ISOPT = 1, MF = 10 IS ILLEGAL AND CAN BE REPLACED BY..
350	C 11 FOR NONSTIFF METHOD, USER-SUPPLIED FULL JACOBIAN.
351	C 12 FOR NONSTIFF METHOD, INTERNALLY GENERATED FULL JACOBIAN.
352	C 14 FOR NONSTIFF METHOD, USER-SUPPLIED BANDED JACOBIAN.
353	C 15 FOR NONSTIFF METHOD, INTERNALLY GENERATED BANDED JACOBIAN.
354	C NOTE THAT THE MAIN PROGRAM MUST DECLARE ARRAYS Y, RWORK, IWORK, AND
355	C POSSIBLY ATOL AND RTOL, AS WELL AS NEQ, IOPT, AND PAR IF ISOPT = 1.
356	C
357	C E. THE OUTPUT FROM THE FIRST CALL (OR ANY CALL) IS..
358	C Y = ARRAY OF COMPUTED VALUES OF Y(T) VECTOR.
359	C T = CORRESPONDING VALUE OF INDEPENDENT VARIABLE (NORMALLY TOUT).
360	C ISTATE = 2 IF ODESSA WAS SUCCESSFUL, NEGATIVE OTHERWISE.
361	C -1 MEANS EXCESS WORK DONE ON THIS CALL (PERHAPS WRONG MF).
362	C -2 MEANS EXCESS ACCURACY REQUESTED (TOLERANCES TOO SMALL).
363	C -3 MEANS ILLEGAL INPUT DETECTED (SEE PRINTED MESSAGE).
364	C -4 MEANS REPEATED ERROR TEST FAILURES (CHECK ALL INPUTS).
365	C -5 MEANS REPEATED CONVERGENCE FAILURES (PERHAPS BAD JACOBIAN
366	C SUPPLIED OR WRONG CHOICE OF MF OR TOLERANCES).
367	C -6 MEANS ERROR WEIGHT BECAME ZERO DURING PROBLEM. (SOLUTION
368	C COMPONENT I,J VANISHED, AND ATOL OR ATOL(I,J) = 0.0)
369	C
370	C F. TO CONTINUE THE INTEGRATION AFTER A SUCCESSFUL RETURN, SIMPLY
371	C RESET TOUT AND CALL ODESSA AGAIN. NO OTHER PARAMETERS NEED BE RESET.
372	C----------------------------------------------------------------------
373	C EXAMPLE PROBLEM.
374	C
375	C THE FOLLOWING IS A SIMPLE EXAMPLE PROBLEM, WITH THE CODING
376	C NEEDED FOR ITS SOLUTION BY ODESSA. THE PROBLEM IS FROM CHEMICAL
377	C KINETICS, AND CONSISTS OF THE FOLLOWING THREE RATE EQUATIONS..
378	C DY1/DT = -PAR(1)Y1 + PAR(2)Y2*Y3 ; PAR(1) = .04, PAR(2) = 1.E4
379	C DY2/DT = PAR(1)Y1 - PAR(2)Y2Y3 - PAR(3)Y2**2 ; PAR(3) = 3.E7
380	C DY3/DT = PAR(3)Y2*2
381	C ON THE INTERVAL FROM T = 0.0 TO T = 4.E10, WITH INITIAL CONDITIONS
382	C Y1 = 1.0, Y2 = Y3 = 0, AND S(I,J) = 0, I = 1,3, J = 1,3.
383	C THE PROBLEM IS STIFF.
384	C
385	C THE FOLLOWING CODING KppSolveS THIS PROBLEM WITH ODESSA, USING
386	C MF = 21 AND PRINTING RESULTS AT T = .4, 4., ..., 4.E10.
387	C IT USES ITOL = 4 AND ATOL MUCH SMALLER FOR Y2 THAN Y1 OR Y3,
388	C BECAUSE Y2 HAS MUCH SMALLER VALUES. LESS STRINGENT TOLERANCES
389	C ARE ASSIGNED FOR THE SENSITIVITIES TO ACHIEVE GREATER EFFICIENCY.
390	C AT THE END OF THE RUN, STATISTICAL QUANTITIES OF INTEREST ARE
391	C PRINTED (SEE OPTIONAL OUTPUTS IN THE FULL DESCRIPTION BELOW).
392	C
393	C DOUBLE PRECISION ATOL, RWORK, RTOL, T, TOUT, Y, PAR
394	C EXTERNAL FEX, JEX, DFEX
395	C DIMENSION Y(3,4), PAR(3), ATOL(3,4), RTOL(3,4), RWORK(130),
396	C 1 IWORK(27), NEQ(2), IOPT(3)
397	C N = 3
398	C NPAR = 3
399	C NEQ(1) = N
400	C NEQ(2) = NPAR
401	C NSV = NPAR+1
402	C DO 10 I = 1,N
403	C DO 10 J = 1,NSV
404	C 10 Y(I,J) = 0.0D0
405	C Y(1,1) = 1.0D0
406	C PAR(1) = 0.04D0
407	C PAR(2) = 1.0D4
408	C PAR(3) = 3.0D7
409	C T = 0.D0
410	C TOUT = .4D0
411	C ITOL = 4
412	C ATOL(1,1) = 1.D-6
413	C ATOL(2,1) = 1.D-10
414	C ATOL(3,1) = 1.D-6
415	C DO 20 I = 1,N
416	C RTOL(I,1) = 1.D-4
417	C DO 15 J = 2,NSV
418	C RTOL(I,J) = 1.D-3
419	C 15 ATOL(I,J) = 1.D2 * ATOL(I,1)
420	C 20 CONTINUE
421	C ITASK = 1
422	C ISTATE = 1
423	C IOPT(1) = 0
424	C IOPT(2) = 1
425	C IOPT(3) = 1
426	C LRW = 130
427	C LIW = 27
428	C MF = 21
429	C DO 60 IOUT = 1,12
430	C CALL ODESSA(FEX,DFEX,NEQ,Y,PAR,T,TOUT,ITOL,RTOL,ATOL,
431	C 1 ITASK,ISTATE, IOPT,RWORK,LRW,IWORK,LIW,JEX,MF)
432	C WRITE(6,30)T,Y(1,1),Y(2,1),Y(3,1)
433	C 30 FORMAT(1X,7H AT T =,E12.4,6H Y =,3E14.6)
434	C DO 50 J = 2,NSV
435	C JPAR = J-1
436	C WRITE(6,40)JPAR,Y(1,J),Y(2,J),Y(3,J)
437	C 40 FORMAT(20X,2HS(,I1,3H) =,3E14.6)
438	C 50 CONTINUE
439	C IF (ISTATE .LT. 0) GO TO 80
440	C 60 TOUT = TOUT*10.D0
441	C WRITE(6,70)IWORK(11),IWORK(12),IWORK(13),IWORK(19)
442	C 70 FORMAT(1X,/,12H NO. STEPS =,I4,11H NO. F-S =,I4,11H NO. J-S =,
443	C 1 I4,12H NO. DF-S =,I4)
444	C STOP
445	C 80 WRITE(6,90)ISTATE
446	C 90 FORMAT(///22H ERROR HALT.. ISTATE =,I3)
447	C STOP
448	C END
449	C
450	C SUBROUTINE FEX (NEQ, T, Y, PAR, YDOT)
451	C DOUBLE PRECISION T, Y, YDOT, PAR
452	C DIMENSION Y(3), YDOT(3), PAR(3)
453	C YDOT(1) = -PAR(1)Y(1) + PAR(2)Y(2)*Y(3)
454	C YDOT(3) = PAR(3)Y(2)Y(2)
455	C YDOT(2) = -YDOT(1) - YDOT(3)
456	C RETURN
457	C END
458	C
459	C SUBROUTINE JEX (NEQ, T, Y, PAR, ML, MU, PD, NRPD)
460	C DOUBLE PRECISION PD, T, Y, PAR
461	C DIMENSION Y(3), PD(NRPD,3), PAR(3)
462	C PD(1,1) = -PAR(1)
463	C PD(1,2) = PAR(2)*Y(3)
464	C PD(1,3) = PAR(2)*Y(2)
465	C PD(2,1) = PAR(1)
466	C PD(2,3) = -PD(1,3)
467	C PD(3,2) = 2.D0PAR(3)Y(2)
468	C PD(2,2) = -PD(1,2) - PD(3,2)
469	C RETURN
470	C END
471	C
472	C SUBROUTINE DFEX (NEQ, T, Y, PAR, DFDP, JPAR)
473	C DOUBLE PRECISION T, Y, PAR, DFDP
474	C DIMENSION Y(3), PAR(3), DFDP(3)
475	C GO TO (1,2,3), JPAR
476	C 1 DFDP(1) = -Y(1)
477	C DFDP(2) = Y(1)
478	C RETURN
479	C 2 DFDP(1) = Y(2)*Y(3)
480	C DFDP(2) = -Y(2)*Y(3)
481	C RETURN
482	C 3 DFDP(2) = -Y(2)*Y(2)
483	C DFDP(3) = Y(2)*Y(2)
484	C RETURN
485	C END
486	C
487	C THE OUTPUT OF THIS PROGRAM (ON A DATA GENERAL MV-8000 IN
488	C DOUBLE PRECISION IS AS FOLLOWS:
489	C
490	C AT T = .4000E+00 Y = .985173E+00 .338641E-04 .147930E-01
491	C S(1) = -.355914E+00 .390261E-03 .355524E+00
492	C S(2) = .955150E-07 -.213065E-09 -.953019E-07
493	C S(3) = -.158466E-10 -.529012E-12 .163756E-10
494	C AT T = .4000E+01 Y = .905516E+00 .224044E-04 .944615E-01
495	C S(1) = -.187621E+01 .179197E-03 .187603E+01
496	C S(2) = .296093E-05 -.583104E-09 -.296034E-05
497	C S(3) = -.493267E-09 -.276246E-12 .493544E-09
498	C AT T = .4000E+02 Y = .715848E+00 .918628E-05 .284143E+00
499	C S(1) = -.424730E+01 .459360E-04 .424726E+01
500	C S(2) = .137294E-04 -.235815E-09 -.137291E-04
501	C S(3) = -.228818E-08 -.113803E-12 .228829E-08
502	C AT T = .4000E+03 Y = .450526E+00 .322299E-05 .549471E+00
503	C S(1) = -.595837E+01 .354310E-05 .595836E+01
504	C S(2) = .227380E-04 -.226041E-10 -.227380E-04
505	C S(3) = -.378971E-08 -.499501E-13 .378976E-08
506	C AT T = .4000E+04 Y = .183185E+00 .894131E-06 .816814E+00
507	C S(1) = -.475006E+01 -.599504E-05 .475007E+01
508	C S(2) = .188089E-04 .231330E-10 -.188089E-04
509	C S(3) = -.313478E-08 -.187575E-13 .313480E-08
510	C AT T = .4000E+05 Y = .389733E-01 .162133E-06 .961027E+00
511	C S(1) = -.157477E+01 -.276199E-05 .157477E+01
512	C S(2) = .628668E-05 .110026E-10 -.628670E-05
513	C S(3) = -.104776E-08 -.453588E-14 .104776E-08
514	C AT T = .4000E+06 Y = .493609E-02 .198411E-07 .995064E+00
515	C S(1) = -.236244E+00 -.458262E-06 .236244E+00
516	C S(2) = .944669E-06 .183193E-11 -.944671E-06
517	C S(3) = -.157441E-09 -.635990E-15 .157442E-09
518	C AT T = .4000E+07 Y = .516087E-03 .206540E-08 .999484E+00
519	C S(1) = -.256277E-01 -.509808E-07 .256278E-01
520	C S(2) = .102506E-06 .203905E-12 -.102506E-06
521	C S(3) = -.170825E-10 -.684002E-16 .170826E-10
522	C AT T = .4000E+08 Y = .519314E-04 .207736E-09 .999948E+00
523	C S(1) = -.259316E-02 -.518029E-08 .259316E-02
524	C S(2) = .103726E-07 .207209E-13 -.103726E-07
525	C S(3) = -.172845E-11 -.691450E-17 .172845E-11
526	C AT T = .4000E+09 Y = .544710E-05 .217885E-10 .999995E+00
527	C S(1) = -.271637E-03 -.541849E-09 .271638E-03
528	C S(2) = .108655E-08 .216739E-14 -.108655E-08
529	C S(3) = -.180902E-12 -.723615E-18 .180902E-12
530	C AT T = .4000E+10 Y = .446748E-06 .178699E-11 .100000E+01
531	C S(1) = -.322322E-04 -.842541E-10 .322323E-04
532	C S(2) = .128929E-09 .337016E-15 -.128929E-09
533	C S(3) = -.209715E-13 -.838859E-19 .209715E-13
534	C AT T = .4000E+11 Y = -.363960E-07 -.145584E-12 .100000E+01
535	C S(1) = -.164109E-06 -.429604E-11 .164113E-06
536	C S(2) = .656436E-12 .171842E-16 -.656451E-12
537	C S(3) = -.689361E-15 -.275745E-20 .689363E-15
538	C
539	C NO. STEPS = 340 NO. F-S = 412 NO. J-S = 343 NO. DF-S =1023
540	C----------------------------------------------------------------------
541	C FULL DESCRIPTION OF USER INTERFACE TO ODESSA.
542	C
543	C THE USER INTERFACE TO ODESSA CONSISTS OF THE FOLLOWING PARTS.
544	C
545	C I. THE CALL SEQUENCE TO SUBROUTINE ODESSA, WHICH IS A DRIVER
546	C ROUTINE FOR THE KppSolveR. THIS INCLUDES DESCRIPTIONS OF BOTH
547	C THE CALL SEQUENCE ARGUMENTS AND OF USER-SUPPLIED ROUTINES.
548	C FOLLOWING THESE DESCRIPTIONS IS A DESCRIPTION OF
549	C OPTIONAL INPUTS AVAILABLE THROUGH THE CALL SEQUENCE, AND THEN
550	C A DESCRIPTION OF OPTIONAL OUTPUTS (IN THE WORK ARRAYS).
551	C
552	C II. DESCRIPTIONS OF OTHER ROUTINES IN THE ODESSA PACKAGE THAT MAY
553	C BE (OPTIONALLY) CALLED BY THE USER. THESE PROVIDE THE ABILITY
554	C TO ALTER ERROR MESSAGE HANDLING, SAVE AND RESTORE THE INTERNAL
555	C COMMON, AND OBTAIN SPECIFIED DERIVATIVES OF THE SOLUTION Y(T).
556	C
557	C III. DESCRIPTIONS OF COMMON BLOCKS TO BE DECLARED IN OVERLAY
558	C OR SIMILAR ENVIRONMENTS, OR TO BE SAVED WHEN DOING AN INTERRUPT
559	C OF THE PROBLEM AND CONTINUED SOLUTION LATER.
560	C
561	C IV. DESCRIPTION OF TWO SUBROUTINES IN THE ODESSA PACKAGE, EITHER OF
562	C WHICH THE USER MAY REPLACE WITH HIS OWN VERSION, IF DESIRED.
563	C THESE RELATE TO THE MEASUREMENT OF ERRORS.
564	C
565	C V. GENERAL REMARKS WHICH HIGHLIGHT DIFFERENCES BETWEEN THE LSODE
566	C PACKAGE AND THE ODESSA PACKAGE.
567	C----------------------------------------------------------------------
568	C PART I. CALL SEQUENCE.
569	C
570	C THE CALL SEQUENCE PARAMETERS USED FOR INPUT ONLY ARE..
571	C F, DF, NEQ, PAR, TOUT, ITOL, RTOL, ATOL, ITASK, IOPT, LRW, LIW,
572	C JAC, MF,
573	C AND THOSE USED FOR BOTH INPUT AND OUTPUT ARE
574	C Y, T, ISTATE.
575	C THE WORK ARRAYS RWORK AND IWORK ARE ALSO USED FOR CONDITIONAL AND
576	C OPTIONAL INPUTS AND OPTIONAL OUTPUTS. (THE TERM OUTPUT HERE REFERS
577	C TO THE RETURN FROM SUBROUTINE ODESSA TO THE USER-S CALLING PROGRAM.)
578	C
579	C THE LEGALITY OF INPUT PARAMETERS WILL BE THOROUGHLY CHECKED ON THE
580	C INITIAL CALL FOR THE PROBLEM, BUT NOT CHECKED THEREAFTER UNLESS A
581	C CHANGE IN INPUT PARAMETERS IS FLAGGED BY ISTATE = 3 ON INPUT.
582	C
583	C THE DESCRIPTIONS OF THE CALL ARGUMENTS ARE AS FOLLOWS.
584	C
585	C F = THE NAME OF THE USER-SUPPLIED SUBROUTINE DEFINING THE
586	C ODE MODEL. THIS SYSTEM MUST BE PUT IN THE FIRST-ORDER
587	C FORM DY/DT = F(Y,T;P), WHERE F IS A VECTOR-VALUED FUNCTION
588	C OF THE SCALAR T AND VECTORS Y, AND PAR. SUBROUTINE F IS TO
589	C COMPUTE THE FUNCTION F. IT IS TO HAVE THE FORM..
590	C SUBROUTINE F (NEQ, T, Y, PAR, YDOT)
591	C DOUBLE PRECISION T, Y, PAR, YDOT
592	C DIMENSION Y(1), PAR(1), YDOT(1)
593	C WHERE NEQ, T, Y, AND PAR ARE INPUT, AND YDOT = F(Y,T;P)
594	C IS OUTPUT. Y AND YDOT ARE ARRAYS OF LENGTH N (= NEQ(1)).
595	C (IN THE DIMENSION STATEMENT ABOVE, 1 IS A DUMMY
596	C DIMENSION.. IT CAN BE REPLACED BY ANY VALUE.)
597	C F SHOULD NOT ALTER ARRAY Y, OR PAR(1),...,PAR(NPAR).
598	C F MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM.
599	C
600	C SUBROUTINE F MAY ACCESS USER-DEFINED QUANTITIES IN
601	C NEQ(2),... AND PAR(NPAR+1),... IF NEQ IS AN ARRAY
602	C (DIMENSIONED IN F) AND PAR HAS LENGTH EXCEEDING NPAR.
603	C SEE THE DESCRIPTIONS OF NEQ AND PAR BELOW.
604	C
605	C DF = THE NAME OF THE USER-SUPPLIED ROUTINE (IDF = 1) TO COMPUTE
606	C THE INHOMOGENEITY MATRIX, DF/DP, AS A FUNCTION OF THE SCALAR
607	C T, AND THE VECTORS Y, AND PAR. IT IS TO HAVE THE FORM
608	C SUBROUTINE DF (NEQ, T, Y, PAR, DFDP, JPAR)
609	C DOUBLE PRECISION T, Y, PAR, DFDP
610	C DIMENSION Y(1), PAR(1), DFDP(1)
611	C GO TO (1,2,...,NPAR) JPAR
612	C 1 DFDP(1) = DF(1)/DP(1)
613	C .
614	C DFDP(I) = DF(I)/DP(1)
615	C .
616	C DFDP(N) = DF(N)/DP(1)
617	C RETURN
618	C 2 DFDP(1) = DF(1)/DP(2)
619	C .
620	C DFDP(I) = DF(I)/DP(2)
621	C .
622	C DFDP(N) = DF(N)/DP(2)
623	C .
624	C RETURN
625	C . .
626	C . .
627	C NPAR DFDP(1) = DF(1)/DP(NPAR)
628	C .
629	C DFDP(I) = DF(I)/DP(NPAR)
630	C .
631	C DFDP(N) = DF(N)/DP(NPAR)
632	C RETURN
633	C END
634	C WHERE NEQ, T, Y, PAR, AND JPAR ARE INPUT AND THE VECTOR
635	C DFDP(*,JPAR) IS TO BE LOADED WITH THE PARTIAL DERIVATIVES
636	C DF(Y,T;PAR)/DP(JPAR) ON OUTPUT. ONLY NONZERO ELEMENTS NEED
637	C BE LOADED. T, Y, AND PAR HAVE THE SAME MEANING AS IN
638	C SUBROUTINE F. (IN THE DIMENSION STATEMENT ABOVE, 1 IS A DUMMY
639	C DIMENSION.. IT CAN BE REPLACED BY ANY VALUE).
640	C
641	C DFDP(*,JPAR) IS PRESET TO ZERO BY THE KppSolveR, SO THAT ONLY
642	C THE NONZERO ELEMENTS NEED BE LOADED BY DF. SUBROUTINE DF
643	C MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM IF USED.
644	C IF IDF = 0 (OR ISOPT = 0), A DUMMY ARGUMENT CAN BE USED.
645	C
646	C SUBROUTINE DF MAY ACCESS USER-DEFINED QUANTITIES IN
647	C NEQ(2),... AND PAR(NPAR+1),... IF NEQ IS AN ARRAY
648	C (DIMENSIONED IN DF) AND PAR HAS A LENGTH EXCEEDING NPAR.
649	C SEE THE DESCRIPTIONS OF NEQ AND PAR (BELOW).
650	C
651	C NEQ = THE SIZE OF THE ODE SYSTEM (NUMBER OF FIRST ORDER ORDINARY
652	C DIFFERENTIAL EQUATIONS (N) IN THE MODEL). USED ONLY FOR
653	C INPUT. NEQ MAY NOT BE CHANGED DURING THE PROBLEM.
654	C
655	C FOR ISOPT = 0, NEQ IS NORMALLY A SCALAR. HOWEVER, NEQ MAY
656	C BE AN ARRAY, WITH NEQ(1) SET TO THE SYSTEM SIZE (N), IN WHICH
657	C CASE THE ODESSA PACKAGE ACCESSES ONLY NEQ(1). HOWEVER,
658	C THIS PARAMETER IS PASSED AS THE NEQ ARGUMENT IN ALL CALLS
659	C TO F, DF, AND JAC. HENCE, IF IT IS AN ARRAY, LOCATIONS
660	C NEQ(2),... MAY BE USED TO STORE OTHER INTEGER DATA AND PASS
661	C IT TO F, DF, AND/OR JAC. FOR ISOPT = 1, NPAR MUST BE LOADED
662	C INTO NEQ(2), AND IS NOT ALLOWED TO CHANGE DURING THE PROBLEM.
663	C IN THESE CASES, SUBROUTINES F, DF, AND/OR JAC MUST INCLUDE
664	C NEQ IN A DIMENSION STATEMENT.
665	C
666	C Y = A REAL ARRAY FOR THE VECTOR OF DEPENDENT VARIABLES, OF
667	C DIMENSION (N) BY (NPAR+1). USED FOR BOTH INPUT AND
668	C OUTPUT ON THE FIRST CALL (ISTATE = 1), AND ONLY FOR
669	C OUTPUT ON OTHER CALLS. ON THE FIRST CALL, Y MUST CONTAIN
670	C THE VECTORS OF INITIAL VALUES. ON OUTPUT, Y CONTAINS THE
671	C COMPUTED SOLUTION VECTORS, EVALUATED AT T.
672	C
673	C PAR = A REAL ARRAY FOR THE VECTOR OF CONSTANT MODEL PARAMETERS
674	C OF INTEREST IN THE SENSITIVITY ANALYSIS, OF LENGTH NPAR
675	C OR MORE. PAR IS PASSED AS AN ARGUMENT IN ALL CALLS TO F,
676	C DF, AND JAC. HENCE LOCATIONS PAR(NPAR+1),... MAY BE USED
677	C TO STORE OTHER REAL DATA AND PASS IT TO F, DF, AND/OR JAC.
678	C LOCATIONS PAR(1),...,PAR(NPAR) ARE USED AS INPUT ONLY,
679	C AND MUST NOT BE CHANGED DURING THE PROBLEM.
680	C
681	C T = THE INDEPENDENT VARIABLE. ON INPUT, T IS USED ONLY ON THE
682	C FIRST CALL, AS THE INITIAL POINT OF THE INTEGRATION.
683	C ON OUTPUT, AFTER EACH CALL, T IS THE VALUE AT WHICH A
684	C COMPUTED SOLUTION Y IS EVALUATED (USUALLY THE SAME AS TOUT).
685	C ON AN ERROR RETURN, T IS THE FARTHEST POINT REACHED.
686	C
687	C TOUT = THE NEXT VALUE OF T AT WHICH A COMPUTED SOLUTION IS DESIRED.
688	C USED ONLY FOR INPUT.
689	C
690	C WHEN STARTING THE PROBLEM (ISTATE = 1), TOUT MAY BE EQUAL
691	C TO T FOR ONE CALL, THEN SHOULD .NE. T FOR THE NEXT CALL.
692	C FOR THE INITIAL T, AN INPUT VALUE OF TOUT .NE. T IS USED
693	C IN ORDER TO DETERMINE THE DIRECTION OF THE INTEGRATION
694	C (I.E. THE ALGEBRAIC SIGN OF THE STEP SIZES) AND THE ROUGH
695	C SCALE OF THE PROBLEM. INTEGRATION IN EITHER DIRECTION
696	C (FORWARD OR BACKWARD IN T) IS PERMITTED.
697	C
698	C IF ITASK = 2 OR 5 (ONE-STEP MODES), TOUT IS IGNORED AFTER
699	C THE FIRST CALL (I.E. THE FIRST CALL WITH TOUT .NE. T).
700	C OTHERWISE, TOUT IS REQUIRED ON EVERY CALL.
701	C
702	C IF ITASK = 1, 3, OR 4, THE VALUES OF TOUT NEED NOT BE
703	C MONOTONE, BUT A VALUE OF TOUT WHICH BACKS UP IS LIMITED
704	C TO THE CURRENT INTERNAL T INTERVAL, WHOSE ENDPOINTS ARE
705	C TCUR - HU AND TCUR (SEE OPTIONAL OUTPUTS, BELOW, FOR
706	C TCUR AND HU).
707	C
708	C ITOL = AN INDICATOR FOR THE TYPE OF ERROR CONTROL. SEE
709	C DESCRIPTION BELOW UNDER ATOL. USED ONLY FOR INPUT.
710	C
711	C RTOL = A RELATIVE ERROR TOLERANCE PARAMETER, EITHER A SCALAR OR
712	C AN ARRAY OF SPACE (N) BY (NPAR+1). SEE DESCRIPTION BELOW
713	C UNDER ATOL. INPUT ONLY.
714	C
715	C ATOL = AN ABSOLUTE ERROR TOLERANCE PARAMETER, EITHER A SCALAR OR
716	C AN ARRAY OF SPACE (N) BY (NPAR+1). INPUT ONLY.
717	C
718	C THE INPUT PARAMETERS ITOL, RTOL, AND ATOL DETERMINE
719	C THE ERROR CONTROL PERFORMED BY THE KppSolveR. THE KppSolveR WILL
720	C CONTROL THE VECTOR E = (E(I,J)) OF ESTIMATED LOCAL ERRORS
721	C IN Y, ACCORDING TO AN INEQUALITY OF THE FORM
722	C RMS-NORM OF ( E(I,J)/EWT(I,J) ) .LE. 1,
723	C WHERE EWT(I,J) = RTOL(I,J)*ABS(Y(I,J)) + ATOL(I,J),
724	C AND THE RMS-NORM (ROOT-MEAN-SQUARE NORM) HERE IS
725	C RMS-NORM(V) = SQRT ( (1/N) * SUM (V(I,J)**2) ); I =1,...,N.
726	C HERE EWT = (EWT(I,J)) IS A VECTOR OF WEIGHTS WHICH MUST
727	C ALWAYS BE POSITIVE, AND THE VALUES OF RTOL AND ATOL SHOULD
728	C ALL BE NON-NEGATIVE. THE FOLLOWING TABLE GIVES THE TYPES
729	C (SCALAR/ARRAY) OF RTOL AND ATOL, AND THE CORRESPONDING FORM
730	C OF EWT(I,J).
731	C
732	C ITOL RTOL ATOL EWT(I,J)
733	C 1 SCALAR SCALAR RTOL*ABS(Y(I,J)) + ATOL
734	C 2 SCALAR ARRAY RTOL*ABS(Y(I,J)) + ATOL(I,J)
735	C 3 ARRAY SCALAR RTOL(I,J)*ABS(Y(I,J)) + ATOL
736	C 4 ARRAY ARRAY RTOL(I,J)*ABS(Y(I,J)) + ATOL(I,J)
737	C
738	C WHEN EITHER OF THESE PARAMETERS IS A SCALAR, IT NEED NOT
739	C BE DIMENSIONED IN THE USER-S CALLING PROGRAM.
740	C
741	C THE TOTAL NUMBER OF ERROR TEST FAILURES DUE TO THE SENSITIVITY
742	C ANALYSIS, AND WHICH REQUIRE AN INTEGRATION STEP TO BE
743	C REPEATED, ARE ACCUMULATED IN THE LAST NPAR+1 LOCATIONS OF THE
744	C INTEGER WORK ARRAY IWORK (SEE OPTIONAL OUTPUTS BELOW).
745	C THIS INFORMATION MAY BE OF VALUE IN DETERMINING APPROPRIATE
746	C ERROR TOLERANCES TO BE APPLIED TO THE SENSITIVITY FUNCTIONS.
747	C
748	C IF NONE OF THE ABOVE CHOICES (WITH ITOL, RTOL, AND ATOL
749	C FIXED THROUGHOUT THE PROBLEM) IS SUITABLE, MORE GENERAL
750	C ERROR CONTROLS CAN BE OBTAINED BY SUBSTITUTING
751	C USER-SUPPLIED ROUTINES FOR THE SETTING OF EWT AND/OR FOR
752	C THE NORM CALCULATION. SEE PART IV BELOW.
753	C
754	C IF GLOBAL ERRORS ARE TO BE ESTIMATED BY MAKING A REPEATED
755	C RUN ON THE SAME PROBLEM WITH SMALLER TOLERANCES, THEN ALL
756	C COMPONENTS OF RTOL AND ATOL (I.E. OF EWT) SHOULD BE SCALED
757	C DOWN UNIFORMLY.
758	C
759	C ITASK = AN INDEX SPECIFYING THE TASK TO BE PERFORMED.
760	C INPUT ONLY. ITASK HAS THE FOLLOWING VALUES AND MEANINGS.
761	C 1 MEANS NORMAL COMPUTATION OF OUTPUT VALUES OF Y(T) AT
762	C T = TOUT (BY OVERSHOOTING AND INTERPOLATING).
763	C 2 MEANS TAKE ONE STEP ONLY AND RETURN.
764	C 3 MEANS STOP AT THE FIRST INTERNAL MESH POINT AT OR
765	C BEYOND T = TOUT AND RETURN.
766	C 4 MEANS NORMAL COMPUTATION OF OUTPUT VALUES OF Y(T) AT
767	C T = TOUT BUT WITHOUT OVERSHOOTING T = TCRIT.
768	C TCRIT MUST BE INPUT AS RWORK(1). TCRIT MAY BE EQUAL TO
769	C OR BEYOND TOUT, BUT NOT BEHIND IT IN THE DIRECTION OF
770	C INTEGRATION. THIS OPTION IS USEFUL IF THE PROBLEM
771	C HAS A SINGULARITY AT OR BEYOND T = TCRIT.
772	C 5 MEANS TAKE ONE STEP, WITHOUT PASSING TCRIT, AND RETURN.
773	C TCRIT MUST BE INPUT AS RWORK(1).
774	C
775	C NOTE.. IF ITASK = 4 OR 5 AND THE KppSolveR REACHES TCRIT
776	C (WITHIN ROUNDOFF), IT WILL RETURN T = TCRIT (EXACTLY) TO
777	C INDICATE THIS (UNLESS ITASK = 4 AND TOUT COMES BEFORE TCRIT,
778	C IN WHICH CASE ANSWERS AT T = TOUT ARE RETURNED FIRST).
779	C
780	C ISTATE = AN INDEX USED FOR INPUT AND OUTPUT TO SPECIFY THE
781	C THE STATE OF THE CALCULATION.
782	C
783	C ON INPUT, THE VALUES OF ISTATE ARE AS FOLLOWS.
784	C 1 MEANS THIS IS THE FIRST CALL FOR THE PROBLEM
785	C (INITIALIZATIONS WILL BE DONE). SEE NOTE BELOW.
786	C 2 MEANS THIS IS NOT THE FIRST CALL, AND THE CALCULATION
787	C IS TO CONTINUE NORMALLY, WITH NO CHANGE IN ANY INPUT
788	C PARAMETERS EXCEPT POSSIBLY TOUT AND ITASK.
789	C (IF ITOL, RTOL, AND/OR ATOL ARE CHANGED BETWEEN CALLS
790	C WITH ISTATE = 2, THE NEW VALUES WILL BE USED BUT NOT
791	C TESTED FOR LEGALITY.)
792	C 3 MEANS THIS IS NOT THE FIRST CALL, AND THE
793	C CALCULATION IS TO CONTINUE NORMALLY, BUT WITH
794	C A CHANGE IN INPUT PARAMETERS OTHER THAN
795	C TOUT AND ITASK. CHANGES ARE ALLOWED IN
796	C ITOL, RTOL, ATOL, IOPT, LRW, LIW, MF, ML, MU,
797	C AND ANY OF THE OPTIONAL INPUTS EXCEPT H0.
798	C (SEE IWORK DESCRIPTION FOR ML AND MU.)
799	C NOTE.. A PRELIMINARY CALL WITH TOUT = T IS NOT COUNTED
800	C AS A FIRST CALL HERE, AS NO INITIALIZATION OR CHECKING OF
801	C INPUT IS DONE. (SUCH A CALL IS SOMETIMES USEFUL FOR THE
802	C PURPOSE OF OUTPUTTING THE INITIAL CONDITIONS.)
803	C THUS THE FIRST CALL FOR WHICH TOUT .NE. T REQUIRES
804	C ISTATE = 1 ON INPUT.
805	C
806	C ON OUTPUT, ISTATE HAS THE FOLLOWING VALUES AND MEANINGS.
807	C 1 MEANS NOTHING WAS DONE, AS TOUT WAS EQUAL TO T WITH
808	C ISTATE = 1 ON INPUT. (HOWEVER, AN INTERNAL COUNTER WAS
809	C SET TO DETECT AND PREVENT REPEATED CALLS OF THIS TYPE.)
810	C 2 MEANS THE INTEGRATION WAS PERFORMED SUCCESSFULLY.
811	C -1 MEANS AN EXCESSIVE AMOUNT OF WORK (MORE THAN MXSTEP
812	C STEPS) WAS DONE ON THIS CALL, BEFORE COMPLETING THE
813	C REQUESTED TASK, BUT THE INTEGRATION WAS OTHERWISE
814	C SUCCESSFUL AS FAR AS T. (MXSTEP IS AN OPTIONAL INPUT
815	C AND IS NORMALLY 500.) TO CONTINUE, THE USER MAY
816	C SIMPLY RESET ISTATE TO A VALUE .GT. 1 AND CALL AGAIN
817	C (THE EXCESS WORK STEP COUNTER WILL BE RESET TO 0).
818	C IN ADDITION, THE USER MAY INCREASE MXSTEP TO AVOID
819	C THIS ERROR RETURN (SEE BELOW ON OPTIONAL INPUTS).
820	C -2 MEANS TOO MUCH ACCURACY WAS REQUESTED FOR THE PRECISION
821	C OF THE MACHINE BEING USED. THIS WAS DETECTED BEFORE
822	C COMPLETING THE REQUESTED TASK, BUT THE INTEGRATION
823	C WAS SUCCESSFUL AS FAR AS T. TO CONTINUE, THE TOLERANCE
824	C PARAMETERS MUST BE RESET, AND ISTATE MUST BE SET
825	C TO 3. THE OPTIONAL OUTPUT TOLSF MAY BE USED FOR THIS
826	C PURPOSE. (NOTE.. IF THIS CONDITION IS DETECTED BEFORE
827	C TAKING ANY STEPS, THEN AN ILLEGAL INPUT RETURN
828	C (ISTATE = -3) OCCURS INSTEAD.)
829	C -3 MEANS ILLEGAL INPUT WAS DETECTED, BEFORE TAKING ANY
830	C INTEGRATION STEPS. SEE WRITTEN MESSAGE FOR DETAILS.
831	C NOTE.. IF THE KppSolveR DETECTS AN INFINITE LOOP OF CALLS
832	C TO THE KppSolveR WITH ILLEGAL INPUT, IT WILL CAUSE
833	C THE RUN TO STOP.
834	C -4 MEANS THERE WERE REPEATED ERROR TEST FAILURES ON
835	C ONE ATTEMPTED STEP, BEFORE COMPLETING THE REQUESTED
836	C TASK, BUT THE INTEGRATION WAS SUCCESSFUL AS FAR AS T.
837	C THE PROBLEM MAY HAVE A SINGULARITY, OR THE INPUT
838	C MAY BE INAPPROPRIATE.
839	C -5 MEANS THERE WERE REPEATED CONVERGENCE TEST FAILURES ON
840	C ONE ATTEMPTED STEP, BEFORE COMPLETING THE REQUESTED
841	C TASK, BUT THE INTEGRATION WAS SUCCESSFUL AS FAR AS T.
842	C THIS MAY BE CAUSED BY AN INACCURATE JACOBIAN MATRIX,
843	C IF ONE IS BEING USED.
844	C -6 MEANS EWT(I,J) BECAME ZERO FOR SOME I,J DURING THE
845	C INTEGRATION. PURE RELATIVE ERROR CONTROL (ATOL(I,J)=0.0)
846	C WAS REQUESTED ON A VARIABLE WHICH HAS NOW VANISHED.
847	C THE INTEGRATION WAS SUCCESSFUL AS FAR AS T.
848	C
849	C NOTE.. SINCE THE NORMAL OUTPUT VALUE OF ISTATE IS 2,
850	C IT DOES NOT NEED TO BE RESET FOR NORMAL CONTINUATION.
851	C ALSO, SINCE A NEGATIVE INPUT VALUE OF ISTATE WILL BE
852	C REGARDED AS ILLEGAL, A NEGATIVE OUTPUT VALUE REQUIRES THE
853	C USER TO CHANGE IT, AND POSSIBLY OTHER INPUTS, BEFORE
854	C CALLING THE KppSolveR AGAIN.
855	C
856	C IOPT = AN INTEGER ARRAY FLAG TO SPECIFY WHETHER OR NOT ANY OPTIONAL
857	C INPUTS ARE BEING USED ON THIS CALL. INPUT ONLY.
858	C THE OPTIONAL INPUTS ARE LISTED SEPARATELY BELOW.
859	C IOPT(1) = 0 MEANS NO OPTIONAL INPUTS FOR THE KppSolveR WILL BE
860	C USED. DEFAULT VALUES WILL BE USED IN ALL CASES.
861	C = 1 MEANS ONE OR MORE OPTIONAL INPUTS FOR THE
862	C KppSolveR ARE BEING USED.
863	C NOTE : IOPT(1) IS INDEPENDENT OF ISOPT AND IDF.
864	C IOPT(2) = 0 MEANS NO SENSITIVITY ANALYSIS WILL BE PERFORMED.
865	C = 1 MEANS A SENSITIVITY ANALYSIS WILL BE PERFORMED.
866	C NOTE : IOPT(2) IS RENAMED TO ISOPT IN ODESSA.
867	C = 0 MEANS DF/DP WILL BE CALCULATED BY FINITE
868	C DIFFERENCE WITHIN ODESSA.
869	C IOPT(3) = 1 MEANS DF/DP WILL BE CALCULATED BY A USER-SUPPLIED
870	C ROUTINE.
871	C NOTE : IOPT(3) IS RENAMED TO IDF IN ODESSA.
872	C IF IDF = 1, THE USER MUST SUPPLY A
873	C SUBROUTINE DF (THE NAME IS ARBITRARY) AS
874	C DESCRIBED BELOW UNDER DF. FOR IDF = 0,
875	C A DUMMY ARGUMENT CAN BE USED.
876	C
877	C RWORK = A REAL WORKING ARRAY (DOUBLE PRECISION).
878	C FOR ISOPT = 0, THE LENGTH OF RWORK MUST BE AT LEAST..
879	C 20 + NYH(MAXORD + 1) + 3NEQ + LWM
880	C FOR ISOPT = 1, THE LENGTH OF RWORK MUST BE AT LEAST..
881	C 20 + NYH(MAXORD + 1) + 2NYH + LWM + N
882	C WHERE..
883	C NYH = THE TOTAL NUMBER OF DEPENDENT VARIABLES;
884	C (= N IF ISOPT = 0, AND N*(NPAR+1) IF ISOPT = 1).
885	C MAXORD = 12 (IF METH = 1) OR 5 (IF METH = 2) (UNLESS A
886	C SMALLER VALUE IS GIVEN AS AN OPTIONAL INPUT),
887	C LWM = 0 IF MITER = 0,
888	C LWM = N**2 + 2 IF MITER IS 1 OR 2,
889	C LWM = N + 2 IF MITER = 3, AND
890	C LWM = (2ML+MU+1)N + 2 IF MITER IS 4 OR 5.
891	C (SEE THE MF DESCRIPTION FOR METH AND MITER.)
892	C
893	C THE FIRST 20 WORDS OF RWORK ARE RESERVED FOR CONDITIONAL
894	C AND OPTIONAL INPUTS AND OPTIONAL OUTPUTS.
895	C
896	C THE FOLLOWING WORD IN RWORK IS A CONDITIONAL INPUT..
897	C RWORK(1) = TCRIT = CRITICAL VALUE OF T WHICH THE KppSolveR
898	C IS NOT TO OVERSHOOT. REQUIRED IF ITASK IS
899	C 4 OR 5, AND IGNORED OTHERWISE. (SEE ITASK.)
900	C
901	C LRW = THE LENGTH OF THE ARRAY RWORK, AS DECLARED BY THE USER.
902	C (THIS WILL BE CHECKED BY THE KppSolveR.)
903	C
904	C IWORK = AN INTEGER WORK ARRAY. THE LENGTH MUST BE AT LEAST..
905	C 20 IF MITER = 0 OR 3 (MF = 10, 13, 20, 23), OR
906	C 20 + N OTHERWISE (MF = 11, 12, 14, 15, 21, 22, 24, 25).
907	C FOR ISOPT = 0, OR..
908	C 21 + N + NPAR
909	C FOR ISOPT = 1.
910	C THE FIRST FEW WORDS OF IWORK ARE USED FOR CONDITIONAL AND
911	C OPTIONAL INPUTS AND OPTIONAL OUTPUTS.
912	C
913	C THE FOLLOWING 2 WORDS IN IWORK ARE CONDITIONAL INPUTS..
914	C IWORK(1) = ML THESE ARE THE LOWER AND UPPER
915	C IWORK(2) = MU HALF-BANDWIDTHS, RESPECTIVELY, OF THE
916	C BANDED JACOBIAN, EXCLUDING THE MAIN DIAGONAL.
917	C THE BAND IS DEFINED BY THE MATRIX LOCATIONS
918	C (I,J) WITH I-ML .LE. J .LE. I+MU. ML AND MU
919	C MUST SATISFY 0 .LE. ML,MU .LE. NEQ-1.
920	C THESE ARE REQUIRED IF MITER IS 4 OR 5, AND
921	C IGNORED OTHERWISE. ML AND MU MAY IN FACT BE
922	C THE BAND PARAMETERS FOR A MATRIX TO WHICH
923	C DF/DY IS ONLY APPROXIMATELY EQUAL.
924	*
925	C
926	C LIW = THE LENGTH OF THE ARRAY IWORK, AS DECLARED BY THE USER.
927	C (THIS WILL BE CHECKED BY THE KppSolveR.)
928	C
929	C NOTE.. THE WORK ARRAYS MUST NOT BE ALTERED BETWEEN CALLS TO ODESSA
930	C FOR THE SAME PROBLEM, EXCEPT POSSIBLY FOR THE CONDITIONAL AND
931	C OPTIONAL INPUTS, AND EXCEPT FOR THE LAST 2*NYH + N WORDS OF RWORK.
932	C THE LATTER SPACE IS USED FOR INTERNAL SCRATCH SPACE, AND SO IS
933	C AVAILABLE FOR USE BY THE USER OUTSIDE ODESSA BETWEEN CALLS, IF
934	C DESIRED (BUT NOT FOR USE BY F, DF, OR JAC).
935	C
936	C JAC = THE NAME OF THE USER-SUPPLIED ROUTINE (MITER = 1 OR 4) TO
937	C COMPUTE THE JACOBIAN MATRIX, DF/DY, AS A FUNCTION OF THE
938	C SCALAR T AND THE VECTORS Y, AND PAR. IT IS TO HAVE THE FORM
939	C SUBROUTINE JAC (NEQ, T, Y, PAR, ML, MU, PD, NROWPD)
940	C DOUBLE PRECISION T, Y, PAR, PD
941	C DIMENSION Y(1), PAR(1), PD(NROWPD,1)
942	C WHERE NEQ, T, Y, PAR, ML, MU, AND NROWPD ARE INPUT AND THE
943	C ARRAY PD IS TO BE LOADED WITH PARTIAL DERIVATIVES (ELEMENTS
944	C OF THE JACOBIAN MATRIX) ON OUTPUT. PD MUST BE GIVEN A FIRST
945	C DIMENSION OF NROWPD. T, Y, AND PAR HAVE THE SAME MEANING AS
946	C IN SUBROUTINE F. (IN THE DIMENSION STATEMENT ABOVE, 1 IS A
947	C DUMMY DIMENSION.. IT CAN BE REPLACED BY ANY VALUE.)
948	C IN THE FULL MATRIX CASE (MITER = 1), ML AND MU ARE
949	C IGNORED, AND THE JACOBIAN IS TO BE LOADED INTO PD IN
950	C COLUMNWISE MANNER, WITH DF(I)/DY(J) LOADED INTO PD(I,J).
951	C IN THE BAND MATRIX CASE (MITER = 4), THE ELEMENTS
952	C WITHIN THE BAND ARE TO BE LOADED INTO PD IN COLUMNWISE
953	C MANNER, WITH DIAGONAL LINES OF DF/DY LOADED INTO THE ROWS
954	C OF PD. THUS DF(I)/DY(J) IS TO BE LOADED INTO PD(I-J+MU+1,J).
955	C ML AND MU ARE THE HALF-BANDWIDTH PARAMETERS (SEE IWORK).
956	C THE LOCATIONS IN PD IN THE TWO TRIANGULAR AREAS WHICH
957	C CORRESPOND TO NONEXISTENT MATRIX ELEMENTS CAN BE IGNORED
958	C OR LOADED ARBITRARILY, AS THEY ARE OVERWRITTEN BY ODESSA.
959	C PD IS PRESET TO ZERO BY THE KppSolveR, SO THAT ONLY THE
960	C NONZERO ELEMENTS NEED BE LOADED BY JAC. EACH CALL TO JAC IS
961	C PRECEDED BY A CALL TO F WITH THE SAME ARGUMENTS NEQ, T, Y,
962	C AND PAR. THUS TO GAIN SOME EFFICIENCY, INTERMEDIATE
963	C QUANTITIES SHARED BY BOTH CALCULATIONS MAY BE SAVED IN A
964	C USER COMMON BLOCK BY F AND NOT RECOMPUTED BY JAC, IF
965	C DESIRED. ALSO, JAC MAY ALTER THE Y ARRAY, IF DESIRED.
966	C JAC MUST BE DECLARED EXTERNAL IN THE CALLING PROGRAM.
967	C SUBROUTINE JAC MAY ACCESS USER-DEFINED QUANTITIES IN
968	C NEQ(2),... AND PAR(NPAR+1),.... SEE THE DESCRIPTIONS OF
969	C NEQ (ABOVE) AND PAR (BELOW).
970	C
971	C MF = THE METHOD FLAG. USED ONLY FOR INPUT. THE LEGAL VALUES OF
972	C MF ARE 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24, AND 25.
973	C MF HAS DECIMAL DIGITS METH AND MITER.. MF = 10*METH + MITER.
974	C METH INDICATES THE BASIC LINEAR MULTISTEP METHOD..
975	C METH = 1 MEANS THE IMPLICIT ADAMS METHOD.
976	*
977	C METH = 2 MEANS THE METHOD BASED ON BACKWARD
978	C DIFFERENTIATION FORMULAS (BDF-S).
979	C MITER INDICATES THE CORRECTOR ITERATION METHOD..
980	C MITER = 0 MEANS FUNCTIONAL ITERATION (NO JACOBIAN MATRIX
981	C IS INVOLVED).
982	C MITER = 1 MEANS CHORD ITERATION WITH A USER-SUPPLIED
983	C FULL (NEQ BY NEQ) JACOBIAN.
984	C MITER = 2 MEANS CHORD ITERATION WITH AN INTERNALLY
985	C GENERATED (DIFFERENCE QUOTIENT) FULL JACOBIAN
986	C (USING NEQ EXTRA CALLS TO F PER DF/DY VALUE).
987	C MITER = 3 MEANS CHORD ITERATION WITH AN INTERNALLY
988	C GENERATED DIAGONAL JACOBIAN APPROXIMATION.
989	C (USING 1 EXTRA CALL TO F PER DF/DY EVALUATION).
990	C MITER = 4 MEANS CHORD ITERATION WITH A USER-SUPPLIED
991	C BANDED JACOBIAN.
992	C MITER = 5 MEANS CHORD ITERATION WITH AN INTERNALLY
993	C GENERATED BANDED JACOBIAN (USING ML+MU+1 EXTRA
994	C CALLS TO F PER DF/DY EVALUATION).
995	C IF MITER = 1 OR 4, THE USER MUST SUPPLY A SUBROUTINE JAC
996	C (THE NAME IS ARBITRARY) AS DESCRIBED ABOVE UNDER JAC.
997	C FOR OTHER VALUES OF MITER, A DUMMY ARGUMENT CAN BE USED.
998	C
999	C IF A SENSITIVITY ANLYSIS IS DESIRED (ISOPT = 1), MITER = 0
1000	C AND 3 ARE DISALLOWED. IN THESE CASES, THE USER IS RECOMMENDED
1001	C TO SUPPLY AN ANALYTICAL JACOBIAN (MITER = 1 OR 4) AND AN
1002	C ANALYTICAL INHOMOGENEITY MATRIX (IDF = 1).
1003	C----------------------------------------------------------------------
1004	C OPTIONAL INPUTS.
1005	C
1006	C THE FOLLOWING IS A LIST OF THE OPTIONAL INPUTS PROVIDED FOR IN THE
1007	C CALL SEQUENCE. (SEE ALSO PART II.) FOR EACH SUCH INPUT VARIABLE,
1008	C THIS TABLE LISTS ITS NAME AS USED IN THIS DOCUMENTATION, ITS
1009	C LOCATION IN THE CALL SEQUENCE, ITS MEANING, AND THE DEFAULT VALUE.
1010	C THE USE OF ANY OF THESE INPUTS REQUIRES IOPT(1) = 1, AND IN THAT
1011	C CASE ALL OF THESE INPUTS ARE EXAMINED. A VALUE OF ZERO FOR ANY
1012	C OF THESE OPTIONAL INPUTS WILL CAUSE THE DEFAULT VALUE TO BE USED.
1013	C THUS TO USE A SUBSET OF THE OPTIONAL INPUTS, SIMPLY PRELOAD
1014	C LOCATIONS 5 TO 10 IN RWORK AND IWORK TO 0.0 AND 0 RESPECTIVELY, AND
1015	C THEN SET THOSE OF INTEREST TO NONZERO VALUES.
1016	C
1017	C NAME LOCATION MEANING AND DEFAULT VALUE
1018	C
1019	C H0 RWORK(5) THE STEP SIZE TO BE ATTEMPTED ON THE FIRST STEP.
1020	C THE DEFAULT VALUE IS DETERMINED BY THE KppSolveR.
1021	C
1022	C HMAX RWORK(6) THE MAXIMUM ABSOLUTE STEP SIZE ALLOWED.
1023	C THE DEFAULT VALUE IS INFINITE.
1024	C
1025	C HMIN RWORK(7) THE MINIMUM ABSOLUTE STEP SIZE ALLOWED.
1026	C THE DEFAULT VALUE IS 0. (THIS LOWER BOUND IS NOT
1027	C ENFORCED ON THE FINAL STEP BEFORE REACHING TCRIT
1028	C WHEN ITASK = 4 OR 5.)
1029	C
1030	C MAXORD IWORK(5) THE MAXIMUM ORDER TO BE ALLOWED. THE DEFAULT
1031	C VALUE IS 12 IF METH = 1, AND 5 IF METH = 2.
1032	C IF MAXORD EXCEEDS THE DEFAULT VALUE, IT WILL
1033	C BE REDUCED TO THE DEFAULT VALUE.
1034	C IF MAXORD IS CHANGED DURING THE PROBLEM, IT MAY
1035	C CAUSE THE CURRENT ORDER TO BE REDUCED.
1036	C
1037	C MXSTEP IWORK(6) MAXIMUM NUMBER OF (INTERNALLY DEFINED) STEPS
1038	C ALLOWED DURING ONE CALL TO THE KppSolveR.
1039	C THE DEFAULT VALUE IS 500.
1040	C
1041	C MXHNIL IWORK(7) MAXIMUM NUMBER OF MESSAGES PRINTED (PER PROBLEM)
1042	C WARNING THAT T + H = T ON A STEP (H = STEP SIZE).
1043	C THIS MUST BE POSITIVE TO RESULT IN A NON-DEFAULT
1044	C VALUE. THE DEFAULT VALUE IS 10.
1045	C----------------------------------------------------------------------
1046	C OPTIONAL OUTPUTS.
1047	C
1048	C AS OPTIONAL ADDITIONAL OUTPUT FROM ODESSA, THE VARIABLES LISTED
1049	C BELOW ARE QUANTITIES RELATED TO THE PERFORMANCE OF ODESSA
1050	C WHICH ARE AVAILABLE TO THE USER. THESE ARE COMMUNICATED BY WAY OF
1051	C THE WORK ARRAYS, BUT ALSO HAVE INTERNAL MNEMONIC NAMES AS SHOWN.
1052	C EXCEPT WHERE STATED OTHERWISE, ALL OF THESE OUTPUTS ARE DEFINED
1053	C ON ANY SUCCESSFUL RETURN FROM ODESSA, AND ON ANY RETURN WITH
1054	C ISTATE = -1, -2, -4, -5, OR -6. ON AN ILLEGAL INPUT RETURN
1055	C (ISTATE = -3), THEY WILL BE UNCHANGED FROM THEIR EXISTING VALUES
1056	C (IF ANY), EXCEPT POSSIBLY FOR TOLSF, LENRW, AND LENIW.
1057	C ON ANY ERROR RETURN, OUTPUTS RELEVANT TO THE ERROR WILL BE DEFINED,
1058	C AS NOTED BELOW.
1059	C
1060	C NAME LOCATION MEANING
1061	C
1062	C HU RWORK(11) THE STEP SIZE IN T LAST USED (SUCCESSFULLY).
1063	C
1064	C HCUR RWORK(12) THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP.
1065	C
1066	C TCUR RWORK(13) THE CURRENT VALUE OF THE INDEPENDENT VARIABLE
1067	C WHICH THE KppSolveR HAS ACTUALLY REACHED, I.E. THE
1068	C CURRENT INTERNAL MESH POINT IN T. ON OUTPUT, TCUR
1069	C WILL ALWAYS BE AT LEAST AS FAR AS THE ARGUMENT
1070	C T, BUT MAY BE FARTHER (IF INTERPOLATION WAS DONE).
1071	C
1072	C TOLSF RWORK(14) A TOLERANCE SCALE FACTOR, GREATER THAN 1.0,
1073	C COMPUTED WHEN A REQUEST FOR TOO MUCH ACCURACY WAS
1074	C DETECTED (ISTATE = -3 IF DETECTED AT THE START OF
1075	C THE PROBLEM, ISTATE = -2 OTHERWISE). IF ITOL IS
1076	C LEFT UNALTERED BUT RTOL AND ATOL ARE UNIFORMLY
1077	C SCALED UP BY A FACTOR OF TOLSF FOR THE NEXT CALL,
1078	C THEN THE KppSolveR IS DEEMED LIKELY TO SUCCEED.
1079	C (THE USER MAY ALSO IGNORE TOLSF AND ALTER THE
1080	C TOLERANCE PARAMETERS IN ANY OTHER WAY APPROPRIATE.)
1081	C
1082	C NST IWORK(11) THE NUMBER OF STEPS TAKEN FOR THE PROBLEM SO FAR.
1083	C
1084	C NFE IWORK(12) THE NUMBER OF F EVALUATIONS FOR THE PROBLEM SO FAR.
1085	C
1086	C NJE IWORK(13) THE NUMBER OF JACOBIAN EVALUATIONS (AND OF MATRIX
1087	C LU DECOMPOSITIONS IF ISOPT = 0) FOR THE PROBLEM SO
1088	C FAR. IF ISOPT = 1, THE NUMBER OF LU DECOMPOSITIONS
1089	C IS EQUAL TO NJE - NSPE (SEE BELOW).
1090	C
1091	C NQU IWORK(14) THE METHOD ORDER LAST USED (SUCCESSFULLY).
1092	C
1093	C NQCUR IWORK(15) THE ORDER TO BE ATTEMPTED ON THE NEXT STEP.
1094	C
1095	C IMXER IWORK(16) THE INDEX OF THE COMPONENT OF LARGEST MAGNITUDE IN
1096	C THE WEIGHTED LOCAL ERROR VECTOR (E(I,J)/EWT(I,J)),
1097	C ON AN ERROR RETURN WITH ISTATE = -4 OR -5.
1098	C
1099	C LENRW IWORK(17) THE LENGTH OF RWORK ACTUALLY REQUIRED.
1100	C THIS IS DEFINED ON NORMAL RETURNS AND ON AN ILLEGAL
1101	C INPUT RETURN FOR INSUFFICIENT STORAGE.
1102	C
1103	C LENIW IWORK(18) THE LENGTH OF IWORK ACTUALLY REQUIRED.
1104	C THIS IS DEFINED ON NORMAL RETURNS AND ON AN ILLEGAL
1105	C INPUT RETURN FOR INSUFFICIENT STORAGE.
1106	C
1107	C NDFE IWORK(19) THE NUMBER OF DF/DP (VECTOR) EVALUATIONS.
1108	C
1109	C NSPE IWORK(20) THE NUMBER OF CALLS TO SUBROUTINE SPRIME. EACH CALL
1110	C TO SPRIME REQUIRES A JACOBIAN EVALUATION, BUT NOT
1111	C AN LU DECOMPOSITION.
1112	C
1113	C THE FOLLOWING ARRAYS ARE SEGMENTS OF THE RWORK AND IWORK ARRAYS
1114	C WHICH MAY ALSO BE OF INTEREST TO THE USER AS OPTIONAL OUTPUTS.
1115	C FOR EACH ARRAY, THE TABLE BELOW GIVES ITS INTERNAL NAME, ITS BASE
1116	C ADDRESS IN RWORK OR IWORK, AND ITS DESCRIPTION.
1117	C
1118	C NAME BASE ADDRESS DESCRIPTION
1119	C
1120	C YH 21 IN RWORK THE NORDSIECK HISTORY ARRAY, OF SIZE NYH BY
1121	C (NQCUR + 1). FOR J = 0,1,...,NQCUR, COLUMN J+1
1122	C OF YH CONTAINS HCUR**J/FACTORIAL(J) TIMES
1123	C THE J-TH DERIVATIVE OF THE INTERPOLATING
1124	C POLYNOMIAL CURRENTLY REPRESENTING THE SOLUTION,
1125	C EVALUATED AT T = TCUR.
1126	C
1127	C ACOR LENRW-NYH+1 ARRAY OF SIZE NYH USED FOR THE ACCUMULATED
1128	C IN RWORK CORRECTIONS ON EACH STEP, SCALED ON OUTPUT
1129	C TO REPRESENT THE ESTIMATED LOCAL ERROR IN Y
1130	C ON THE LAST STEP. THIS IS THE VECTOR E IN
1131	C THE DESCRIPTION OF THE ERROR CONTROL.
1132	C IT IS DEFINED ONLY ON A SUCCESSFUL RETURN
1133	C FROM ODESSA.
1134	C NRS LENIW-NPAR ARRAY OF SIZE NPAR+1, USED TO STORE THE
1135	C IN IWORK ACCUMULATED NUMBER OF REPEATED STEPS DUE TO
1136	C THE SENSITIVITY ANALYSIS..
1137	C NRS(1) = TOTAL NUMBER OF REPEATED STEPS,
1138	C NRS(2),... = NUMBER OF REPEATED STEPS DUE TO
1139	C MODEL PARAMETER 1,...
1140	C
1141	C----------------------------------------------------------------------
1142	C PART II. OTHER ROUTINES CALLABLE.
1143	C
1144	C THE FOLLOWING ARE OPTIONAL CALLS WHICH THE USER MAY MAKE TO
1145	C GAIN ADDITIONAL CAPABILITIES IN CONJUNCTION WITH ODESSA.
1146	C (THE ROUTINES XSETUN AND XSETF ARE DESIGNED TO CONFORM TO THE
1147	C SLATEC ERROR HANDLING PACKAGE.)
1148	C
1149	C FORM OF CALL FUNCTION
1150	C CALL XSETUN(LUN) SET THE LOGICAL UNIT NUMBER, LUN, FOR
1151	C OUTPUT OF MESSAGES FROM ODESSA, IF
1152	C THE DEFAULT IS NOT DESIRED.
1153	C THE DEFAULT VALUE OF LUN IS 6.
1154	C
1155	C CALL XSETF(MFLAG) SET A FLAG TO CONTROL THE PRINTING OF
1156	C MESSAGES BY ODESSA..
1157	C MFLAG = 0 MEANS DO NOT PRINT. (DANGER..
1158	C THIS RISKS LOSING VALUABLE INFORMATION.)
1159	C MFLAG = 1 MEANS PRINT (THE DEFAULT).
1160	C
1161	C EITHER OF THE ABOVE CALLS MAY BE MADE AT
1162	C ANY TIME AND WILL TAKE EFFECT IMMEDIATELY.
1163	C
1164	C CALL SVCOM (RSAV, ISAV) STORE IN RSAV AND ISAV THE CONTENTS
1165	C OF THE INTERNAL COMMON BLOCKS USED BY
1166	C ODESSA (SEE PART III BELOW).
1167	C RSAV MUST BE A REAL ARRAY OF LENGTH 222
1168	C OR MORE, AND ISAV MUST BE AN INTEGER
1169	C ARRAY OF LENGTH 54 OR MORE.
1170	C
1171	C CALL RSCOM (RSAV, ISAV) RESTORE, FROM RSAV AND ISAV, THE CONTENTS
1172	C OF THE INTERNAL COMMON BLOCKS USED BY
1173	C ODESSA. PRESUMES A PRIOR CALL TO SVCOM
1174	C WITH THE SAME ARGUMENTS.
1175	C
1176	C SVCOM AND RSCOM ARE USEFUL IF
1177	C INTERRUPTING A RUN AND RESTARTING
1178	C LATER, OR ALTERNATING BETWEEN TWO OR
1179	C MORE PROBLEMS KppSolveD WITH ODESSA.
1180	C
1181	C CALL INTDY(,,,,,) PROVIDE DERIVATIVES OF Y, OF VARIOUS
1182	C (SEE BELOW) ORDERS, AT A SPECIFIED POINT T, IF
1183	C DESIRED. IT MAY BE CALLED ONLY AFTER
1184	C A SUCCESSFUL RETURN FROM ODESSA.
1185	C
1186	C THE DETAILED INSTRUCTIONS FOR USING INTDY ARE AS FOLLOWS.
1187	C THE FORM OF THE CALL IS..
1188	C
1189	C CALL INTDY (T, K, RWORK(21), NYH, DKY, IFLAG)
1190	C
1191	C THE INPUT PARAMETERS ARE..
1192	C
1193	C T = VALUE OF INDEPENDENT VARIABLE WHERE ANSWERS ARE DESIRED
1194	C (NORMALLY THE SAME AS THE T LAST RETURNED BY ODESSA).
1195	C FOR VALID RESULTS, T MUST LIE BETWEEN TCUR - HU AND TCUR.
1196	C (SEE OPTIONAL OUTPUTS FOR TCUR AND HU.)
1197	C K = INTEGER ORDER OF THE DERIVATIVE DESIRED. K MUST SATISFY
1198	C 0 .LE. K .LE. NQCUR, WHERE NQCUR IS THE CURRENT ORDER
1199	C (SEE OPTIONAL OUTPUTS). THE CAPABILITY CORRESPONDING
1200	C TO K = 0, I.E. COMPUTING Y(T), IS ALREADY PROVIDED
1201	C BY ODESSA DIRECTLY. SINCE NQCUR .GE. 1, THE FIRST
1202	C DERIVATIVE DY/DT IS ALWAYS AVAILABLE WITH INTDY.
1203	C RWORK(21) = THE BASE ADDRESS OF THE HISTORY ARRAY YH.
1204	C NYH = COLUMN LENGTH OF YH, EQUAL TO THE TOTAL NUMBER OF
1205	C DEPENDENT VARIABLES. IF ISOPT = 0, NYH = N. IF ISOPT = 1,
1206	C NYH = N * (NPAR + 1).
1207	C
1208	C THE OUTPUT PARAMETERS ARE..
1209	C
1210	C DKY = A REAL ARRAY OF LENGTH NYH CONTAINING THE COMPUTED VALUE
1211	C OF THE K-TH DERIVATIVE OF Y(T).
1212	C IFLAG = INTEGER FLAG, RETURNED AS 0 IF K AND T WERE LEGAL,
1213	C -1 IF K WAS ILLEGAL, AND -2 IF T WAS ILLEGAL.
1214	C ON AN ERROR RETURN, A MESSAGE IS ALSO WRITTEN.
1215	C----------------------------------------------------------------------
1216	C PART III. COMMON BLOCKS.
1217	C
1218	C IF ODESSA IS TO BE USED IN AN OVERLAY SITUATION, THE USER
1219	C MUST DECLARE, IN THE PRIMARY OVERLAY, THE VARIABLES IN..
1220	C (1) THE CALL SEQUENCE TO ODESSA,
1221	C (2) THE THREE INTERNAL COMMON BLOCKS
1222	C /ODE001/ OF LENGTH 258 (219 DOUBLE PRECISION WORDS
1223	C FOLLOWED BY 39 INTEGER WORDS),
1224	C /ODE002/ OF LENGTH 14 (3 DOUBLE PRECISION WORDS FOLLOWED
1225	C BY 11 INTEGER WORDS),
1226	C /EH0001/ OF LENGTH 2 (INTEGER WORDS).
1227	C
1228	C IF ODESSA IS USED ON A SYSTEM IN WHICH THE CONTENTS OF INTERNAL
1229	C COMMON BLOCKS ARE NOT PRESERVED BETWEEN CALLS, THE USER SHOULD
1230	C DECLARE THE ABOVE THREE COMMON BLOCKS IN HIS MAIN PROGRAM TO INSURE
1231	C THAT THEIR CONTENTS ARE PRESERVED.
1232	C
1233	C IF THE SOLUTION OF A GIVEN PROBLEM BY ODESSA IS TO BE INTERRUPTED
1234	C AND THEN LATER CONTINUED, SUCH AS WHEN RESTARTING AN INTERRUPTED RUN
1235	C OR ALTERNATING BETWEEN TWO OR MORE PROBLEMS, THE USER SHOULD SAVE,
1236	C FOLLOWING THE RETURN FROM THE LAST ODESSA CALL PRIOR TO THE
1237	C INTERRUPTION, THE CONTENTS OF THE CALL SEQUENCE VARIABLES AND THE
1238	C INTERNAL COMMON BLOCKS, AND LATER RESTORE THESE VALUES BEFORE THE
1239	C NEXT ODESSA CALL FOR THAT PROBLEM. TO SAVE AND RESTORE THE COMMON
1240	C BLOCKS, USE SUBROUTINES SVCOM AND RSCOM (SEE PART II ABOVE).
1241	C
1242	C----------------------------------------------------------------------
1243	C PART IV. OPTIONALLY REPLACEABLE KppSolveR ROUTINES.
1244	C
1245	C BELOW ARE DESCRIPTIONS OF TWO ROUTINES IN THE ODESSA PACKAGE WHICH
1246	C RELATE TO THE MEASUREMENT OF ERRORS. EITHER ROUTINE CAN BE
1247	C REPLACED BY A USER-SUPPLIED VERSION, IF DESIRED. HOWEVER, SINCE SUCH
1248	C A REPLACEMENT MAY HAVE A MAJOR IMPACT ON PERFORMANCE, IT SHOULD BE
1249	C DONE ONLY WHEN ABSOLUTELY NECESSARY, AND ONLY WITH GREAT CAUTION.
1250	C (NOTE.. THE MEANS BY WHICH THE PACKAGE VERSION OF A ROUTINE IS
1251	C SUPERSEDED BY THE USER-S VERSION MAY BE SYSTEM-DEPENDENT.)
1252	C
1253	C (A) EWSET.
1254	C THE FOLLOWING SUBROUTINE IS CALLED JUST BEFORE EACH INTERNAL
1255	C INTEGRATION STEP, AND SETS THE ARRAY OF ERROR WEIGHTS, EWT, AS
1256	C DESCRIBED UNDER ITOL/RTOL/ATOL ABOVE..
1257	C SUBROUTINE EWSET (NYH, ITOL, RTOL, ATOL, YCUR, EWT)
1258	C WHERE NEQ, ITOL, RTOL, AND ATOL ARE AS IN THE ODESSA CALL SEQUENCE,
1259	C YCUR CONTAINS THE CURRENT DEPENDENT VARIABLE VECTOR, AND
1260	C EWT IS THE ARRAY OF WEIGHTS SET BY EWSET.
1261	C
1262	C IF THE USER SUPPLIES THIS SUBROUTINE, IT MUST RETURN IN EWT(I)
1263	C (I = 1,...,NYH) A POSITIVE QUANTITY SUITABLE FOR COMPARING ERRORS
1264	C IN Y(I) TO. THE EWT ARRAY RETURNED BY EWSET IS PASSED TO THE
1265	C VNORM ROUTINE (SEE BELOW), AND ALSO USED BY ODESSA IN THE COMPUTATION
1266	C OF THE OPTIONAL OUTPUT IMXER, THE DIAGONAL JACOBIAN APPROXIMATION,
1267	C AND THE INCREMENTS FOR DIFFERENCE QUOTIENT JACOBIANS.
1268	C
1269	C IN THE USER-SUPPLIED VERSION OF EWSET, IT MAY BE DESIRABLE TO USE
1270	C THE CURRENT VALUES OF DERIVATIVES OF Y. DERIVATIVES UP TO ORDER NQ
1271	C ARE AVAILABLE FROM THE HISTORY ARRAY YH, DESCRIBED ABOVE UNDER
1272	C OPTIONAL OUTPUTS. IN EWSET, YH IS IDENTICAL TO THE YCUR ARRAY,
1273	C EXTENDED TO NQ + 1 COLUMNS WITH A COLUMN LENGTH OF NYH AND SCALE
1274	C FACTORS OF H**J/FACTORIAL(J). ON THE FIRST CALL FOR THE PROBLEM,
1275	C GIVEN BY NST = 0, NQ IS 1 AND H IS TEMPORARILY SET TO 1.0.
1276	C THE QUANTITIES NQ, NYH, H, AND NST CAN BE OBTAINED BY INCLUDING
1277	C IN EWSET THE STATEMENTS..
1278	C DOUBLE PRECISION H, RLS
1279	C COMMON /ODE001/ RLS(219),ILS(39)
1280	C NQ = ILS(35)
1281	C NYH = ILS(14)
1282	C NST = ILS(36)
1283	C H = RLS(213)
1284	C THUS, FOR EXAMPLE, THE CURRENT VALUE OF DY/DT CAN BE OBTAINED AS
1285	C YCUR(NYH+I)/H (I=1,...,N) (AND THE DIVISION BY H IS
1286	C UNNECESSARY WHEN NST = 0).
1287	C
1288	C (B) VNORM.
1289	C THE FOLLOWING IS A REAL FUNCTION ROUTINE WHICH COMPUTES THE WEIGHTED
1290	C ROOT-MEAN-SQUARE NORM OF A VECTOR V..
1291	C D = VNORM (LV, V, W)
1292	C WHERE..
1293	C LV = THE LENGTH OF THE VECTOR,
1294	C V = REAL ARRAY OF LENGTH N CONTAINING THE VECTOR,
1295	C W = REAL ARRAY OF LENGTH N CONTAINING WEIGHTS,
1296	C D = SQRT( (1/N) * SUM(V(I)W(I))*2 ).
1297	C VNORM IS CALLED WITH LV = N AND WITH W(I) = 1.0/EWT(I), WHERE
1298	C EWT IS AS SET BY SUBROUTINE EWSET.
1299	C
1300	C IF THE USER SUPPLIES THIS FUNCTION, IT SHOULD RETURN A NON-NEGATIVE
1301	C VALUE OF VNORM SUITABLE FOR USE IN THE ERROR CONTROL IN ODESSA.
1302	C NONE OF THE ARGUMENTS SHOULD BE ALTERED BY VNORM.
1303	C FOR EXAMPLE, A USER-SUPPLIED VNORM ROUTINE MIGHT..
1304	C -SUBSTITUTE A MAX-NORM OF (V(I)*W(I)) FOR THE RMS-NORM, OR
1305	C -IGNORE SOME COMPONENTS OF V IN THE NORM, WITH THE EFFECT OF
1306	C SUPPRESSING THE ERROR CONTROL ON THOSE COMPONENTS OF Y.
1307	C----------------------------------------------------------------------
1308	C OTHER ROUTINES IN THE ODESSA PACKAGE.
1309	C
1310	C IN ADDITION TO SUBROUTINE ODESSA, THE ODESSA PACKAGE INCLUDES THE
1311	C FOLLOWING SUBROUTINES AND FUNCTION ROUTINES..
1312	C INTDY COMPUTES AN INTERPOLATED VALUE OF THE Y VECTOR AT T = TOUT.
1313	C STODE IS THE CORE INTEGRATOR, WHICH DOES ONE STEP OF THE
1314	C INTEGRATION AND THE ASSOCIATED ERROR CONTROL.
1315	C STESA MANAGES THE SOLUTION OF THE SENSITIVITY FUNCTIONS.
1316	C CFODE SETS ALL METHOD COEFFICIENTS AND TEST CONSTANTS.
1317	C PREPJ COMPUTES AND PREPROCESSES THE JACOBIAN MATRIX J = DF/DY
1318	C AND THE NEWTON ITERATION MATRIX P = I - HL0J.
1319	C IT IS ALSO CALLED BY SPRIME (WITH JOPT = 1) TO JUST
1320	C COMPUTE THE JACOBIAN MATRIX.
1321	C PREPDF COMPUTES THE INHOMOGENEITY MATRIX DF/DP.
1322	C SPRIME DEFINES THE SYSTEM OF SENSITIVITY EQUATIONS.
1323	C SOLSY MANAGES SOLUTION OF LINEAR SYSTEM IN CHORD ITERATION.
1324	C EWSET SETS THE ERROR WEIGHT VECTOR EWT BEFORE EACH STEP.
1325	C VNORM COMPUTES THE WEIGHTED R.M.S. NORM OF A VECTOR.
1326	C SVCOM AND RSCOM ARE USER-CALLABLE ROUTINES TO SAVE AND RESTORE,
1327	C RESPECTIVELY, THE CONTENTS OF THE INTERNAL COMMON BLOCKS.
1328	C DGEFA AND DGESL ARE ROUTINES FROM LINPACK FOR SOLVING FULL
1329	C SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS.
1330	C DGBFA AND DGBSL ARE ROUTINES FROM LINPACK FOR SOLVING BANDED
1331	C LINEAR SYSTEMS.
1332	C DAXPY, DSCAL, IDAMAX, AND DDOT ARE BASIC LINEAR ALGEBRA MODULES
1333	C (BLAS) USED BY THE ABOVE LINPACK ROUTINES.
1334	C D1MACH COMPUTES THE UNIT ROUNDOFF IN A MACHINE-INDEPENDENT MANNER.
1335	C XERR, XSETUN, AND XSETF HANDLE THE PRINTING OF ALL ERROR
1336	C MESSAGES AND WARNINGS.
1337	C NOTE.. VNORM, IDAMAX, DDOT, AND D1MACH ARE FUNCTION ROUTINES.
1338	C ALL THE OTHERS ARE SUBROUTINES.
1339	C
1340	C THE FORTRAN GENERIC INTRINSIC FUNCTIONS USED BY ODESSA ARE..
1341	C ABS, MAX, MIN, REAL, MOD, SIGN, SQRT, AND WRITE
1342	C
1343	C A BLOCK DATA SUBPROGRAM IS ALSO INCLUDED WITH THE PACKAGE,
1344	C FOR LOADING SOME OF THE VARIABLES IN INTERNAL COMMON.
1345	C
1346	C----------------------------------------------------------------------
1347	C PART V. GENERAL REMARKS
1348	C
1349	C THIS SECTION HIGHLIGHTS THE BASIC DIFFERENCES BETWEEN THE ORIGINAL
1350	C LSODE PACKAGE AND THE ODESSA MODIFICATION. THIS IS PROVIDED AS A
1351	C SERVICE TO EXPERIENCED LSODE USERS TO EXPEDITE FAMILIARIZATION WITH
1352	C ODESSA.
1353	C
1354	C (A). ORIGINAL SUBROUTINES AND FUNCTIONS.
1355	C
1356	C OF THE ORIGINAL 22 SUBROUTINES AND FUNCTIONS USED IN THE LSODE
1357	C PACKAGE, ALL ARE USED BY ODESSA, WITH THE FOLLOWING HAVING BEEN
1358	C MODIFIED..
1359	C
1360	C LSODE THE ORIGINAL DRIVER SUBROUTINE FOR THE LSODE PACKAGE IS
1361	C EXTENSIVELY MODIFIED AND RENAMED ODESSA, WHICH NOW
1362	C CONTAINS A CALL TO SPRIME TO ESTABLISH INITIAL CONDITIONS
1363	C FOR THE SENSITIVITY CALCULATIONS.
1364	C
1365	C STODE THE ONE STEP INTEGRATOR IS SLIGHTLY MODIFIED AND RETAINS
1366	C ITS ORIGINAL NAME. IT NOW CONTAINS THE CALL TO STESA,
1367	C AND ALSO CALLS SPRIME IF KFLAG .LE. -3.
1368	C
1369	C PREPJ ALSO NAMED PREPJ IN ODESSA IS SLIGHTLY MODIFIED TO ALLOW
1370	C FOR THE CALCULATION OF JACOBIAN WITH NO PREPROCESSING
1371	C (JOPT = 1).
1372	C
1373	C (B). NEW SUBROUTINES.
1374	C
1375	C IN ADDITION TO THE CHANGES NOTED ABOVE, THREE NEW SUBROUTINES
1376	C HAVE BEEN INTRODUCED (SEE STESA, SPRIME, AND PREPDF AS DESCRIBED
1377	C IN PART IV. ABOVE).
1378	C
1379	C (C). COMMON BLOCKS.
1380	C
1381	C /LS0001/ RETAINS THE SAME LENGTH AND IS RENAMED /ODE001/;
1382	C HOWEVER THE REAL ARRAY ROWNS(209) IS SHORTENED TO A
1383	C LENGTH OF (173) REAL WORDS, ALLOWING THE REMOVAL OF
1384	C TESCO(3,12) WHICH IS NOW PASSED FROM STODE TO STESA.
1385	C IN ADDITION, THE INTEGER ARRAY IOWNS(6) IS SHORTENED
1386	C TO A LENGTH OF (4) INTEGER WORDS, ALLOWING THE REMOVAL
1387	C OF IALTH AND LMAX WHICH ARE NOW PASSED FROM STODE TO
1388	C STESA.
1389	C
1390	C /ODE002/ ADDED COMMON BLOCK FOR VARIABLES IMPORTANT TO
1391	C SENSITIVITY ANALYSIS (SEE PART III. ABOVE). A BLOCK
1392	C DATA PROGRAM IS NOT REQUIRED FOR THIS COMMON BLOCK.
1393	C
1394	C SVCOM,RSCOM THESE TWO SUBROUTINES ARE MODIFIED TO HANDLE
1395	C COMMON BLOCK /ODE002/ AS WELL.
1396	C
1397	C (D). OPTIONAL INPUTS.
1398	C
1399	C THE FULL SET OF OPTIONAL INPUTS AVAILABLE IN LSODE IS ALSO
1400	C AVAILABLE IN ODESSA, WITH THE EXCEPTION THAT THE NUMBER OF ODE'S
1401	C IN THE MODEL (NEQ(1)), MAY NOT BE CHANGED DURING THE PROBLEM.
1402	C IN ODESSA, NYH NOW REFERS TO THE TOTAL NUMBER OF FIRST-ORDER
1403	C ODE'S (MODEL AND SENSITIVITY EQUATIONS) WHICH IS EQUAL TO
1404	C NEQ(1) IF ISOPT = 0, OR NEQ(1)*(NEQ(2)+1) IF ISOPT = 1.
1405	C NEQ(1), NEQ(2), AND NYH ARE NOT ALLOWED TO CHANGE DURING
1406	C THE COURSE OF AN INTEGRATION.
1407	C
1408	C (E). OPTIONAL OUTPUTS.
1409	C
1410	C THE FULL SET OF OPTIONAL OUTPUTS AVAILABLE IN LSODE IS ALSO
1411	C AVAILABLE IN ODESSA. IN ADDITION, IWORK(19) AND IWORK(20) ARE
1412	C LOADED WITH NDFE AND NSPE, RESPECTIVELY, UPON OUTPUT. THE TOTAL
1413	C NUMBER OF LU DECOMPOSITIONS OF THE PROCESSED JACOBIAN IS EQUAL
1414	C TO NJE - NSPE.
1415	C-----------------------------------------------------------------------
1416	SUBROUTINE KPP_ODESSA (F, DF, NEQ, Y, PAR, T, TOUT,
1417	1 ITOL, RTOL, ATOL,
1418	1 ITASK, ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF)
1419	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1420	LOGICAL IHIT
1421	EXTERNAL F, DF, JAC, PREPJ, SOLSY, PREPDF
1422	DIMENSION NEQ(), Y(), PAR(), RTOL(), ATOL(), IOPT(),
1423	1 RWORK(LRW), IWORK(LIW), MORD(2)
1424	C-----------------------------------------------------------------------
1425	C THIS IS THE SEPTEMBER 1, 1986 VERSION OF ODESSA..
1426	C AN ORDINARY DIFFERENTIAL EQUATION KppSolveR WITH EXPLICIT SIMULTANEOUS
1427	C SENSITIVITY ANALYSIS.
1428	C
1429	C THIS PACKAGE IS A MODIFICATION OF THE AUGUST 13, 1981 VERSION OF
1430	C LSODE.. LIVERMORE KppSolveR FOR ORDINARY DIFFERENTIAL EQUATIONS.
1431	C THIS VERSION IS IN DOUBLE PRECISION.
1432	C
1433	C ODESSA KppSolveS FOR THE FIRST-ORDER SENSITIVITY COEFFICIENTS..
1434	C DY(I)/DP, FOR A SINGLE PARAMETER, OR,
1435	C DY(I)/DP(J), FOR MULTIPLE PARAMETERS,
1436	C ASSOCIATED WITH A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS..
1437	C DY(T)/DT = F(Y,T;P).
1438	C-----------------------------------------------------------------------
1439	C REFERENCES...
1440	C
1441	C 1. JORGE R. LEIS AND MARK A. KRAMER, THE SIMULTANEOUS SOLUTION AND
1442	C EXPLICIT SENSITIVITY ANALYSIS OF SYSTEMS DESCRIBED BY ORDINARY
1443	C DIFFERENTIAL EQUATIONS, SUBMITTED TO ACM TRANS. MATH. SOFTWARE,
1444	C (1985).
1445	C
1446	C 2. JORGE R. LEIS AND MARK A. KRAMER, ODESSA - AN ORDINARY
1447	C DIFFERENTIAL EQUATION KppSolveR WITH EXPLICIT SIMULTANEOUS
1448	C SENSITIVITY ANALYSIS, SUBMITTED TO ACM TRANS. MATH. SOFTWARE.
1449	C (1985).
1450	C
1451	C 3. ALAN C. HINDMARSH, LSODE AND LSODI, TWO NEW INITIAL VALUE
1452	C ORDINARY DIFFERENTIAL EQUATION KppSolveRS, ACM-SIGNUM NEWSLETTER,
1453	C VOL. 15, NO. 4 (1980), PP. 10-11.
1454	C-----------------------------------------------------------------------
1455	C THE FOLLOWING INTERNAL COMMON BLOCKS CONTAIN
1456	C (A) VARIABLES WHICH ARE LOCAL TO ANY SUBROUTINE BUT WHOSE VALUES MUST
1457	C BE PRESERVED BETWEEN CALLS TO THE ROUTINE (OWN VARIABLES), AND
1458	C (B) VARIABLES WHICH ARE COMMUNICATED BETWEEN SUBROUTINES.
1459	C THE STRUCTURE OF THE BLOCKS ARE AS FOLLOWS.. ALL REAL VARIABLES ARE
1460	C LISTED FIRST, FOLLOWED BY ALL INTEGERS. WITHIN EACH TYPE, THE
1461	C VARIABLES ARE GROUPED WITH THOSE LOCAL TO SUBROUTINE ODESSA FIRST,
1462	C THEN THOSE LOCAL TO SUBROUTINE STODE, AND FINALLY THOSE USED
1463	C FOR COMMUNICATION. THE BLOCKS ARE DECLARED IN SUBROUTINES ODESSA
1464	C INTDY, STODE, STESA, PREPJ, PREPDF, AND SOLSY. GROUPS OF VARIABLES
1465	C ARE REPLACED BY DUMMY ARRAYS IN THE COMMON DECLARATIONS IN ROUTINES
1466	C WHERE THOSE VARIABLES ARE NOT USED.
1467	C-----------------------------------------------------------------------
1468	COMMON /ODE001/ TRET, ROWNS(173),
1469	1 TESCO(3,12), CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
1470	2 ILLIN, INIT, LYH, LEWT, LACOR, LSAVF, LWM, LIWM,
1471	3 MXSTEP, MXHNIL, NHNIL, NTREP, NSLAST, NYH, IOWNS(4),
1472	4 IALTH, LMAX, ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH,
1473	5 MITER, MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
1474	COMMON /ODE002/ DUPS, DSMS, DDNS,
1475	1 NPAR, LDFDP, LNRS,
1476	2 ISOPT, NSV, NDFE, NSPE, IDF, IERSP, JOPT, KFLAGS
1477	PARAMETER (ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0,FOUR=4.0D0)
1478	DATA MORD(1),MORD(2)/12,5/, MXSTP0/500/, MXHNL0/10/
1479	C-----------------------------------------------------------------------
1480	C BLOCK A.
1481	C THIS CODE BLOCK IS EXECUTED ON EVERY CALL.
1482	C IT TESTS ISTATE AND ITASK FOR LEGALITY AND BRANCHES APPROPIATELY.
1483	C IF ISTATE .GT. 1 BUT THE FLAG INIT SHOWS THAT INITIALIZATION HAS
1484	C NOT YET BEEN DONE, AN ERROR RETURN OCCURS.
1485	C IF ISTATE = 1 AND TOUT = T, JUMP TO BLOCK G AND RETURN IMMEDIATELY.
1486	C-----------------------------------------------------------------------
1487	IF (ISTATE .LT. 1 .OR. ISTATE .GT. 3) GO TO 601
1488	IF (ITASK .LT. 1 .OR. ITASK .GT. 5) GO TO 602
1489	IF (ISTATE .EQ. 1) GO TO 10
1490	IF (INIT .EQ. 0) GO TO 603
1491	IF (ISTATE .EQ. 2) GO TO 200
1492	GO TO 20
1493	10 INIT = 0
1494	IF (TOUT .EQ. T) GO TO 430
1495	20 NTREP = 0
1496	C-----------------------------------------------------------------------
1497	C BLOCK B.
1498	C THE NEXT CODE BLOCK IS EXECUTED FOR THE INITIAL CALL (ISTATE = 1),
1499	C OR FOR A CONTINUATION CALL WITH PARAMETER CHANGES (ISTATE = 3).
1500	C IT CONTAINS CHECKING OF ALL INPUTS AND VARIOUS INITIALIZATIONS.
1501	C
1502	C FIRST CHECK LEGALITY OF THE NON-OPTIONAL INPUTS NEQ, ITOL, IOPT,
1503	C MF, ML, AND MU.
1504	C-----------------------------------------------------------------------
1505	IF (NEQ(1) .LE. 0) GO TO 604
1506	IF (ISTATE .EQ. 1) GO TO 25
1507	IF (NEQ(1) .NE. N) GO TO 605
1508	25 N = NEQ(1)
1509	IF (ITOL .LT. 1 .OR. ITOL .GT. 4) GO TO 606
1510	DO 26 I = 1,3
1511	26 IF (IOPT(I) .LT. 0 .OR. IOPT(I) .GT. 1) GO TO 607
1512	ISOPT = IOPT(2)
1513	IDF = IOPT(3)
1514	NYH = N
1515	NSV = 1
1516	METH = MF/10
1517	MITER = MF - 10*METH
1518	IF (METH .LT. 1 .OR. METH .GT. 2) GO TO 608
1519	IF (MITER .LT. 0 .OR. MITER .GT. 5) GO TO 608
1520	IF (MITER .LE. 3) GO TO 30
1521	ML = IWORK(1)
1522	MU = IWORK(2)
1523	IF (ML .LT. 0 .OR. ML .GE. N) GO TO 609
1524	IF (MU .LT. 0 .OR. MU .GE. N) GO TO 610
1525	30 IF (ISOPT .EQ. 0) GO TO 32
1526	C CHECK LEGALITY OF THE NON-OPTIONAL INPUTS ISOPT, NPAR.
1527	C COMPUTE NUMBER OF SOLUTION VECTORS AND TOTAL NUMBER OF EQUATIONS.
1528	IF (NEQ(2) .LE. 0) GO TO 628
1529	IF (ISTATE .EQ. 1) GO TO 31
1530	IF (NEQ(2) .NE. NPAR) GO TO 629
1531	31 NPAR = NEQ(2)
1532	NSV = NPAR + 1
1533	NYH = NSV * N
1534	IF (MITER .EQ. 0 .OR. MITER .EQ. 3) GO TO 630
1535	C NEXT PROCESS AND CHECK THE OPTIONAL INPUTS. --------------------------
1536	32 IF (IOPT(1) .EQ. 1) GO TO 40
1537	MAXORD = MORD(METH)
1538	MXSTEP = MXSTP0
1539	MXHNIL = MXHNL0
1540	IF (ISTATE .EQ. 1) H0 = ZERO
1541	HMXI = ZERO
1542	HMIN = ZERO
1543	GO TO 60
1544	40 MAXORD = IWORK(5)
1545	IF (MAXORD .LT. 0) GO TO 611
1546	IF (MAXORD .EQ. 0) MAXORD = 100
1547	MAXORD = MIN(MAXORD,MORD(METH))
1548	MXSTEP = IWORK(6)
1549	IF (MXSTEP .LT. 0) GO TO 612
1550	IF (MXSTEP .EQ. 0) MXSTEP = MXSTP0
1551	MXHNIL = IWORK(7)
1552	IF (MXHNIL .LT. 0) GO TO 613
1553	IF (MXHNIL .EQ. 0) MXHNIL = MXHNL0
1554	IF (ISTATE .NE. 1) GO TO 50
1555	H0 = RWORK(5)
1556	IF ((TOUT - T)*H0 .LT. ZERO) GO TO 614
1557	50 HMAX = RWORK(6)
1558	IF (HMAX .LT. ZERO) GO TO 615
1559	HMXI = ZERO
1560	IF (HMAX .GT. ZERO) HMXI = ONE/HMAX
1561	HMIN = RWORK(7)
1562	IF (HMIN .LT. ZERO) GO TO 616
1563	C-----------------------------------------------------------------------
1564	C SET WORK ARRAY POINTERS AND CHECK LENGTHS LRW AND LIW.
1565	C POINTERS TO SEGMENTS OF RWORK AND IWORK ARE NAMED BY PREFIXING L TO
1566	C THE NAME OF THE SEGMENT. E.G., THE SEGMENT YH STARTS AT RWORK(LYH).
1567	C SEGMENTS OF RWORK (IN ORDER) ARE DENOTED YH, WM, EWT, SAVF, ACOR.
1568	C WORK SPACE FOR DFDP IS CONTAINED IN ACOR.
1569	C-----------------------------------------------------------------------
1570	60 LYH = 21
1571	LWM = LYH + (MAXORD + 1)*NYH
1572	IF (MITER .EQ. 0) LENWM = 0
1573	IF (MITER .EQ. 1 .OR. MITER .EQ. 2) LENWM = N*N + 2
1574	IF (MITER .EQ. 3) LENWM = N + 2
1575	IF (MITER .GE. 4) LENWM = (2ML + MU + 1)N + 2
1576	LEWT = LWM + LENWM
1577	LSAVF = LEWT + NYH
1578	LACOR = LSAVF + N
1579	LDFDP = LACOR + N
1580	LENRW = LACOR + NYH - 1
1581	IWORK(17) = LENRW
1582	LIWM = 1
1583	LENIW = 20 + N
1584	IF (MITER .EQ. 0 .OR. MITER .EQ. 3) LENIW = 20
1585	LNRS = LENIW + LIWM
1586	IF (ISOPT .EQ. 1) LENIW = LNRS + NPAR
1587	IWORK(18) = LENIW
1588	IF (LENRW .GT. LRW) GO TO 617
1589	IF (LENIW .GT. LIW) GO TO 618
1590	C CHECK RTOL AND ATOL FOR LEGALITY. ------------------------------------
1591	RTOLI = RTOL(1)
1592	ATOLI = ATOL(1)
1593	DO 70 I = 1,NYH
1594	IF (ITOL .GE. 3) RTOLI = RTOL(I)
1595	IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) ATOLI = ATOL(I)
1596	IF (RTOLI .LT. ZERO) GO TO 619
1597	IF (ATOLI .LT. ZERO) GO TO 620
1598	70 CONTINUE
1599	IF (ISTATE .EQ. 1) GO TO 100
1600	C IF ISTATE = 3, SET FLAG TO SIGNAL PARAMETER CHANGES TO STODE. --------
1601	JSTART = -1
1602	IF (NQ .LE. MAXORD) GO TO 90
1603	C MAXORD WAS REDUCED BELOW NQ. COPY YH(*,MAXORD+2) INTO SAVF. ---------
1604	DO 80 I = 1,N
1605	80 RWORK(I+LSAVF-1) = RWORK(I+LWM-1)
1606	C RELOAD WM(1) = RWORK(LWM), SINCE LWM MAY HAVE CHANGED. ---------------
1607	90 IF (MITER .GT. 0) RWORK(LWM) = SQRT(UROUND)
1608	GO TO 200
1609	C-----------------------------------------------------------------------
1610	C BLOCK C.
1611	C THE NEXT BLOCK IS FOR THE INITIAL CALL ONLY (ISTATE = 1).
1612	C IT CONTAINS ALL REMAINING INITIALIZATIONS, THE INITIAL CALL TO F,
1613	C THE INITIAL CALL TO SPRIME IF ISOPT = 1,
1614	C AND THE CALCULATION OF THE INITIAL STEP SIZE.
1615	C THE ERROR WEIGHTS IN EWT ARE INVERTED AFTER BEING LOADED.
1616	C-----------------------------------------------------------------------
1617	100 UROUND = D1MACH(4)
1618	TN = T
1619	IF (ITASK .NE. 4 .AND. ITASK .NE. 5) GO TO 105
1620	TCRIT = RWORK(1)
1621	IF ((TCRIT - TOUT)*(TOUT - T) .LT. ZERO) GO TO 625
1622	IF (H0 .NE. ZERO .AND. (T + H0 - TCRIT)*H0 .GT. ZERO)
1623	1 H0 = TCRIT - T
1624	105 JSTART = 0
1625	IF (MITER .GT. 0) RWORK(LWM) = SQRT(UROUND)
1626	NHNIL = 0
1627	NST = 0
1628	NJE = 0
1629	NSLAST = 0
1630	HU = ZERO
1631	NQU = 0
1632	CCMAX = 0.3D0
1633	MAXCOR = 3
1634	IF (ISOPT .EQ. 1) MAXCOR = 4
1635	MSBP = 20
1636	MXNCF = 10
1637	C INITIAL CALL TO F. (LF0 POINTS TO YH(1,2) AND LOADS IN VALUES).
1638	LF0 = LYH + NYH
1639	CALL F (NEQ, T, Y, PAR, RWORK(LF0))
1640	NFE = 1
1641	DUPS = ZERO
1642	DSMS = ZERO
1643	DDNS = ZERO
1644	NDFE = 0
1645	NSPE = 0
1646	IF (ISOPT .EQ. 0) GO TO 114
1647	C INITIALIZE COUNTS FOR REPEATED STEPS DUE TO SENSITIVITY ANALYSIS.
1648	DO 110 J = 1,NSV
1649	110 IWORK(J + LNRS - 1) = 0
1650	C LOAD THE INITIAL VALUE VECTOR IN YH. ---------------------------------
1651	114 DO 115 I = 1,NYH
1652	115 RWORK(I+LYH-1) = Y(I)
1653	C LOAD AND INVERT THE EWT ARRAY. (H IS TEMPORARILY SET TO ONE.) -------
1654	NQ = 1
1655	H = ONE
1656	CALL EWSET (NYH, ITOL, RTOL, ATOL, RWORK(LYH), RWORK(LEWT))
1657	DO 120 I = 1,NYH
1658	IF (RWORK(I+LEWT-1) .LE. ZERO) GO TO 621
1659	120 RWORK(I+LEWT-1) = ONE/RWORK(I+LEWT-1)
1660	IF (ISOPT .EQ. 0) GO TO 125
1661	C CALL SPRIME TO LOAD FIRST-ORDER SENSITIVITY DERIVATIVES INTO
1662	C REMAINING YH(*,2) POSITIONS.
1663	CALL SPRIME (NEQ, Y, RWORK(LYH), NYH, N, NSV, RWORK(LWM),
1664	1 IWORK(LIWM), RWORK(LEWT), RWORK(LF0), RWORK(LACOR),
1665	2 RWORK(LDFDP), PAR, F, JAC, DF, PREPJ, PREPDF)
1666	IF (IERSP .EQ. -1) GO TO 631
1667	IF (IERSP .EQ. -2) GO TO 632
1668	C-----------------------------------------------------------------------
1669	C THE CODING BELOW COMPUTES THE STEP SIZE, H0, TO BE ATTEMPTED ON THE
1670	C FIRST STEP, UNLESS THE USER HAS SUPPLIED A VALUE FOR THIS.
1671	C FIRST CHECK THAT TOUT - T DIFFERS SIGNIFICANTLY FROM ZERO.
1672	C A SCALAR TOLERANCE QUANTITY TOL IS COMPUTED, AS MAX(RTOL(I))
1673	C IF THIS IS POSITIVE, OR MAX(ATOL(I)/ABS(Y(I))) OTHERWISE, ADJUSTED
1674	C SO AS TO BE BETWEEN 100*UROUND AND 1.0E-3. ONLY THE ORIGINAL
1675	C SOLUTION VECTOR IS CONSIDERED IN THIS CALCULATION (ISOPT = 0 OR 1).
1676	C THEN THE COMPUTED VALUE H0 IS GIVEN BY..
1677	C NEQ
1678	C H02 = TOL / ( W0-2 + (1/NEQ) * SUM ( F(I)/YWT(I) )**2 )
1679	C 1
1680	C WHERE W0 = MAX ( ABS(T), ABS(TOUT) ),
1681	C F(I) = I-TH COMPONENT OF INITIAL VALUE OF F,
1682	C YWT(I) = EWT(I)/TOL (A WEIGHT FOR Y(I)).
1683	C THE SIGN OF H0 IS INFERRED FROM THE INITIAL VALUES OF TOUT AND T.
1684	C-----------------------------------------------------------------------
1685	125 IF (H0 .NE. ZERO) GO TO 180
1686	TDIST = ABS(TOUT - T)
1687	W0 = MAX(ABS(T),ABS(TOUT))
1688	IF (TDIST .LT. TWOUROUNDW0) GO TO 622
1689	TOL = RTOL(1)
1690	IF (ITOL .LE. 2) GO TO 140
1691	DO 130 I = 1,N
1692	130 TOL = MAX(TOL,RTOL(I))
1693	140 IF (TOL .GT. ZERO) GO TO 160
1694	ATOLI = ATOL(1)
1695	DO 150 I = 1,N
1696	IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) ATOLI = ATOL(I)
1697	AYI = ABS(Y(I))
1698	IF (AYI .NE. ZERO) TOL = MAX(TOL,ATOLI/AYI)
1699	150 CONTINUE
1700	160 TOL = MAX(TOL,100.0D0*UROUND)
1701	TOL = MIN(TOL,0.001D0)
1702	SUM = VNORM (N, RWORK(LF0), RWORK(LEWT))
1703	SUM = ONE/(TOLW0W0) + TOLSUM*2
1704	H0 = ONE/SQRT(SUM)
1705	H0 = MIN(H0,TDIST)
1706	H0 = SIGN(H0,TOUT-T)
1707	C ADJUST H0 IF NECESSARY TO MEET HMAX BOUND. ---------------------------
1708	180 RH = ABS(H0)*HMXI
1709	IF (RH .GT. ONE) H0 = H0/RH
1710	C LOAD H WITH H0 AND SCALE YH(*,2) BY H0. ------------------------------
1711	H = H0
1712	DO 190 I = 1,NYH
1713	190 RWORK(I+LF0-1) = H0*RWORK(I+LF0-1)
1714	GO TO 270
1715	C-----------------------------------------------------------------------
1716	C BLOCK D.
1717	C THE NEXT CODE BLOCK IS FOR CONTINUATION CALLS ONLY (ISTATE = 2 OR 3)
1718	C AND IS TO CHECK STOP CONDITIONS BEFORE TAKING A STEP.
1719	C-----------------------------------------------------------------------
1720	200 NSLAST = NST
1721	GO TO (210, 250, 220, 230, 240), ITASK
1722	210 IF ((TN - TOUT)*H .LT. ZERO) GO TO 250
1723	CALL INTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG)
1724	IF (IFLAG .NE. 0) GO TO 627
1725	T = TOUT
1726	GO TO 420
1727	220 TP = TN - HU(ONE + 100.0D0UROUND)
1728	IF ((TP - TOUT)*H .GT. ZERO) GO TO 623
1729	IF ((TN - TOUT)*H .LT. ZERO) GO TO 250
1730	GO TO 400
1731	230 TCRIT = RWORK(1)
1732	IF ((TN - TCRIT)*H .GT. ZERO) GO TO 624
1733	IF ((TCRIT - TOUT)*H .LT. ZERO) GO TO 625
1734	IF ((TN - TOUT)*H .LT. ZERO) GO TO 245
1735	CALL INTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG)
1736	IF (IFLAG .NE. 0) GO TO 627
1737	T = TOUT
1738	GO TO 420
1739	240 TCRIT = RWORK(1)
1740	IF ((TN - TCRIT)*H .GT. ZERO) GO TO 624
1741	245 HMX = ABS(TN) + ABS(H)
1742	IHIT = ABS(TN - TCRIT) .LE. 100.0D0UROUNDHMX
1743	IF (IHIT) GO TO 400
1744	TNEXT = TN + H(ONE + FOURUROUND)
1745	IF ((TNEXT - TCRIT)*H .LE. ZERO) GO TO 250
1746	H = (TCRIT - TN)(ONE - FOURUROUND)
1747	IF (ISTATE .EQ. 2) JSTART = -2
1748	C-----------------------------------------------------------------------
1749	C BLOCK E.
1750	C THE NEXT BLOCK IS NORMALLY EXECUTED FOR ALL CALLS AND CONTAINS
1751	C THE CALL TO THE ONE-STEP CORE INTEGRATOR STODE.
1752	C
1753	C THIS IS A LOOPING POINT FOR THE INTEGRATION STEPS.
1754	C
1755	C FIRST CHECK FOR TOO MANY STEPS BEING TAKEN, UPDATE EWT (IF NOT AT
1756	C START OF PROBLEM), CHECK FOR TOO MUCH ACCURACY BEING REQUESTED, AND
1757	C CHECK FOR H BELOW THE ROUNDOFF LEVEL IN T.
1758	C TOLSF IS CALCULATED CONSIDERING ALL SOLUTION VECTORS.
1759	C-----------------------------------------------------------------------
1760	250 CONTINUE
1761	IF ((NST-NSLAST) .GE. MXSTEP) GO TO 500
1762	CALL EWSET (NYH, ITOL, RTOL, ATOL, RWORK(LYH), RWORK(LEWT))
1763	DO 260 I = 1,NYH
1764	IF (RWORK(I+LEWT-1) .LE. ZERO) GO TO 510
1765	260 RWORK(I+LEWT-1) = ONE/RWORK(I+LEWT-1)
1766	270 TOLSF = UROUND*VNORM (NYH, RWORK(LYH), RWORK(LEWT))
1767	IF (TOLSF .LE. ONE) GO TO 280
1768	TOLSF = TOLSF*2.0D0
1769	IF (NST .EQ. 0) GO TO 626
1770	GO TO 520
1771	280 IF (ADDX(TN,H) .NE. TN) GO TO 290
1772	NHNIL = NHNIL + 1
1773	IF (NHNIL .GT. MXHNIL) GO TO 290
1774	CALL XERR ('ODESSA - WARNING..INTERNAL T (=R1) AND H (=R2) ARE',
1775	1 101, 1, 0, 0, 0, 0, ZERO, ZERO)
1776	CALL XERR ('SUCH THAT IN THE MACHINE, T + H = T ON THE NEXT STEP',
1777	1 101, 1, 0, 0, 0, 0, ZERO, ZERO)
1778	CALL XERR ('(H = STEP SIZE). KppSolveR WILL CONTINUE ANYWAY',
1779	1 101, 1, 0, 0, 0, 2, TN, H)
1780	IF (NHNIL .LT. MXHNIL) GO TO 290
1781	CALL XERR ('ODESSA - ABOVE WARNING HAS BEEN ISSUED I1 TIMES.',
1782	1 102, 1, 0, 0, 0, 0, ZERO, ZERO)
1783	CALL XERR ('IT WILL NOT BE ISSUED AGAIN FOR THIS PROBLEM',
1784	1 102, 1, 1, MXHNIL, 0, 0, ZERO,ZERO)
1785	290 CONTINUE
1786	C-----------------------------------------------------------------------
1787	C CALL STODE(NEQ,Y,YH,NYH,YH,WM,IWM,EWT,SAVF,ACOR,PAR,NRS,
1788	C 1 F,JAC,DF,PREPJ,PREPDF,SOLSY)
1789	C-----------------------------------------------------------------------
1790	CALL STODE (NEQ, Y, RWORK(LYH), NYH, RWORK(LYH), RWORK(LWM),
1791	1 IWORK(LIWM), RWORK(LEWT), RWORK(LSAVF), RWORK(LACOR),
1792	2 PAR, IWORK(LNRS), F, JAC, DF, PREPJ, PREPDF, SOLSY)
1793	KGO = 1 - KFLAG
1794	GO TO (300, 530, 540, 633), KGO
1795	C-----------------------------------------------------------------------
1796	C BLOCK F.
1797	C THE FOLLOWING BLOCK HANDLES THE CASE OF A SUCCESSFUL RETURN FROM THE
1798	C CORE INTEGRATOR (KFLAG = 0). TEST FOR STOP CONDITIONS.
1799	C-----------------------------------------------------------------------
1800	300 INIT = 1
1801	GO TO (310, 400, 330, 340, 350), ITASK
1802	C ITASK = 1. IF TOUT HAS BEEN REACHED, INTERPOLATE. -------------------
1803	310 IF ((TN - TOUT)*H .LT. ZERO) GO TO 250
1804	CALL INTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG)
1805	T = TOUT
1806	GO TO 420
1807	C ITASK = 3. JUMP TO EXIT IF TOUT WAS REACHED. ------------------------
1808	330 IF ((TN - TOUT)*H .GE. ZERO) GO TO 400
1809	GO TO 250
1810	C ITASK = 4. SEE IF TOUT OR TCRIT WAS REACHED. ADJUST H IF NECESSARY.
1811	340 IF ((TN - TOUT)*H .LT. ZERO) GO TO 345
1812	CALL INTDY (TOUT, 0, RWORK(LYH), NYH, Y, IFLAG)
1813	T = TOUT
1814	GO TO 420
1815	345 HMX = ABS(TN) + ABS(H)
1816	IHIT = ABS(TN - TCRIT) .LE. 100.0D0UROUNDHMX
1817	IF (IHIT) GO TO 400
1818	TNEXT = TN + H(ONE + FOURUROUND)
1819	IF ((TNEXT - TCRIT)*H .LE. ZERO) GO TO 250
1820	H = (TCRIT - TN)(ONE - FOURUROUND)
1821	JSTART = -2
1822	GO TO 250
1823	C ITASK = 5. SEE IF TCRIT WAS REACHED AND JUMP TO EXIT. ---------------
1824	350 HMX = ABS(TN) + ABS(H)
1825	IHIT = ABS(TN - TCRIT) .LE. 100.0D0UROUNDHMX
1826	C-----------------------------------------------------------------------
1827	C BLOCK G.
1828	C THE FOLLOWING BLOCK HANDLES ALL SUCCESSFUL RETURNS FROM ODESSA.
1829	C IF ITASK .NE. 1, Y IS LOADED FROM YH AND T IS SET ACCORDINGLY.
1830	C ISTATE IS SET TO 2, THE ILLEGAL INPUT COUNTER IS ZEROED, AND THE
1831	C OPTIONAL OUTPUTS ARE LOADED INTO THE WORK ARRAYS BEFORE RETURNING.
1832	C IF ISTATE = 1 AND TOUT = T, THERE IS A RETURN WITH NO ACTION TAKEN,
1833	C EXCEPT THAT IF THIS HAS HAPPENED REPEATEDLY, THE RUN IS TERMINATED.
1834	C-----------------------------------------------------------------------
1835	400 DO 410 I = 1,NYH
1836	410 Y(I) = RWORK(I+LYH-1)
1837	T = TN
1838	IF (ITASK .NE. 4 .AND. ITASK .NE. 5) GO TO 420
1839	IF (IHIT) T = TCRIT
1840	420 ISTATE = 2
1841	ILLIN = 0
1842	RWORK(11) = HU
1843	RWORK(12) = H
1844	RWORK(13) = TN
1845	IWORK(11) = NST
1846	IWORK(12) = NFE
1847	IWORK(13) = NJE
1848	IWORK(14) = NQU
1849	IWORK(15) = NQ
1850	IF (ISOPT .EQ. 0) RETURN
1851	IWORK(19) = NDFE
1852	IWORK(20) = NSPE
1853	RETURN
1854	430 NTREP = NTREP + 1
1855	IF (NTREP .LT. 5) RETURN
1856	CALL XERR ('ODESSA -- REPEATED CALLS WITH ISTATE = 1 AND
1857	1TOUT = T (=R1)', 301, 1, 0, 0, 0, 1, T, ZERO)
1858	GO TO 800
1859	C-----------------------------------------------------------------------
1860	C BLOCK H.
1861	C THE FOLLOWING BLOCK HANDLES ALL UNSUCCESSFUL RETURNS OTHER THAN
1862	C THOSE FOR ILLEGAL INPUT. FIRST THE ERROR MESSAGE ROUTINE IS CALLED.
1863	C IF THERE WAS AN ERROR TEST OR CONVERGENCE TEST FAILURE, IMXER IS SET.
1864	C THEN Y IS LOADED FROM YH, T IS SET TO TN, AND THE ILLEGAL INPUT
1865	C COUNTER ILLIN IS SET TO 0. THE OPTIONAL OUTPUTS ARE LOADED INTO
1866	C THE WORK ARRAYS BEFORE RETURNING.
1867	C-----------------------------------------------------------------------
1868	C THE MAXIMUM NUMBER OF STEPS WAS TAKEN BEFORE REACHING TOUT. ----------
1869	500 CALL XERR ('ODESSA - AT CURRENT T (=R1), MXSTEP (=I1) STEPS',
1870	1 201, 1, 0, 0, 0, 0, ZERO,ZERO)
1871	CALL XERR ('TAKEN ON THIS CALL BEFORE REACHING TOUT',
1872	1 201, 1, 1, MXSTEP, 0, 1, TN, ZERO)
1873	ISTATE = -1
1874	GO TO 580
1875	C EWT(I) .LE. 0.0 FOR SOME I (NOT AT START OF PROBLEM). ----------------
1876	510 EWTI = RWORK(LEWT+I-1)
1877	CALL XERR ('ODESSA - AT T (=R1), EWT(I1) HAS BECOME R2 .LE. 0.',
1878	1 202, 1, 1, I, 0, 2, TN, EWTI)
1879	ISTATE = -6
1880	GO TO 580
1881	C TOO MUCH ACCURACY REQUESTED FOR MACHINE PRECISION. -------------------
1882	520 CALL XERR ('ODESSA - AT T (=R1), TOO MUCH ACCURACY REQUESTED',
1883	1 203, 1, 0, 0, 0, 0, ZERO,ZERO)
1884	CALL XERR ('FOR PRECISION OF MACHINE.. SEE TOLSF (=R2)',
1885	1 203, 1, 0, 0, 0, 2, TN, TOLSF)
1886	RWORK(14) = TOLSF
1887	ISTATE = -2
1888	GO TO 580
1889	C KFLAG = -1. ERROR TEST FAILED REPEATEDLY OR WITH ABS(H) = HMIN. -----
1890	530 CALL XERR ('ODESSA - AT T(=R1) AND STEP SIZE H(=R2), THE ERROR',
1891	1 204, 1, 0, 0, 0, 0, ZERO, ZERO)
1892	CALL XERR ('TEST FAILED REPEATEDLY OR WITH ABS(H) = HMIN',
1893	1 204, 1, 0, 0, 0, 2, TN, H)
1894	ISTATE = -4
1895	GO TO 560
1896	C KFLAG = -2. CONVERGENCE FAILED REPEATEDLY OR WITH ABS(H) = HMIN. ----
1897	540 CALL XERR ('ODESSA - AT T (=R1) AND STEP SIZE H (=R2), THE',
1898	1 205, 1, 0, 0, 0, 0, ZERO,ZERO)
1899	CALL XERR ('CORRECTOR CONVERGENCE FAILED REPEATEDLY',
1900	1 205, 1, 0, 0, 0, 0, ZERO,ZERO)
1901	CALL XERR ('OR WITH ABS(H) = HMIN',
1902	1 205, 1, 0, 0, 0, 2, TN, H)
1903	ISTATE = -5
1904	C COMPUTE IMXER IF RELEVANT. -------------------------------------------
1905	560 BIG = ZERO
1906	IMXER = 1
1907	DO 570 I = 1,NYH
1908	SIZE = ABS(RWORK(I+LACOR-1)*RWORK(I+LEWT-1))
1909	IF (BIG .GE. SIZE) GO TO 570
1910	BIG = SIZE
1911	IMXER = I
1912	570 CONTINUE
1913	IWORK(16) = IMXER
1914	C SET Y VECTOR, T, ILLIN, AND OPTIONAL OUTPUTS. ------------------------
1915	580 DO 590 I = 1,NYH
1916	590 Y(I) = RWORK(I+LYH-1)
1917	T = TN
1918	ILLIN = 0
1919	RWORK(11) = HU
1920	RWORK(12) = H
1921	RWORK(13) = TN
1922	IWORK(11) = NST
1923	IWORK(12) = NFE
1924	IWORK(13) = NJE
1925	IWORK(14) = NQU
1926	IWORK(15) = NQ
1927	IF (ISOPT .EQ. 0) RETURN
1928	IWORK(19) = NDFE
1929	IWORK(20) = NSPE
1930	RETURN
1931	C-----------------------------------------------------------------------
1932	C BLOCK I.
1933	C THE FOLLOWING BLOCK HANDLES ALL ERROR RETURNS DUE TO ILLEGAL INPUT
1934	C (ISTATE = -3), AS DETECTED BEFORE CALLING THE CORE INTEGRATOR.
1935	C FIRST THE ERROR MESSAGE ROUTINE IS CALLED. THEN IF THERE HAVE BEEN
1936	C 5 CONSECUTIVE SUCH RETURNS JUST BEFORE THIS CALL TO THE KppSolveR,
1937	C THE RUN IS HALTED.
1938	C-----------------------------------------------------------------------
1939	601 CALL XERR ('ODESSA - ISTATE (=I1) ILLEGAL',
1940	1 1, 1, 1, ISTATE, 0, 0, ZERO,ZERO)
1941	GO TO 700
1942	602 CALL XERR ('ODESSA - ITASK (=I1) ILLEGAL',
1943	1 2, 1, 1, ITASK, 0, 0, ZERO,ZERO)
1944	GO TO 700
1945	603 CALL XERR ('ODESSA - ISTATE .GT. 1 BUT ODESSA NOT INITIALIZED',
1946	1 3, 1, 0, 0, 0, 0, ZERO,ZERO)
1947	GO TO 700
1948	604 CALL XERR ('ODESSA - NEQ (=I1) .LT. 1',
1949	1 4, 1, 1, NEQ(1), 0, 0, ZERO,ZERO)
1950	GO TO 700
1951	605 CALL XERR ('ODESSA - ISTATE = 3 AND NEQ CHANGED. (I1 TO I2)',
1952	1 5, 1, 2, N, NEQ(1), 0, ZERO,ZERO)
1953	GO TO 700
1954	606 CALL XERR ('ODESSA - ITOL (=I1) ILLEGAL',
1955	1 6, 1, 1, ITOL, 0, 0, ZERO,ZERO)
1956	GO TO 700
1957	607 CALL XERR ('ODESSA - IOPT (=I1) ILLEGAL',
1958	1 7, 1, 1, IOPT, 0, 0, ZERO,ZERO)
1959	GO TO 700
1960	608 CALL XERR('ODESSA - MF (=I1) ILLEGAL',
1961	1 8, 1, 1, MF, 0, 0, ZERO,ZERO)
1962	GO TO 700
1963	609 CALL XERR('ODESSA - ML (=I1) ILLEGAL.. .LT.0 OR .GE.NEQ (=I2)',
1964	1 9, 1, 2, ML, NEQ(1), 0, ZERO,ZERO)
1965	GO TO 700
1966	610 CALL XERR('ODESSA - MU (=I1) ILLEGAL.. .LT.0 OR .GE.NEQ (=I2)',
1967	1 10, 1, 2, MU, NEQ(1), 0, ZERO,ZERO)
1968	GO TO 700
1969	611 CALL XERR('ODESSA - MAXORD (=I1) .LT. 0',
1970	1 11, 1, 1, MAXORD, 0, 0, ZERO,ZERO)
1971	GO TO 700
1972	612 CALL XERR('ODESSA - MXSTEP (=I1) .LT. 0',
1973	1 12, 1, 1, MXSTEP, 0, 0, ZERO,ZERO)
1974	GO TO 700
1975	613 CALL XERR('ODESSA - MXHNIL (=I1) .LT. 0',
1976	1 13, 1, 1, MXHNIL, 0, 0, ZERO,ZERO)
1977	GO TO 700
1978	614 CALL XERR('ODESSA - TOUT (=R1) BEHIND T (=R2)',
1979	1 14, 1, 0, 0, 0, 2, TOUT, T)
1980	CALL XERR('INTEGRATION DIRECTION IS GIVEN BY H0 (=R1)',
1981	1 14, 1, 0, 0, 0, 1, H0, ZERO)
1982	GO TO 700
1983	615 CALL XERR('ODESSA - HMAX (=R1) .LT. 0.0',
1984	1 15, 1, 0, 0, 0, 1, HMAX, ZERO)
1985	GO TO 700
1986	616 CALL XERR('ODESSA - HMIN (=R1) .LT. 0.0',
1987	1 16, 1, 0, 0, 0, 1, HMIN, ZERO)
1988	GO TO 700
1989	617 CALL XERR('ODESSA - RWORK LENGTH NEEDED, LENRW (=I1), EXCEEDS
1990	1 LRW (=I2)', 17, 1, 2, LENRW, LRW, 0, ZERO,ZERO)
1991	GO TO 700
1992	618 CALL XERR('ODESSA - IWORK LENGTH NEEDED, LENIW (=I1), EXCEEDS
1993	1 LIW (=I2)', 18, 1, 2, LENIW, LIW, 0, ZERO,ZERO)
1994	GO TO 700
1995	619 CALL XERR('ODESSA - RTOL(I1) IS R1 .LT. 0.0',
1996	1 19, 1, 1, I, 0, 1, RTOLI, ZREO)
1997	GO TO 700
1998	620 CALL XERR('ODESSA - ATOL(I1) IS R1 .LT. 0.0',
1999	1 20, 1, 1, I, 0, 1, ATOLI, ZERO)
2000	GO TO 700
2001	*
2002	621 EWTI = RWORK(LEWT+I-1)
2003	CALL XERR('ODESSA - EWT(I1) IS R1 .LE. 0.0',
2004	1 21, 1, 1, I, 0, 1, EWTI, ZERO)
2005	GO TO 700
2006	622 CALL XERR('ODESSA - TOUT (=R1) TOO CLOSE TO T(=R2) TO START
2007	1 INTEGRATION', 22, 1, 0, 0, 0, 2, TOUT, T)
2008	GO TO 700
2009	623 CALL XERR('ODESSA - ITASK = I1 AND TOUT (=R1) BEHIND TCUR - HU
2010	1 (= R2)', 23, 1, 1, ITASK, 0, 2, TOUT, TP)
2011	GO TO 700
2012	624 CALL XERR('ODESSA - ITASK = 4 OR 5 AND TCRIT (=R1) BEHIND TCUR
2013	1 (=R2)', 24, 1, 0, 0, 0, 2, TCRIT, TN)
2014	GO TO 700
2015	625 CALL XERR('ODESSA - ITASK = 4 OR 5 AND TCRIT (=R1) BEHIND TOUT
2016	1 (=R2)', 25, 1, 0, 0, 0, 2, TCRIT, TOUT)
2017	GO TO 700
2018	626 CALL XERR('ODESSA - AT START OF PROBLEM, TOO MUCH ACCURACY',
2019	1 26, 1, 0, 0, 0, 0, ZERO,ZERO)
2020	CALL XERR('REQUESTED FOR PRECISION OF MACHINE. SEE TOLSF (=R1)',
2021	1 26, 1, 0, 0, 0, 1, TOLSF, ZERO)
2022	RWORK(14) = TOLSF
2023	GO TO 700
2024	627 CALL XERR('ODESSA - TROUBLE FROM INTDY. ITASK = I1, TOUT = R1',
2025	1 27, 1, 1, ITASK, 0, 1, TOUT, ZERO)
2026	GO TO 700
2027	C ERROR STATEMENTS ASSOCIATED WITH SENSITIVITY ANALYSIS.
2028	628 CALL XERR('ODESSA - NPAR (=I1) .LT. 1',
2029	1 28, 1, 1, NPAR, 0, 0, ZERO,ZERO)
2030	GO TO 700
2031	629 CALL XERR('ODESSA - ISTATE = 3 AND NPAR CHANGED (I1 TO I2)',
2032	1 29, 1, 2, NP, NPAR, 0, ZERO,ZERO)
2033	GO TO 700
2034	630 CALL XERR('ODESSA - MITER (=I1) ILLEGAL',
2035	1 30, 1, 1, MITER, 0, 0, ZERO,ZERO)
2036	GO TO 700
2037	631 CALL XERR('ODESSA - TROUBLE IN SPRIME (IERPJ)',
2038	1 31, 1, 0, 0, 0, 0, ZERO,ZERO)
2039	GO TO 700
2040	632 CALL XERR('ODESSA - TROUBLE IN SPRIME (MITER)',
2041	1 32, 1, 0, 0, 0, 0, ZERO,ZERO)
2042	GO TO 700
2043	633 CALL XERR('ODESSA - FATAL ERROR IN STODE (KFLAG = -3)',
2044	1 33, 2, 0, 0, 0, 0, ZERO,ZERO)
2045	GO TO 801
2046	C
2047	700 IF (ILLIN .EQ. 5) GO TO 710
2048	ILLIN = ILLIN + 1
2049	ISTATE = -3
2050	RETURN
2051	710 CALL XERR('ODESSA - REPEATED OCCURRENCES OF ILLEGAL INPUT',
2052	1 302, 1, 0, 0, 0, 0, ZERO,ZERO)
2053	C
2054	800 CALL XERR('ODESSA - RUN ABORTED.. APPARENT INFINITE LOOP',
2055	1 303, 2, 0, 0, 0, 0, ZERO,ZERO)
2056	RETURN
2057	801 CALL XERR('ODESSA - RUN ABORTED',
2058	1 304, 2, 0, 0, 0, 0, ZERO,ZERO)
2059	RETURN
2060	C-------------------- END OF SUBROUTINE ODESSA -------------------------
2061	END
2062	DOUBLE PRECISION FUNCTION ADDX(A,B)
2063	DOUBLE PRECISION A,B
2064	C
2065	C THIS FUNCTION IS NECESSARY TO FORCE OPTIMIZING COMPILERS TO
2066	C EXECUTE AND STORE A SUM, FOR SUCCESSFUL EXECUTION OF THE
2067	C TEST A + B = B.
2068	C
2069	ADDX = A + B
2070	RETURN
2071	C-------------------- END OF FUNCTION SUM ------------------------------
2072	END
2073	SUBROUTINE SPRIME (NEQ, Y, YH, NYH, NROW, NCOL, WM, IWM,
2074	1 EWT, SAVF, FTEM, DFDP, PAR, F, JAC, DF, PJAC, PDF)
2075	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2076	DIMENSION NEQ(), Y(), YH(NROW,NCOL,), WM(), IWM(*),
2077	1 EWT(), SAVF(), FTEM(), DFDP(NROW,), PAR(*)
2078	EXTERNAL F, JAC, DF, PJAC, PDF
2079	PARAMETER (ONE=1.0D0,ZERO=0.0D0)
2080	COMMON /ODE001/ ROWND, ROWNS(173),
2081	1 RDUM1(37),EL0, H, RDUM2(6),
2082	2 IOWND1(14), IOWNS(4),
2083	3 IDUM1(3), IERPJ, IDUM2(6),
2084	4 MITER, IDUM3(4), N, IDUM4(5)
2085	COMMON /ODE002/ RDUM3(3),
2086	1 IOWND2(3), IDUM5, NSV, IDUM6, NSPE, IDUM7, IERSP, JOPT, IDUM8
2087	C-----------------------------------------------------------------------
2088	C SPRIME IS CALLED BY ODESSA TO INITIALIZE THE YH ARRAY. IT IS ALSO
2089	C CALLED BY STODE TO REEVALUATE FIRST ORDER DERIVATIVES WHEN KFLAG
2090	C .LE. -3. SPRIME COMPUTES THE FIRST DERIVATIVES OF THE SENSITIVITY
2091	C COEFFICIENTS WITH RESPECT TO THE INDEPENDENT VARIABLE T...
2092	C
2093	C SPRIME = D(DY/DP)/DT = JAC*DY/DP + DF/DP
2094	C WHERE JAC = JACOBIAN MATRIX
2095	C DY/DP = SENSITIVITY MATRIX
2096	C DF/DP = INHOMOGENEITY MATRIX
2097	C THIS ROUTINE USES THE COMMON VARIABLES EL0, H, IERPJ, MITER, N,
2098	C NSV, NSPE, IERSP, JOPT
2099	C-----------------------------------------------------------------------
2100	C CALL PREPJ WITH JOPT = 1.
2101	C IF MITER = 2 OR 5, EL0 IS TEMPORARILY SET TO -1.0 AND H IS
2102	C TEMPORARILY SET TO 1.0D0.
2103	C-----------------------------------------------------------------------
2104	NSPE = NSPE + 1
2105	JOPT = 1
2106	IF (MITER .EQ. 1 .OR. MITER .EQ. 4) GO TO 10
2107	HTEMP = H
2108	ETEMP = EL0
2109	H = ONE
2110	EL0 = -ONE
2111	10 CALL PJAC (NEQ, Y, YH, NYH, WM, IWM, EWT, SAVF, FTEM,
2112	1 PAR, F, JAC, JOPT)
2113	IF (IERPJ .NE. 0) GO TO 300
2114	JOPT = 0
2115	IF (MITER .EQ. 1 .OR. MITER .EQ. 4) GO TO 20
2116	H = HTEMP
2117	EL0 = ETEMP
2118	C-----------------------------------------------------------------------
2119	C CALL PREPDF AND LOAD DFDP(*,JPAR).
2120	C-----------------------------------------------------------------------
2121	20 DO 30 J = 2,NSV
2122	JPAR = J - 1
2123	CALL PDF (NEQ, Y, WM, SAVF, FTEM, DFDP(1,JPAR), PAR,
2124	1 F, DF, JPAR)
2125	30 CONTINUE
2126	C-----------------------------------------------------------------------
2127	C COMPUTE JACDY/DP AND STORE RESULTS IN YH(,*,2).
2128	C-----------------------------------------------------------------------
2129	GO TO (40,40,310,100,100) MITER
2130	C THE JACOBIAN IS FULL.------------------------------------------------
2131	C FOR EACH ROW OF THE JACOBIAN..
2132	C 40 DO 70 IROW = 1,N
2133	C AND EACH COLUMN OF THE SENSITIVITY MATRIX..
2134	C DO 60 J = 2,NSV
2135	C SUM = ZERO
2136	C TAKE THE VECTOR DOT PRODUCT..
2137	C DO 50 I = 1,N
2138	C IPD = IROW + N*(I-1) + 2
2139	C SUM = SUM + WM(IPD)*YH(I,J,1)
2140	C 50 CONTINUE
2141	C YH(IROW,J,2) = SUM
2142	C 60 CONTINUE
2143	C 70 CONTINUE
2144	40 CONTINUE
2145	C FOR EACH COLUMN OF THE SENSITIVITY MATRIX..
2146	DO 60 J = 2,NSV
2147	CALL Jac_SP_Vec( WM(3), YH(1,J,1), YH(1,J,2) )
2148	60 CONTINUE
2149	GO TO 200
2150	C THE JACOBIAN IS BANDED.-----------------------------------------------
2151	100 ML = IWM(1)
2152	MU = IWM(2)
2153	ICOUNT = 1
2154	MBAND = ML + MU + 1
2155	MEBAND = MBAND + ML
2156	NMU = N - MU
2157	ML1 = ML + 1
2158	C FOR EACH ROW OF THE JACOBIAN..
2159	DO 160 IROW = 1,N
2160	IF (IROW .GT. ML1) GO TO 110
2161	IPD = MBAND + IROW + 1
2162	IYH = 1
2163	LBAND = MU + IROW
2164	GO TO 120
2165	110 ICOUNT = ICOUNT + 1
2166	IPD = ICOUNT*MEBAND + 2
2167	IYH = IYH + 1
2168	LBAND = LBAND - 1
2169	IF (IROW .LE. NMU) LBAND = MBAND
2170	C AND EACH COLUMN OF THE SENSITIVITY MATRIX..
2171	120 DO 150 J = 2,NSV
2172	SUM = ZERO
2173	I1 = IPD
2174	I2 = IYH
2175	C TAKE THE VECTOR DOT PRODUCT.
2176	DO 140 I = 1,LBAND
2177	SUM = SUM + WM(I1)*YH(I2,J,1)
2178	I1 = I1 + MEBAND - 1
2179	I2 = I2 + 1
2180	140 CONTINUE
2181	YH(IROW,J,2) = SUM
2182	150 CONTINUE
2183	160 CONTINUE
2184	C-----------------------------------------------------------------------
2185	C ADD THE INHOMOGENEITY TERM, I.E., ADD DFDP(,JPAR) TO YH(,JPAR+1,2).
2186	C-----------------------------------------------------------------------
2187	200 DO 220 J = 2,NSV
2188	JPAR = J - 1
2189	DO 210 I = 1,N
2190	YH(I,J,2) = YH(I,J,2) + DFDP(I,JPAR)
2191	210 CONTINUE
2192	220 CONTINUE
2193	RETURN
2194	C-----------------------------------------------------------------------
2195	C ERROR RETURNS.
2196	C-----------------------------------------------------------------------
2197	300 IERSP = -1
2198	RETURN
2199	310 IERSP = -2
2200	RETURN
2201	C------------------------END OF SUBROUTINE SPRIME-----------------------
2202	END
2203	SUBROUTINE PREPJ (NEQ, Y, YH, NYH, WM, IWM, EWT, SAVF, FTEM,
2204	1 PAR, FUNC_CHEM, JAC, JOPT)
2205	C IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2206	INCLUDE 'KPP_ROOT_Parameters.h'
2207	INCLUDE 'KPP_ROOT_Sparse.h'
2208	DIMENSION NEQ(), Y(), YH(NYH,), WM(), IWM(), EWT(),
2209	1 SAVF(), FTEM(), PAR(*)
2210	EXTERNAL FUNC_CHEM, JAC
2211	PARAMETER (ZERO=0.0D0,ONE=1.0D0)
2212	COMMON /ODE001/ ROWND, ROWNS(173),
2213	2 RDUM1(37), EL0, H, RDUM2(4), TN, UROUND,
2214	3 IOWND(14), IOWNS(4),
2215	4 IDUM1(3), IERPJ, IDUM2, JCUR, IDUM3(4),
2216	5 MITER, IDUM4(4), N, IDUM5(2), NFE, NJE, IDUM6
2217	C-----------------------------------------------------------------------
2218	C PREPJ IS CALLED BY STODE TO COMPUTE AND PROCESS THE MATRIX
2219	C P = I - HEL(1)J , WHERE J IS AN APPROXIMATION TO THE JACOBIAN.
2220	C IF ISOPT = 1, PREPJ IS ALSO CALLED BY SPRIME WITH JOPT = 1.
2221	C HERE J IS COMPUTED BY THE USER-SUPPLIED ROUTINE JAC IF
2222	C MITER = 1 OR 4, OR BY FINITE DIFFERENCING IF MITER = 2, 3, OR 5.
2223	C IF MITER = 3, A DIAGONAL APPROXIMATION TO J IS USED.
2224	C J IS STORED IN WM AND REPLACED BY P. IF MITER .NE. 3, P IS THEN
2225	C SUBJECTED TO LU DECOMPOSITION (JOPT = 0) IN PREPARATION FOR LATER
2226	C SOLUTION OF LINEAR SYSTEMS WITH P AS COEFFICIENT MATRIX. THIS IS
2227	C DONE BY DGEFA IF MITER = 1 OR 2, AND BY DGBFA IF MITER = 4 OR 5.
2228	C
2229	C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
2230	C WITH PREPJ USES THE FOLLOWING..
2231	C Y = ARRAY CONTAINING PREDICTED VALUES ON ENTRY.
2232	C FTEM = WORK ARRAY OF LENGTH N (ACOR IN STODE).
2233	C SAVF = ARRAY CONTAINING F EVALUATED AT PREDICTED Y.
2234	C WM = REAL WORK SPACE FOR MATRICES. ON OUTPUT IT CONTAINS THE
2235	C INVERSE DIAGONAL MATRIX IF MITER = 3 AND THE LU DECOMPOSITION
2236	C OF P IF MITER IS 1, 2 , 4, OR 5.
2237	C STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
2238	C WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
2239	C WM(1) = SQRT(UROUND), USED IN NUMERICAL JACOBIAN INCREMENTS.
2240	C WM(2) = H*EL0, SAVED FOR LATER USE IF MITER = 3.
2241	C IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT
2242	C IWM(21), IF MITER IS 1, 2, 4, OR 5. IWM ALSO CONTAINS BAND
2243	C PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5.
2244	C EL0 = EL(1) (INPUT).
2245	C IERPJ = OUTPUT ERROR FLAG, = 0 IF NO TROUBLE, .GT. 0 IF
2246	C P MATRIX FOUND TO BE SINGULAR.
2247	C JCUR = OUTPUT FLAG = 1 TO INDICATE THAT THE JACOBIAN MATRIX
2248	C (OR APPROXIMATION) IS NOW CURRENT.
2249	C JOPT = INPUT JACOBIAN OPTION, = 1 IF JAC IS DESIRED ONLY.
2250	C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, UROUND,
2251	C IERPJ, JCUR, MITER, N, NFE, AND NJE.
2252	C-----------------------------------------------------------------------
2253	NJE = NJE + 1
2254	IERPJ = 0
2255	JCUR = 1
2256	HL0 = H*EL0
2257	GO TO (100, 200, 300, 400, 500), MITER
2258	C IF MITER = 1, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
2259	C 100 LENP = N*N
2260	100 LENP = LU_NONZERO
2261	DO 110 I = 1,LU_NONZERO
2262	110 WM(I+2) = ZERO
2263	CALL JAC (NEQ, TN, Y, PAR, 0, 0, WM(3), N)
2264	IF (JOPT .EQ. 1) RETURN
2265	CON = -HL0
2266	DO 120 I = 1,LU_NONZERO
2267	120 WM(I+2) = WM(I+2)*CON
2268	GO TO 240
2269	C IF MITER = 2, MAKE N CALLS TO F TO APPROXIMATE J. --------------------
2270	200 FAC = VNORM (N, SAVF, EWT)
2271	R0 = 1000.0D0ABS(H)UROUNDREAL(N)FAC
2272	IF (R0 .EQ. ZERO) R0 = ONE
2273	SRUR = WM(1)
2274	J1 = 2
2275	DO 230 J = 1,N
2276	YJ = Y(J)
2277	R = MAX(SRUR*ABS(YJ),R0/EWT(J))
2278	Y(J) = Y(J) + R
2279	FAC = -HL0/R
2280	CALL FUNC_CHEM (NEQ, TN, Y, PAR, FTEM)
2281	DO 220 I = 1,N
2282	220 WM(I+J1) = (FTEM(I) - SAVF(I))*FAC
2283	Y(J) = YJ
2284	J1 = J1 + N
2285	230 CONTINUE
2286	NFE = NFE + N
2287	IF (JOPT .EQ. 1) RETURN
2288	C ADD IDENTITY MATRIX. -------------------------------------------------
2289	240 J = 3
2290	C DO 250 I = 1,N
2291	C WM(J) = WM(J) + ONE
2292	C 250 J = J + (N + 1)
2293	DO 250 I = 1,NVAR
2294	250 WM(2+LU_DIAG(I)) = WM(2+LU_DIAG(I)) + ONE
2295	C DO LU DECOMPOSITION ON P. --------------------------------------------
2296	C CALL DGEFA (WM(3), N, N, IWM(21), IER)
2297	CALL KppDecomp (WM(3), IER)
2298	IF (IER .NE. 0) THEN
2299	IERPJ = 1
2300	PRINT*,"Singular Matrix"
2301	STOP
2302	END IF
2303	RETURN
2304	C IF MITER = 3, CONSTRUCT A DIAGONAL APPROXIMATION TO J AND P. ---------
2305	300 WM(2) = HL0
2306	R = EL0*0.1D0
2307	DO 310 I = 1,N
2308	310 Y(I) = Y(I) + R(HSAVF(I) - YH(I,2))
2309	CALL FUNC_CHEM (NEQ, TN, Y, PAR, WM(3))
2310	NFE = NFE + 1
2311	DO 320 I = 1,N
2312	R0 = H*SAVF(I) - YH(I,2)
2313	DI = 0.1D0R0 - H(WM(I+2) - SAVF(I))
2314	WM(I+2) = 1.0D0
2315	IF (ABS(R0) .LT. UROUND/EWT(I)) GO TO 320
2316	IF (ABS(DI) .EQ. ZERO) GO TO 330
2317	WM(I+2) = 0.1D0*R0/DI
2318	320 CONTINUE
2319	RETURN
2320	330 IERPJ = 1
2321	RETURN
2322	C IF MITER = 4, CALL JAC AND MULTIPLY BY SCALAR. -----------------------
2323	400 ML = IWM(1)
2324	MU = IWM(2)
2325	ML3 = ML + 3
2326	MBAND = ML + MU + 1
2327	MEBAND = MBAND + ML
2328	LENP = MEBAND*N
2329	DO 410 I = 1,LENP
2330	410 WM(I+2) = ZERO
2331	CALL JAC (NEQ, TN, Y, PAR, ML, MU, WM(ML3), MEBAND)
2332	IF (JOPT .EQ. 1) RETURN
2333	CON = -HL0
2334	DO 420 I = 1,LENP
2335	420 WM(I+2) = WM(I+2)*CON
2336	GO TO 570
2337	C IF MITER = 5, MAKE MBAND CALLS TO F TO APPROXIMATE J. ----------------
2338	500 ML = IWM(1)
2339	MU = IWM(2)
2340	MBAND = ML + MU + 1
2341	MBA = MIN(MBAND,N)
2342	MEBAND = MBAND + ML
2343	MEB1 = MEBAND - 1
2344	SRUR = WM(1)
2345	FAC = VNORM (N, SAVF, EWT)
2346	R0 = 1000.0D0ABS(H)UROUNDREAL(N)FAC
2347	IF (R0 .EQ. ZERO) R0 = ONE
2348	DO 560 J = 1,MBA
2349	DO 530 I = J,N,MBAND
2350	YI = Y(I)
2351	R = MAX(SRUR*ABS(YI),R0/EWT(I))
2352	530 Y(I) = Y(I) + R
2353	CALL FUNC_CHEM (NEQ, TN, Y, PAR, FTEM)
2354	DO 550 JJ = J,N,MBAND
2355	Y(JJ) = YH(JJ,1)
2356	YJJ = Y(JJ)
2357	R = MAX(SRUR*ABS(YJJ),R0/EWT(JJ))
2358	FAC = -HL0/R
2359	I1 = MAX(JJ-MU,1)
2360	I2 = MIN(JJ+ML,N)
2361	II = JJ*MEB1 - ML + 2
2362	DO 540 I = I1,I2
2363	540 WM(II+I) = (FTEM(I) - SAVF(I))*FAC
2364	550 CONTINUE
2365	560 CONTINUE
2366	NFE = NFE + MBA
2367	IF (JOPT .EQ. 1) RETURN
2368	C ADD IDENTITY MATRIX. -------------------------------------------------
2369	570 II = MBAND + 2
2370	DO 580 I = 1,N
2371	WM(II) = WM(II) + ONE
2372	580 II = II + MEBAND
2373	C DO LU DECOMPOSITION OF P. --------------------------------------------
2374	CALL DGBFA (WM(3), MEBAND, N, ML, MU, IWM(21), IER)
2375	IF (IER .NE. 0) IERPJ = 1
2376	RETURN
2377	C----------------------- END OF SUBROUTINE PREPJ -----------------------
2378	END
2379	SUBROUTINE PREPDF (NEQ, Y, SRUR, SAVF, FTEM, DFDP, PAR,
2380	1 F, DF, JPAR)
2381	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2382	EXTERNAL F, DF
2383	DIMENSION NEQ(), Y(), SAVF(), FTEM(), DFDP(), PAR()
2384	COMMON /ODE001/ ROWND, ROWNS(173),
2385	1 RDUM1(43), TN, RDUM2,
2386	2 IOWND1(14), IOWNS(4),
2387	3 IDUM1(10), MITER, IDUM2(4), N, IDUM3(2), NFE, IDUM4(2)
2388	COMMON /ODE002/ RDUM3(3),
2389	1 IOWND2(3), IDUM5(2), NDFE, IDUM6, IDF, IDUM7(3)
2390	C-----------------------------------------------------------------------
2391	C PREPDF IS CALLED BY SPRIME AND STESA TO COMPUTE THE INHOMOGENEITY
2392	C VECTORS DF(I)/DP(JPAR). HERE DF/DP IS COMPUTED BY THE USER-SUPPLIED
2393	C ROUTINE DF IF IDF = 1, OR BY FINITE DIFFERENCING IF IDF = 0.
2394	C
2395	C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION WITH
2396	C PREPDF USES THE FOLLOWING..
2397	C Y = REAL ARRAY OF LENGTH NYH CONTAINING DEPENDENT VARIABLES.
2398	C PREPDF USES ONLY THE FIRST N ENTRIES OF Y(*).
2399	C SRUR = SQRT(UROUND) (= WM(1)).
2400	C SAVF = REAL ARRAY OF LENGTH N CONTAINING DERIVATIVES DY/DT.
2401	C FTEM = REAL ARRAY OF LENGTH N USED TO TEMPORARILY STORE DY/DT FOR
2402	C NUMERICAL DIFFERENTIATION.
2403	C DFDP = REAL ARRAY OF LENGTH N USED TO STORE DF(I)/DP(JPAR), I = 1,N.
2404	C PAR = REAL ARRAY OF LENGTH NPAR CONTAINING EQUATION PARAMETERS
2405	C OF INTEREST.
2406	C JPAR = INPUT PARAMETER, 2 .LE. JPAR .LE. NSV, DESIGNATING THE
2407	C APPROPRIATE SOLUTION VECTOR CORRESPONDING TO PAR(JPAR).
2408	C THIS ROUTINE ALSO USES THE COMMON VARIABLES TN, MITER, N, NFE, NDFE,
2409	C AND IDF.
2410	C-----------------------------------------------------------------------
2411	NDFE = NDFE + 1
2412	IDF1 = IDF + 1
2413	GO TO (100, 200), IDF1
2414	C IDF = 0, CALL F TO APPROXIMATE DFDP. ---------------------------------
2415	100 RPAR = PAR(JPAR)
2416	R = MAX(SRUR*ABS(RPAR),SRUR)
2417	PAR(JPAR) = RPAR + R
2418	FAC = 1.0D0/R
2419	CALL F (NEQ, TN, Y, PAR, FTEM)
2420	DO 110 I = 1,N
2421	110 DFDP(I) = (FTEM(I) - SAVF(I))*FAC
2422	PAR(JPAR) = RPAR
2423	NFE = NFE + 1
2424	RETURN
2425	C IDF = 1, CALL USER SUPPLIED DF. --------------------------------------
2426	200 DO 210 I = 1,N
2427	210 DFDP(I) = 0.0D0
2428	CALL DF (NEQ, TN, Y, PAR, DFDP, JPAR)
2429	RETURN
2430	C -------------------- END OF SUBROUTINE PREPDF ------------------------
2431	END
2432	SUBROUTINE INTDY (T, K, YH, NYH, DKY, IFLAG)
2433	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2434	DIMENSION YH(NYH,1), DKY(1)
2435	COMMON /ODE001/ ROWND, ROWNS(173),
2436	2 RDUM1(38),H, RDUM2(2), HU, RDUM3, TN, UROUND,
2437	3 IOWND(14), IOWNS(4),
2438	4 IDUM1(8), L, IDUM2,
2439	5 IDUM3(5), N, NQ, IDUM4(4)
2440	C-----------------------------------------------------------------------
2441	C INTDY COMPUTES INTERPOLATED VALUES OF THE K-TH DERIVATIVE OF THE
2442	C DEPENDENT VARIABLE VECTOR Y, AND STORES IT IN DKY. THIS ROUTINE
2443	C IS CALLED WITHIN THE PACKAGE WITH K = 0 AND T = TOUT, BUT MAY
2444	C ALSO BE CALLED BY THE USER FOR ANY K UP TO THE CURRENT ORDER.
2445	C (SEE DETAILED INSTRUCTIONS IN THE USAGE DOCUMENTATION.)
2446	C-----------------------------------------------------------------------
2447	C THE COMPUTED VALUES IN DKY ARE GOTTEN BY INTERPOLATION USING THE
2448	C NORDSIECK HISTORY ARRAY YH. THIS ARRAY CORRESPONDS UNIQUELY TO A
2449	C VECTOR-VALUED POLYNOMIAL OF DEGREE NQCUR OR LESS, AND DKY IS SET
2450	C TO THE K-TH DERIVATIVE OF THIS POLYNOMIAL AT T.
2451	C THE FORMULA FOR DKY IS..
2452	C Q
2453	C DKY(I) = SUM C(J,K) * (T - TN)*(J-K) H*(-J) YH(I,J+1)
2454	C J=K
2455	C WHERE C(J,K) = J(J-1)...*(J-K+1), Q = NQCUR, TN = TCUR, H = HCUR.
2456	C THE QUANTITIES NQ = NQCUR, L = NQ+1, N = NEQ, TN, AND H ARE
2457	C COMMUNICATED BY COMMON. THE ABOVE SUM IS DONE IN REVERSE ORDER.
2458	C IFLAG IS RETURNED NEGATIVE IF EITHER K OR T IS OUT OF BOUNDS.
2459	C-----------------------------------------------------------------------
2460	IFLAG = 0
2461	IF (K .LT. 0 .OR. K .GT. NQ) GO TO 80
2462	TP = TN - HU(1.0D0 + 100.0D0UROUND)
2463	IF ((T-TP)*(T-TN) .GT. 0.0D0) GO TO 90
2464	C
2465	S = (T - TN)/H
2466	IC = 1
2467	IF (K .EQ. 0) GO TO 15
2468	JJ1 = L - K
2469	DO 10 JJ = JJ1,NQ
2470	10 IC = IC*JJ
2471	15 C = REAL(IC)
2472	DO 20 I = 1,NYH
2473	20 DKY(I) = C*YH(I,L)
2474	IF (K .EQ. NQ) GO TO 55
2475	JB2 = NQ - K
2476	DO 50 JB = 1,JB2
2477	J = NQ - JB
2478	JP1 = J + 1
2479	IC = 1
2480	IF (K .EQ. 0) GO TO 35
2481	JJ1 = JP1 - K
2482	DO 30 JJ = JJ1,J
2483	30 IC = IC*JJ
2484	35 C = REAL(IC)
2485	DO 40 I = 1,NYH
2486	40 DKY(I) = CYH(I,JP1) + SDKY(I)
2487	50 CONTINUE
2488	IF (K .EQ. 0) RETURN
2489	55 R = H**(-K)
2490	DO 60 I = 1,NYH
2491	60 DKY(I) = R*DKY(I)
2492	RETURN
2493	C
2494	80 CALL XERR('INTDY-- K (=I1) ILLEGAL',
2495	1 51, 1, 1, K, 0, 0, ZERO,ZERO)
2496	IFLAG = -1
2497	RETURN
2498	90 CALL XERR ('INTDY-- T (=R1) ILLEGAL',
2499	1 52, 1, 0, 0, 0, 1, T, ZERO)
2500	CALL XERR('T NOT IN INTERVAL TCUR - HU (= R1) TO TCUR (=R2)',
2501	1 52, 1, 0, 0, 0, 2, TP, TN)
2502	IFLAG = -2
2503	RETURN
2504	C----------------------- END OF SUBROUTINE INTDY -----------------------
2505	END
2506	SUBROUTINE STESA (NEQ, Y, NROW, NCOL, YH, WM, IWM, EWT, SAVF,
2507	1 ACOR, PAR, NRS, F, JAC, DF, PJAC, PDF, KppSolve)
2508	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2509	EXTERNAL F, JAC, DF, PJAC, PDF, KppSolve
2510	DIMENSION NEQ(), Y(NROW,), YH(NROW,NCOL,), WM(), IWM(*),
2511	1 EWT(NROW,), SAVF(), ACOR(NROW,), PAR(), NRS(*)
2512	PARAMETER (ONE=1.0D0,ZERO=0.0D0)
2513	COMMON /ODE001/ ROWND, ROWNS(173),
2514	1 TESCO(3,12), RDUM1, EL0, H, RDUM2(4), TN, RDUM3,
2515	2 IOWND1(14), IOWNS(4),
2516	3 IALTH, LMAX, IDUM1, IERPJ, IERSL, JCUR, IDUM2, KFLAG, L, IDUM3,
2517	4 MITER, IDUM4(4), N, NQ, IDUM5, NFE, IDUM6(2)
2518	COMMON /ODE002/ DUPS, DSMS, DDNS,
2519	1 IOWND2(3), IDUM7, NSV, IDUM8(2), IDF, IDUM9, JOPT, KFLAGS
2520	C-----------------------------------------------------------------------
2521	C STESA IS CALLED BY STODE TO PERFORM AN EXPLICIT CALCULATION FOR THE
2522	C FIRST-ORDER SENSITIVITY COEFFICIENTS DY(I)/DP(J), I = 1,N; J = 1,NPAR.
2523	C
2524	C IN ADDITION TO THE VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
2525	C WITH STESA USES THE FOLLOWING..
2526	C Y = AN NROW (=N) BY NCOL (=NSV) REAL ARRAY CONTAINING THE
2527	C CORRECTED DEPENDENT VARIABLES ON OUTPUT..
2528	C Y(I,1) , I = 1,N = STATE VARIABLES (INPUT);
2529	C Y(I,J) , I = 1,N , J = 2,NSV ,
2530	C = SENSITIVITY COEFFICIENTS, DY(I)/DP(J).
2531	C YH = AN N BY NSV BY LMAX REAL ARRAY CONTAINING THE PREDICTED
2532	C DEPENDENT VARIABLES AND THEIR APPROXIMATE SCALED DERIVATIVES.
2533	C SAVF = A REAL ARRAY OF LENGTH N USED TO STORE FIRST DERIVATIVES
2534	C OF DEPENDENT VARIABLES IF MITER = 2 OR 5.
2535	C PAR = A REAL ARRAY OF LENGTH NPAR CONTAINING THE EQUATION
2536	C PARAMETERS OF INTEREST.
2537	C NRS = AN INTEGER ARRAY OF LENGTH NPAR + 1 CONTAINING THE NUMBER
2538	C OF REPEATED STEPS (KFLAGS .LT. 0) DUE TO THE SENSITIVITY
2539	C CALCULATIONS..
2540	C NRS(1) = TOTAL NUMBER OF REPEATED STEPS
2541	C NRS(I) , I = 2,NPAR = NUMBER OF REPEATED STEPS DUE
2542	C TO PARAMETER I.
2543	C NSV = NUMBER OF SOLUTION VECTORS = NPAR + 1.
2544	C KFLAGS = LOCAL ERROR TEST FLAG, = 0 IF TEST PASSES, .LT. 0 IF TEST
2545	C FAILS, AND STEP NEEDS TO BE REPEATED. ERROR TEST IS APPLIED
2546	C TO EACH SOLUTION VECTOR INDEPENDENTLY.
2547	C DUPS, DSMS, DDNS = REAL SCALARS USED FOR COMPUTING RHUP, RHSM, RHDN,
2548	C ON RETURN TO STODE (IALTH .EQ. 1).
2549	C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, TN, IALTH, LMAX,
2550	C IERPJ, IERSL, JCUR, KFLAG, L, MITER, N, NQ, NFE, AND JOPT.
2551	C-----------------------------------------------------------------------
2552	DUPS = ZERO
2553	DSMS = ZERO
2554	DDNS = ZERO
2555	HL0 = H*EL0
2556	EL0I = ONE/EL0
2557	TI2 = ONE/TESCO(2,NQ)
2558	TI3 = ONE/TESCO(3,NQ)
2559	C IF MITER = 2 OR 5 (OR IDF = 0), SUPPLY DERIVATIVES AT CORRECTED
2560	C Y(*,1) VALUES FOR NUMERICAL DIFFERENTIATION IN PJAC AND/OR PDF.
2561	IF (MITER .EQ. 2 .OR. MITER .EQ. 5 .OR. IDF .EQ. 0) GO TO 10
2562	GO TO 15
2563	10 CALL F (NEQ, TN, Y, PAR, SAVF)
2564	NFE = NFE + 1
2565	C IF JCUR = 0, UPDATE THE JACOBIAN MATRIX.
2566	C IF MITER = 5, LOAD CORRECTED Y(,1) VALUES INTO Y(,2).
2567	15 IF (JCUR .EQ. 1) GO TO 30
2568	IF (MITER .NE. 5) GO TO 25
2569	DO 20 I = 1,N
2570	20 Y(I,2) = Y(I,1)
2571	25 CALL PJAC (NEQ, Y, Y(1,2), N, WM, IWM, EWT, SAVF, ACOR(1,2),
2572	1 PAR, F, JAC, JOPT)
2573	IF (IERPJ .NE. 0) RETURN
2574	C-----------------------------------------------------------------------
2575	C THIS IS A LOOPING POINT FOR THE SENSITIVITY CALCULATIONS.
2576	C-----------------------------------------------------------------------
2577	C FOR EACH PARAMETER PAR(*), A SENSITIVITY SOLUTION VECTOR IS COMPUTED
2578	C USING THE SAME STEP SIZE (H) AND ORDER (NQ) AS IN STODE.
2579	C A LOCAL ERROR TEST IS APPLIED INDEPENDENTLY TO EACH SOLUTION VECTOR.
2580	C-----------------------------------------------------------------------
2581	30 DO 100 J = 2,NSV
2582	JPAR = J - 1
2583	C EVALUATE INHOMOGENEITY TERM, TEMPORARILY LOAD INTO Y(*,JPAR+1). ------
2584	CALL PDF(NEQ, Y, WM, SAVF, ACOR(1,J), Y(1,J), PAR,
2585	1 F, DF, JPAR)
2586	C-----------------------------------------------------------------------
2587	C LOAD RHS OF SENSITIVITY SOLUTION (CORRECTOR) EQUATION..
2588	C
2589	C RHS = DY/DP - EL(1)HD(DY/DP)/DT + EL(1)HDF/DP
2590	C
2591	C-----------------------------------------------------------------------
2592	DO 40 I = 1,N
2593	40 Y(I,J) = YH(I,J,1) - EL0YH(I,J,2) + HL0Y(I,J)
2594	C-----------------------------------------------------------------------
2595	C KppSolve CORRECTOR EQUATION: THE SOLUTIONS ARE LOCATED IN Y(*,JPAR+1).
2596	C THE EXPLICIT FORMULA IS..
2597	C
2598	C (I - EL(1)HJAC) * DY/DP(CORRECTED) = RHS
2599	C
2600	C-----------------------------------------------------------------------
2601	CALL KppSolve (WM, IWM, Y(1,J), DUM)
2602	IF (IERSL .NE. 0) RETURN
2603	C ESTIMATE LOCAL TRUNCATION ERROR. -------------------------------------
2604	DO 50 I = 1,N
2605	50 ACOR(I,J) = (Y(I,J) - YH(I,J,1))*EL0I
2606	ERR = VNORM(N, ACOR(1,J), EWT(1,J))*TI2
2607	IF (ERR .GT. ONE) GO TO 200
2608	C-----------------------------------------------------------------------
2609	C LOCAL ERROR TEST PASSED. SET KFLAGS TO 0 TO INDICATE THIS.
2610	C IF IALTH = 1, COMPUTE DSMS, DDNS, AND DUPS (IF L .LT. LMAX).
2611	C-----------------------------------------------------------------------
2612	KFLAGS = 0
2613	IF (IALTH .GT. 1) GO TO 100
2614	IF (L .EQ. LMAX) GO TO 70
2615	DO 60 I= 1,N
2616	60 Y(I,J) = ACOR(I,J) - YH(I,J,LMAX)
2617	DUPS = MAX(DUPS,VNORM(N,Y(1,J),EWT(1,J))*TI3)
2618	70 DSMS = MAX(DSMS,ERR)
2619	100 CONTINUE
2620	RETURN
2621	C-----------------------------------------------------------------------
2622	C THIS SECTION IS REACHED IF THE ERROR TOLERANCE FOR SENSITIVITY
2623	C SOLUTION VECTOR JPAR HAS BEEN VIOLATED. KFLAGS IS MADE NEGATIVE TO
2624	C INDICATE THIS. IF KFLAGS = -1, SET KFLAG EQUAL TO ZERO SO THAT KFLAG
2625	C IS SET TO -1 ON RETURN TO STODE BEFORE REPEATING THE STEP.
2626	C INCREMENT NRS(1) (= TOTAL NUMBER OF REPEATED STEPS DUE TO ALL
2627	C SENSITIVITY SOLUTION VECTORS) BY ONE.
2628	C INCREMENT NRS(JPAR+1) (= TOTAL NUMBER OF REPEATED STEPS DUE TO
2629	C SOLUTION VECTOR JPAR+1) BY ONE.
2630	C LOAD DSMS FOR RH CALCULATION IN STODE.
2631	C-----------------------------------------------------------------------
2632	200 KFLAGS = KFLAGS - 1
2633	IF (KFLAGS .EQ. -1) KFLAG = 0
2634	NRS(1) = NRS(1) + 1
2635	NRS(J) = NRS(J) + 1
2636	DSMS = ERR
2637	RETURN
2638	C------------------------ END OF SUBROUTINE STESA ----------------------
2639	END
2640	SUBROUTINE STODE (NEQ, Y, YH, NYH, YH1, WM, IWM, EWT, SAVF, ACOR,
2641	1 PAR, NRS, F, JAC, DF, PJAC, PDF, SLVS)
2642	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2643	EXTERNAL F, JAC, DF, PJAC, PDF, SLVS
2644	DIMENSION NEQ(), Y(), YH(NYH,), YH1(), WM(), IWM(), EWT(*),
2645	1 SAVF(), ACOR(), PAR(), NRS()
2646	PARAMETER (ONE=1.0D0,ZERO=0.0D0)
2647	COMMON /ODE001/ ROWND,
2648	1 CONIT, CRATE, EL(13), ELCO(13,12), HOLD, RMAX,
2649	2 TESCO(3,12), CCMAX, EL0, H, HMIN, HMXI, HU, RC, TN, UROUND,
2650	3 IOWND1(14), IPUP, MEO, NQNYH, NSLP,
2651	4 IALTH, LMAX, ICF, IERPJ, IERSL, JCUR, JSTART, KFLAG, L, METH,
2652	5 MITER, MAXORD, MAXCOR, MSBP, MXNCF, N, NQ, NST, NFE, NJE, NQU
2653	COMMON /ODE002/ DUPS, DSMS, DDNS,
2654	1 IOWND2(3), ISOPT, NSV, NDFE, NSPE, IDF, IERSP, JOPT, KFLAGS
2655	C-----------------------------------------------------------------------
2656	C STODE PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE
2657	C PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS.
2658	C NOTE.. STODE IS INDEPENDENT OF THE VALUE OF THE ITERATION METHOD
2659	C INDICATOR MITER, WHEN THIS IS .NE. 0, AND HENCE IS INDEPENDENT
2660	C OF THE TYPE OF CHORD METHOD USED, OR THE JACOBIAN STRUCTURE.
2661	C FOR ISOPT = 1, STODE CALLS STESA FOR SENSITIVITY CALCULATIONS.
2662	C VARIABLES USED FOR COMMUNICATION WITH STESA ARE DESCRIBED IN STESA.
2663	C COMMUNICATION WITH STODE IS DONE WITH THE FOLLOWING VARIABLES..
2664	C
2665	C NEQ = INTEGER ARRAY CONTAINING PROBLEM SIZE IN NEQ(1), AND
2666	C NUMBER OF PARAMETERS TO BE CONSIDERED IN THE SENSITIVITY
2667	C ANALYSIS NEQ(2) (FOR ISOPT = 1), AND PASSED AS THE
2668	C NEQ ARGUMENT IN ALL CALLS TO F, JAC, AND DF.
2669	C Y = AN ARRAY OF LENGTH .GE. N USED AS THE Y ARGUMENT IN
2670	C ALL CALLS TO F, JAC, AND DF.
2671	C YH = AN NYH BY LMAX ARRAY CONTAINING THE DEPENDENT VARIABLES
2672	C AND THEIR APPROXIMATE SCALED DERIVATIVES, WHERE
2673	C LMAX = MAXORD + 1. YH(I,J+1) CONTAINS THE APPROXIMATE
2674	C J-TH DERIVATIVE OF Y(I), SCALED BY H**J/FACTORIAL(J)
2675	C (J = 0,1,...,NQ). ON ENTRY FOR THE FIRST STEP, THE FIRST
2676	C TWO COLUMNS OF YH MUST BE SET FROM THE INITIAL VALUES.
2677	C NYH = A CONSTANT INTEGER .GE. N, THE FIRST DIMENSION OF YH.
2678	C THE TOTAL NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS..
2679	C NYH = N, ISOPT = 0,
2680	C NYH = N * (NPAR + 1), ISOPT = 1
2681	C YH1 = A ONE-DIMENSIONAL ARRAY OCCUPYING THE SAME SPACE AS YH.
2682	C EWT = AN ARRAY OF LENGTH NYH CONTAINING MULTIPLICATIVE WEIGHTS
2683	C FOR LOCAL ERROR MEASUREMENTS. LOCAL ERRORS IN Y(I) ARE
2684	C COMPARED TO 1.0/EWT(I) IN VARIOUS ERROR TESTS.
2685	C SAVF = AN ARRAY OF WORKING STORAGE, OF LENGTH N.
2686	C ALSO USED FOR INPUT OF YH(*,MAXORD+2) WHEN JSTART = -1
2687	C AND MAXORD .LT. THE CURRENT ORDER NQ.
2688	C ACOR = A WORK ARRAY OF LENGTH NYH, USED FOR THE ACCUMULATED
2689	C CORRECTIONS. ON A SUCCESSFUL RETURN, ACOR(I) CONTAINS
2690	C THE ESTIMATED ONE-STEP LOCAL ERROR IN Y(I).
2691	C WM,IWM = REAL AND INTEGER WORK ARRAYS ASSOCIATED WITH MATRIX
2692	C OPERATIONS IN CHORD ITERATION (MITER .NE. 0).
2693	C PJAC = NAME OF ROUTINE TO EVALUATE AND PREPROCESS JACOBIAN MATRIX
2694	C AND P = I - HEL0JAC, IF A CHORD METHOD IS BEING USED.
2695	C IF ISOPT = 1, PJAC CAN BE CALLED TO CALCULATE JAC BY
2696	C SETTING JOPT = 1.
2697	C SLVS = NAME OF ROUTINE TO KppSolve LINEAR SYSTEM IN CHORD ITERATION.
2698	C CCMAX = MAXIMUM RELATIVE CHANGE IN H*EL0 BEFORE PJAC IS CALLED.
2699	C H = THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP.
2700	C H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING THE
2701	C PROBLEM. H CAN BE EITHER POSITIVE OR NEGATIVE, BUT ITS
2702	C SIGN MUST REMAIN CONSTANT THROUGHOUT THE PROBLEM.
2703	C HMIN = THE MINIMUM ABSOLUTE VALUE OF THE STEP SIZE H TO BE USED.
2704	C HMXI = INVERSE OF THE MAXIMUM ABSOLUTE VALUE OF H TO BE USED.
2705	C HMXI = 0.0 IS ALLOWED AND CORRESPONDS TO AN INFINITE HMAX.
2706	C HMIN AND HMXI MAY BE CHANGED AT ANY TIME, BUT WILL NOT
2707	C TAKE EFFECT UNTIL THE NEXT CHANGE OF H IS CONSIDERED.
2708	C TN = THE INDEPENDENT VARIABLE. TN IS UPDATED ON EACH STEP TAKEN.
2709	C JSTART = AN INTEGER USED FOR INPUT ONLY, WITH THE FOLLOWING
2710	C VALUES AND MEANINGS..
2711	C 0 PERFORM THE FIRST STEP.
2712	C .GT.0 TAKE A NEW STEP CONTINUING FROM THE LAST.
2713	C -1 TAKE THE NEXT STEP WITH A NEW VALUE OF H, MAXORD,
2714	C N, METH, OR MITER.
2715	C -2 TAKE THE NEXT STEP WITH A NEW VALUE OF H,
2716	C BUT WITH OTHER INPUTS UNCHANGED.
2717	C ON RETURN, JSTART IS SET TO 1 TO FACILITATE CONTINUATION.
2718	C KFLAG = A COMPLETION CODE WITH THE FOLLOWING MEANINGS..
2719	C 0 THE STEP WAS SUCCESFUL.
2720	C -1 THE REQUESTED ERROR COULD NOT BE ACHIEVED.
2721	C -2 CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED.
2722	C -3 FATAL ERROR IN PJAC, OR SLVS, (OR STESA).
2723	C A RETURN WITH KFLAG = -1 OR -2 MEANS EITHER
2724	C ABS(H) = HMIN OR 10 CONSECUTIVE FAILURES OCCURRED.
2725	C ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF TN AND
2726	C THE YH ARRAY ARE AS OF THE BEGINNING OF THE LAST
2727	C STEP, AND H IS THE LAST STEP SIZE ATTEMPTED.
2728	C MAXORD = THE MAXIMUM ORDER OF INTEGRATION METHOD TO BE ALLOWED.
2729	C MAXCOR = THE MAXIMUM NUMBER OF CORRECTOR ITERATIONS ALLOWED.
2730	C (= 3, IF ISOPT = 0)
2731	C (= 4, IF ISOPT = 1)
2732	C MSBP = MAXIMUM NUMBER OF STEPS BETWEEN PJAC CALLS (MITER .GT. 0).
2733	C IF ISOPT = 1, PJAC IS CALLED AT LEAST ONCE EVERY STEP.
2734	C MXNCF = MAXIMUM NUMBER OF CONVERGENCE FAILURES ALLOWED.
2735	C METH/MITER = THE METHOD FLAGS. SEE DESCRIPTION IN DRIVER.
2736	C N = THE NUMBER OF FIRST-ORDER MODEL DIFFERENTIAL EQUATIONS.
2737	C-----------------------------------------------------------------------
2738	KFLAG = 0
2739	KFLAGS = 0
2740	TOLD = TN
2741	NCF = 0
2742	IERPJ = 0
2743	IERSL = 0
2744	JCUR = 0
2745	ICF = 0
2746	IF (JSTART .GT. 0) GO TO 200
2747	IF (JSTART .EQ. -1) GO TO 100
2748	IF (JSTART .EQ. -2) GO TO 160
2749	C-----------------------------------------------------------------------
2750	C ON THE FIRST CALL, THE ORDER IS SET TO 1, AND OTHER VARIABLES ARE
2751	C INITIALIZED. RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE INCREASED
2752	C IN A SINGLE STEP. IT IS INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL
2753	C INITIAL H, BUT THEN IS NORMALLY EQUAL TO 10. IF A FAILURE
2754	C OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 2
2755	C FOR THE NEXT INCREASE.
2756	C THESE COMPUTATIONS CONSIDER ONLY THE ORIGINAL SOLUTION VECTOR.
2757	C THE SENSITIVITY SOLUTION VECTORS ARE CONSIDERED IN STESA (ISOPT = 1).
2758	C-----------------------------------------------------------------------
2759	LMAX = MAXORD + 1
2760	NQ = 1
2761	L = 2
2762	IALTH = 2
2763	RMAX = 10000.0D0
2764	RC = ZERO
2765	EL0 = ONE
2766	CRATE = 0.7D0
2767	DELP = ZERO
2768	HOLD = H
2769	MEO = METH
2770	NSLP = 0
2771	IPUP = MITER
2772	IRET = 3
2773	GO TO 140
2774	C-----------------------------------------------------------------------
2775	C THE FOLLOWING BLOCK HANDLES PRELIMINARIES NEEDED WHEN JSTART = -1.
2776	C IPUP IS SET TO MITER TO FORCE A MATRIX UPDATE.
2777	C IF AN ORDER INCREASE IS ABOUT TO BE CONSIDERED (IALTH = 1),
2778	C IALTH IS RESET TO 2 TO POSTPONE CONSIDERATION ONE MORE STEP.
2779	C IF THE CALLER HAS CHANGED METH, CFODE IS CALLED TO RESET
2780	C THE COEFFICIENTS OF THE METHOD.
2781	C IF THE CALLER HAS CHANGED MAXORD TO A VALUE LESS THAN THE CURRENT
2782	C ORDER NQ, NQ IS REDUCED TO MAXORD, AND A NEW H CHOSEN ACCORDINGLY.
2783	C IF H IS TO BE CHANGED, YH MUST BE RESCALED.
2784	C IF H OR METH IS BEING CHANGED, IALTH IS RESET TO L = NQ + 1
2785	C TO PREVENT FURTHER CHANGES IN H FOR THAT MANY STEPS.
2786	C-----------------------------------------------------------------------
2787	100 IPUP = MITER
2788	LMAX = MAXORD + 1
2789	IF (IALTH .EQ. 1) IALTH = 2
2790	IF (METH .EQ. MEO) GO TO 110
2791	CALL CFODE (METH, ELCO, TESCO)
2792	MEO = METH
2793	IF (NQ .GT. MAXORD) GO TO 120
2794	IALTH = L
2795	IRET = 1
2796	GO TO 150
2797	110 IF (NQ .LE. MAXORD) GO TO 160
2798	120 NQ = MAXORD
2799	L = LMAX
2800	DO 125 I = 1,L
2801	125 EL(I) = ELCO(I,NQ)
2802	NQNYH = NQ*NYH
2803	RC = RC*EL(1)/EL0
2804	EL0 = EL(1)
2805	CONIT = 0.5D0/REAL(NQ+2)
2806	DDN = VNORM (N, SAVF, EWT)/TESCO(1,L)
2807	EXDN = ONE/REAL(L)
2808	RHDN = ONE/(1.3D0DDN*EXDN + 0.0000013D0)
2809	RH = MIN(RHDN,ONE)
2810	IREDO = 3
2811	IF (H .EQ. HOLD) GO TO 170
2812	RH = MIN(RH,ABS(H/HOLD))
2813	H = HOLD
2814	GO TO 175
2815	C-----------------------------------------------------------------------
2816	C CFODE IS CALLED TO GET ALL THE INTEGRATION COEFFICIENTS FOR THE
2817	C CURRENT METH. THEN THE EL VECTOR AND RELATED CONSTANTS ARE RESET
2818	C WHENEVER THE ORDER NQ IS CHANGED, OR AT THE START OF THE PROBLEM.
2819	C-----------------------------------------------------------------------
2820	140 CALL CFODE (METH, ELCO, TESCO)
2821	150 DO 155 I = 1,L
2822	155 EL(I) = ELCO(I,NQ)
2823	NQNYH = NQ*NYH
2824	RC = RC*EL(1)/EL0
2825	EL0 = EL(1)
2826	CONIT = 0.5D0/REAL(NQ+2)
2827	GO TO (160, 170, 200), IRET
2828	C-----------------------------------------------------------------------
2829	C IF H IS BEING CHANGED, THE H RATIO RH IS CHECKED AGAINST
2830	C RMAX, HMIN, AND HMXI, AND THE YH ARRAY RESCALED. IALTH IS SET TO
2831	C L = NQ + 1 TO PREVENT A CHANGE OF H FOR THAT MANY STEPS, UNLESS
2832	C FORCED BY A CONVERGENCE OR ERROR TEST FAILURE.
2833	C-----------------------------------------------------------------------
2834	160 IF (H .EQ. HOLD) GO TO 200
2835	RH = H/HOLD
2836	H = HOLD
2837	IREDO = 3
2838	GO TO 175
2839	170 RH = MAX(RH,HMIN/ABS(H))
2840	175 RH = MIN(RH,RMAX)
2841	RH = RH/MAX(ONE,ABS(H)HMXIRH)
2842	R = ONE
2843	DO 180 J = 2,L
2844	R = R*RH
2845	DO 180 I = 1,NYH
2846	180 YH(I,J) = YH(I,J)*R
2847	H = H*RH
2848	RC = RC*RH
2849	IALTH = L
2850	IF (IREDO .EQ. 0) GO TO 690
2851	C-----------------------------------------------------------------------
2852	C THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY
2853	C MULTIPLYING THE YH ARRAY BY THE PASCAL TRIANGLE MATRIX.
2854	C RC IS THE RATIO OF NEW TO OLD VALUES OF THE COEFFICIENT H*EL(1).
2855	C WHEN RC DIFFERS FROM 1 BY MORE THAN CCMAX, IPUP IS SET TO MITER
2856	C TO FORCE PJAC TO BE CALLED, IF A JACOBIAN IS INVOLVED.
2857	C IN ANY CASE, PJAC IS CALLED AT LEAST EVERY MSBP STEPS FOR ISOPT = 0,
2858	C AND AT LEAST ONCE EVERY STEP FOR ISOPT = 1.
2859	C-----------------------------------------------------------------------
2860	200 IF (ABS(RC-ONE) .GT. CCMAX) IPUP = MITER
2861	IF (NST .GE. NSLP+MSBP) IPUP = MITER
2862	TN = TN + H
2863	I1 = NQNYH + 1
2864	DO 215 JB = 1,NQ
2865	I1 = I1 - NYH
2866	DO 210 I = I1,NQNYH
2867	210 YH1(I) = YH1(I) + YH1(I+NYH)
2868	215 CONTINUE
2869	C-----------------------------------------------------------------------
2870	C UP TO MAXCOR CORRECTOR ITERATIONS ARE TAKEN. (= 3, FOR ISOPT = 0;
2871	C = 4, FOR ISOPT = 1). A CONVERGENCE TEST IS MADE ON THE R.M.S. NORM
2872	C OF EACH CORRECTION, WEIGHTED BY THE ERROR WEIGHT VECTOR EWT. THE SUM
2873	C OF THE CORRECTIONS IS ACCUMULATED IN THE VECTOR ACOR(I), I = 1,N.
2874	C (ACOR(I), I = N+1,NYH IS LOADED IN SUBROUTINE STESA (ISOPT = 1).)
2875	C THE YH ARRAY IS NOT ALTERED IN THE CORRECTOR LOOP.
2876	C-----------------------------------------------------------------------
2877	220 M = 0
2878	DO 230 I = 1,N
2879	230 Y(I) = YH(I,1)
2880	CALL F (NEQ, TN, Y, PAR, SAVF)
2881	NFE = NFE + 1
2882	IF (IPUP .LE. 0) GO TO 250
2883	C-----------------------------------------------------------------------
2884	C IF INDICATED, THE MATRIX P = I - HEL(1)J IS REEVALUATED AND
2885	C PREPROCESSED BEFORE STARTING THE CORRECTOR ITERATION. IPUP IS SET
2886	C TO 0 AS AN INDICATOR THAT THIS HAS BEEN DONE.
2887	C-----------------------------------------------------------------------
2888	IPUP = 0
2889	RC = ONE
2890	NSLP = NST
2891	CRATE = 0.7D0
2892	CALL PJAC (NEQ, Y, YH, NYH, WM, IWM, EWT, SAVF, ACOR, PAR,
2893	1 F, JAC, JOPT)
2894	IF (IERPJ .NE. 0) GO TO 430
2895	250 DO 260 I = 1,N
2896	260 ACOR(I) = ZERO
2897	270 IF (MITER .NE. 0) GO TO 350
2898	C-----------------------------------------------------------------------
2899	C IN THE CASE OF FUNCTIONAL ITERATION, UPDATE Y DIRECTLY FROM
2900	C THE RESULT OF THE LAST FUNCTION EVALUATION.
2901	C (IF ISOPT = 1, FUNCTIONAL ITERATION IS NOT ALLOWED.)
2902	C-----------------------------------------------------------------------
2903	DO 290 I = 1,N
2904	SAVF(I) = H*SAVF(I) - YH(I,2)
2905	290 Y(I) = SAVF(I) - ACOR(I)
2906	DEL = VNORM (N, Y, EWT)
2907	DO 300 I = 1,N
2908	Y(I) = YH(I,1) + EL(1)*SAVF(I)
2909	300 ACOR(I) = SAVF(I)
2910	GO TO 400
2911	C-----------------------------------------------------------------------
2912	C IN THE CASE OF THE CHORD METHOD, COMPUTE THE CORRECTOR ERROR,
2913	C AND KppSolve THE LINEAR SYSTEM WITH THAT AS RIGHT-HAND SIDE AND
2914	C P AS COEFFICIENT MATRIX.
2915	C-----------------------------------------------------------------------
2916	350 DO 360 I = 1,N
2917	360 Y(I) = H*SAVF(I) - (YH(I,2) + ACOR(I))
2918	CALL SLVS (WM, IWM, Y, SAVF)
2919	IF (IERSL .LT. 0) GO TO 430
2920	IF (IERSL .GT. 0) GO TO 410
2921	DEL = VNORM (N, Y, EWT)
2922	DO 380 I = 1,N
2923	ACOR(I) = ACOR(I) + Y(I)
2924	380 Y(I) = YH(I,1) + EL(1)*ACOR(I)
2925	C-----------------------------------------------------------------------
2926	C TEST FOR CONVERGENCE. IF M.GT.0, AN ESTIMATE OF THE CONVERGENCE
2927	C RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST.
2928	C-----------------------------------------------------------------------
2929	400 IF (M .NE. 0) CRATE = MAX(0.2D0*CRATE,DEL/DELP)
2930	DCON = DELMIN(ONE,1.5D0CRATE)/(TESCO(2,NQ)*CONIT)
2931	IF (DCON .LE. ONE) GO TO 450
2932	M = M + 1
2933	IF (M .EQ. MAXCOR) GO TO 410
2934	IF (M .GE. 2 .AND. DEL .GT. 2.0D0*DELP) GO TO 410
2935	DELP = DEL
2936	CALL F (NEQ, TN, Y, PAR, SAVF)
2937	NFE = NFE + 1
2938	GO TO 270
2939	C-----------------------------------------------------------------------
2940	C THE CORRECTOR ITERATION FAILED TO CONVERGE IN MAXCOR TRIES.
2941	C IF MITER .NE. 0 AND THE JACOBIAN IS OUT OF DATE, PJAC IS CALLED FOR
2942	C THE NEXT TRY. OTHERWISE THE YH ARRAY IS RETRACTED TO ITS VALUES
2943	C BEFORE PREDICTION, AND H IS REDUCED, IF POSSIBLE. IF H CANNOT BE
2944	C REDUCED OR MXNCF FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -2.
2945	C-----------------------------------------------------------------------
2946	410 IF (MITER .EQ. 0 .OR. JCUR .EQ. 1) GO TO 430
2947	ICF = 1
2948	IPUP = MITER
2949	GO TO 220
2950	430 ICF = 2
2951	NCF = NCF + 1
2952	RMAX = 2.0D0
2953	TN = TOLD
2954	I1 = NQNYH + 1
2955	DO 445 JB = 1,NQ
2956	I1 = I1 - NYH
2957	DO 440 I = I1,NQNYH
2958	440 YH1(I) = YH1(I) - YH1(I+NYH)
2959	445 CONTINUE
2960	IF (IERPJ .LT. 0 .OR. IERSL .LT. 0) GO TO 680
2961	IF (ABS(H) .LE. HMIN*1.00001D0) GO TO 670
2962	IF (NCF .EQ. MXNCF) GO TO 670
2963	RH = 0.25D0
2964	IPUP = MITER
2965	IREDO = 1
2966	GO TO 170
2967	C-----------------------------------------------------------------------
2968	C THE CORRECTOR HAS CONVERGED.
2969	C THE LOCAL ERROR TEST IS MADE AND CONTROL PASSES TO STATEMENT 500
2970	C IF IT FAILS. OTHERWISE, STESA IS CALLED (ISOPT = 1) TO PERFORM
2971	C SENSITIVITY CALCULATIONS AT CURRENT STEP SIZE AND ORDER.
2972	C-----------------------------------------------------------------------
2973	450 CONTINUE
2974	IF (M .EQ. 0) DSM = DEL/TESCO(2,NQ)
2975	IF (M .GT. 0) DSM = VNORM (N, ACOR, EWT)/TESCO(2,NQ)
2976	IF (DSM .GT. ONE) GO TO 500
2977	C
2978	IF (ISOPT .EQ. 0) GO TO 460
2979	C-----------------------------------------------------------------------
2980	C CALL STESA TO PERFORM EXPLICIT SENSITIVITY ANALYSIS.
2981	C IF THE LOCAL ERROR TEST FAILS (WITHIN STESA) FOR ANY SOLUTION VECTOR,
2982	C KFLAGS IS SET .LT. 0 AND CONTROL PASSES TO STATEMENT 500 UPON RETURN.
2983	C IN EITHER CASE, JCUR IS SET TO ZERO TO SIGNAL THAT THE JACOBIAN MAY
2984	C NEED UPDATING LATER.
2985	C-----------------------------------------------------------------------
2986	CALL STESA (NEQ, Y, N, NSV, YH, WM, IWM, EWT, SAVF, ACOR,
2987	1 PAR, NRS, F, JAC, DF, PJAC, PDF, SLVS)
2988	IF (IERPJ .NE. 0 .OR. IERSL .NE. 0) GO TO 680
2989	IF (KFLAGS .LT. 0) GO TO 500
2990	C-----------------------------------------------------------------------
2991	C AFTER A SUCCESSFUL STEP, UPDATE THE YH ARRAY.
2992	C CONSIDER CHANGING H IF IALTH = 1. OTHERWISE DECREASE IALTH BY 1.
2993	C IF IALTH IS THEN 1 AND NQ .LT. MAXORD, THEN ACOR IS SAVED FOR
2994	C USE IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP.
2995	C IF A CHANGE IN H IS CONSIDERED, AN INCREASE OR DECREASE IN ORDER
2996	C BY ONE IS CONSIDERED ALSO. A CHANGE IN H IS MADE ONLY IF IT IS BY A
2997	C FACTOR OF AT LEAST 1.1. IF NOT, IALTH IS SET TO 3 TO PREVENT
2998	C TESTING FOR THAT MANY STEPS.
2999	C-----------------------------------------------------------------------
3000	460 JCUR = 0
3001	KFLAG = 0
3002	IREDO = 0
3003	NST = NST + 1
3004	HU = H
3005	NQU = NQ
3006	DO 470 J = 1,L
3007	DO 470 I = 1,NYH
3008	470 YH(I,J) = YH(I,J) + EL(J)*ACOR(I)
3009	IALTH = IALTH - 1
3010	IF (IALTH .EQ. 0) GO TO 520
3011	IF (IALTH .GT. 1) GO TO 700
3012	IF (L .EQ. LMAX) GO TO 700
3013	DO 490 I = 1,NYH
3014	490 YH(I,LMAX) = ACOR(I)
3015	GO TO 700
3016	C-----------------------------------------------------------------------
3017	C THE ERROR TEST FAILED IN EITHER STODE OR STESA.
3018	C KFLAG KEEPS TRACK OF MULTIPLE FAILURES.
3019	C RESTORE TN AND THE YH ARRAY TO THEIR PREVIOUS VALUES, AND PREPARE
3020	C TO TRY THE STEP AGAIN. COMPUTE THE OPTIMUM STEP SIZE FOR THIS OR
3021	C ONE LOWER ORDER. AFTER 2 OR MORE FAILURES, H IS FORCED TO DECREASE
3022	C BY A FACTOR OF 0.2 OR LESS.
3023	C-----------------------------------------------------------------------
3024	500 KFLAG = KFLAG - 1
3025	JCUR = 0
3026	TN = TOLD
3027	I1 = NQNYH + 1
3028	DO 515 JB = 1,NQ
3029	I1 = I1 - NYH
3030	DO 510 I = I1,NQNYH
3031	510 YH1(I) = YH1(I) - YH1(I+NYH)
3032	515 CONTINUE
3033	RMAX = 2.0D0
3034	IF (ABS(H) .LE. HMIN*1.00001D0) GO TO 660
3035	IF (KFLAG .LE. -3) GO TO 640
3036	IREDO = 2
3037	RHUP = ZERO
3038	GO TO 540
3039	C-----------------------------------------------------------------------
3040	*
3041	C REGARDLESS OF THE SUCCESS OR FAILURE OF THE STEP, FACTORS
3042	C RHDN, RHSM, AND RHUP ARE COMPUTED, BY WHICH H COULD BE MULTIPLIED
3043	C AT ORDER NQ - 1, ORDER NQ, OR ORDER NQ + 1, RESPECTIVELY.
3044	C IN THE CASE OF FAILURE, RHUP = 0.0 TO AVOID AN ORDER INCREASE.
3045	C THE LARGEST OF THESE IS DETERMINED AND THE NEW ORDER CHOSEN
3046	C ACCORDINGLY. IF THE ORDER IS TO BE INCREASED, WE COMPUTE ONE
3047	C ADDITIONAL SCALED DERIVATIVE.
3048	C FOR ISOPT = 1, DUPS AND DSMS ARE LOADED WITH THE LARGEST RMS-NORMS
3049	C OBTAINED BY CONSIDERING SEPARATELY THE SENSITIVITY SOLUTION VECTORS.
3050	C-----------------------------------------------------------------------
3051	520 RHUP = ZERO
3052	IF (L .EQ. LMAX) GO TO 540
3053	DO 530 I = 1,N
3054	530 SAVF(I) = ACOR(I) - YH(I,LMAX)
3055	DUP = VNORM (N, SAVF, EWT)/TESCO(3,NQ)
3056	DUP = MAX(DUP,DUPS)
3057	EXUP = ONE/REAL(L+1)
3058	RHUP = ONE/(1.4D0DUP*EXUP + 0.0000014D0)
3059	540 EXSM = ONE/REAL(L)
3060	DSM = MAX(DSM,DSMS)
3061	RHSM = ONE/(1.2D0DSM*EXSM + 0.0000012D0)
3062	RHDN = ZERO
3063	IF (NQ .EQ. 1) GO TO 560
3064	JPOINT = 1
3065	DO 550 J = 1,NSV
3066	DDN = VNORM (N, YH(JPOINT,L), EWT(JPOINT))/TESCO(1,NQ)
3067	DDNS = MAX(DDNS,DDN)
3068	JPOINT = JPOINT + N
3069	550 CONTINUE
3070	DDN = DDNS
3071	DDNS = ZERO
3072	EXDN = ONE/REAL(NQ)
3073	RHDN = ONE/(1.3D0DDN*EXDN + 0.0000013D0)
3074	560 IF (RHSM .GE. RHUP) GO TO 570
3075	IF (RHUP .GT. RHDN) GO TO 590
3076	GO TO 580
3077	570 IF (RHSM .LT. RHDN) GO TO 580
3078	NEWQ = NQ
3079	RH = RHSM
3080	GO TO 620
3081	580 NEWQ = NQ - 1
3082	RH = RHDN
3083	IF (KFLAG .LT. 0 .AND. RH .GT. ONE) RH = ONE
3084	GO TO 620
3085	590 NEWQ = L
3086	RH = RHUP
3087	IF (RH .LT. 1.1D0) GO TO 610
3088	R = EL(L)/REAL(L)
3089	DO 600 I = 1,NYH
3090	600 YH(I,NEWQ+1) = ACOR(I)*R
3091	GO TO 630
3092	610 IALTH = 3
3093	GO TO 700
3094	620 IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.1D0)) GO TO 610
3095	IF (KFLAG .LE. -2) RH = MIN(RH,0.2D0)
3096	C-----------------------------------------------------------------------
3097	C IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS.
3098	C IN ANY CASE H IS RESET ACCORDING TO RH AND THE YH ARRAY IS RESCALED.
3099	C THEN EXIT FROM 690 IF THE STEP WAS OK, OR REDO THE STEP OTHERWISE.
3100	C-----------------------------------------------------------------------
3101	IF (NEWQ .EQ. NQ) GO TO 170
3102	630 NQ = NEWQ
3103	L = NQ + 1
3104	IRET = 2
3105	GO TO 150
3106	C-----------------------------------------------------------------------
3107	C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES HAVE OCCURED.
3108	C IF 10 FAILURES HAVE OCCURRED, EXIT WITH KFLAG = -1.
3109	C IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE
3110	C YH ARRAY HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST
3111	C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 1. THEN
3112	C H IS REDUCED BY A FACTOR OF 10, AND THE STEP IS RETRIED,
3113	C UNTIL IT SUCCEEDS OR H REACHES HMIN.
3114	C-----------------------------------------------------------------------
3115	640 IF (KFLAG .EQ. -10) GO TO 660
3116	RH = 0.1D0
3117	RH = MAX(HMIN/ABS(H),RH)
3118	H = H*RH
3119	DO 645 I = 1,NYH
3120	645 Y(I) = YH(I,1)
3121	CALL F (NEQ, TN, Y, PAR, SAVF)
3122	NFE = NFE + 1
3123	IF (ISOPT .EQ. 0) GO TO 649
3124	CALL SPRIME (NEQ, Y, YH, NYH, N, NSV, WM, IWM, EWT, SAVF, ACOR,
3125	1 ACOR(N+1), PAR, F, JAC, DF, PJAC, PDF)
3126	IF (IERSP .LT. 0) GO TO 680
3127	DO 646 I = N+1,NYH
3128	646 YH(I,2) = H*YH(I,2)
3129	649 DO 650 I = 1,N
3130	650 YH(I,2) = H*SAVF(I)
3131	IPUP = MITER
3132	IALTH = 5
3133	IF (NQ .EQ. 1) GO TO 200
3134	NQ = 1
3135	L = 2
3136	IRET = 3
3137	GO TO 150
3138	C-----------------------------------------------------------------------
3139	C ALL RETURNS ARE MADE THROUGH THIS SECTION. H IS SAVED IN HOLD
3140	C TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP.
3141	C-----------------------------------------------------------------------
3142	660 KFLAG = -1
3143	GO TO 720
3144	670 KFLAG = -2
3145	GO TO 720
3146	680 KFLAG = -3
3147	GO TO 720
3148	690 RMAX = 10.0D0
3149	700 R = ONE/TESCO(2,NQU)
3150	DO 710 I = 1,NYH
3151	710 ACOR(I) = ACOR(I)*R
3152	720 HOLD = H
3153	JSTART = 1
3154	RETURN
3155	C----------------------- END OF SUBROUTINE STODE -----------------------
3156	END
3157	SUBROUTINE CFODE (METH, ELCO, TESCO)
3158	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3159	DIMENSION ELCO(13,12), TESCO(3,12)
3160	C-----------------------------------------------------------------------
3161	C CFODE IS CALLED BY THE INTEGRATOR ROUTINE TO SET COEFFICIENTS
3162	C NEEDED THERE. THE COEFFICIENTS FOR THE CURRENT METHOD, AS
3163	C GIVEN BY THE VALUE OF METH, ARE SET FOR ALL ORDERS AND SAVED.
3164	C THE MAXIMUM ORDER ASSUMED HERE IS 12 IF METH = 1 AND 5 IF METH = 2.
3165	C (A SMALLER VALUE OF THE MAXIMUM ORDER IS ALSO ALLOWED.)
3166	C CFODE IS CALLED ONCE AT THE BEGINNING OF THE PROBLEM,
3167	C AND IS NOT CALLED AGAIN UNLESS AND UNTIL METH IS CHANGED.
3168	C
3169	C THE ELCO ARRAY CONTAINS THE BASIC METHOD COEFFICIENTS.
3170	C THE COEFFICIENTS EL(I), 1 .LE. I .LE. NQ+1, FOR THE METHOD OF
3171	C ORDER NQ ARE STORED IN ELCO(I,NQ). THEY ARE GIVEN BY A GENETRATING
3172	C POLYNOMIAL, I.E.,
3173	C L(X) = EL(1) + EL(2)X + ... + EL(NQ+1)X**NQ.
3174	C FOR THE IMPLICIT ADAMS METHODS, L(X) IS GIVEN BY
3175	C DL/DX = (X+1)(X+2)...*(X+NQ-1)/FACTORIAL(NQ-1), L(-1) = 0.
3176	C FOR THE BDF METHODS, L(X) IS GIVEN BY
3177	C L(X) = (X+1)(X+2) ... *(X+NQ)/K,
3178	C WHERE K = FACTORIAL(NQ)*(1 + 1/2 + ... + 1/NQ).
3179	C
3180	C THE TESCO ARRAY CONTAINS TEST CONSTANTS USED FOR THE
3181	C LOCAL ERROR TEST AND THE SELECTION OF STEP SIZE AND/OR ORDER.
3182	C AT ORDER NQ, TESCO(K,NQ) IS USED FOR THE SELECTION OF STEP
3183	C SIZE AT ORDER NQ - 1 IF K = 1, AT ORDER NQ IF K = 2, AND AT ORDER
3184	C NQ + 1 IF K = 3.
3185	C-----------------------------------------------------------------------
3186	DIMENSION PC(12)
3187	PARAMETER (ONE=1.0D0,ZERO=0.0D0)
3188	C
3189	GO TO (100, 200), METH
3190	C
3191	100 ELCO(1,1) = ONE
3192	ELCO(2,1) = ONE
3193	TESCO(1,1) = ZERO
3194	TESCO(2,1) = 2.0D0
3195	TESCO(1,2) = ONE
3196	TESCO(3,12) = ZERO
3197	PC(1) = ONE
3198	RQFAC = ONE
3199	DO 140 NQ = 2,12
3200	C-----------------------------------------------------------------------
3201	C THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE POLYNOMIAL
3202	C P(X) = (X+1)(X+2)...*(X+NQ-1).
3203	C INITIALLY, P(X) = 1.
3204	C-----------------------------------------------------------------------
3205	RQ1FAC = RQFAC
3206	RQFAC = RQFAC/REAL(NQ)
3207	NQM1 = NQ - 1
3208	FNQM1 = REAL(NQM1)
3209	NQP1 = NQ + 1
3210	C FORM COEFFICIENTS OF P(X)*(X+NQ-1). ----------------------------------
3211	PC(NQ) = ZERO
3212	DO 110 IB = 1,NQM1
3213	I = NQP1 - IB
3214	110 PC(I) = PC(I-1) + FNQM1*PC(I)
3215	PC(1) = FNQM1*PC(1)
3216	C COMPUTE INTEGRAL, -1 TO 0, OF P(X) AND X*P(X). -----------------------
3217	PINT = PC(1)
3218	XPIN = PC(1)/2.0D0
3219	TSIGN = ONE
3220	DO 120 I = 2,NQ
3221	TSIGN = -TSIGN
3222	PINT = PINT + TSIGN*PC(I)/REAL(I)
3223	120 XPIN = XPIN + TSIGN*PC(I)/REAL(I+1)
3224	C STORE COEFFICIENTS IN ELCO AND TESCO. --------------------------------
3225	ELCO(1,NQ) = PINT*RQ1FAC
3226	ELCO(2,NQ) = ONE
3227	DO 130 I = 2,NQ
3228	130 ELCO(I+1,NQ) = RQ1FAC*PC(I)/REAL(I)
3229	AGAMQ = RQFAC*XPIN
3230	RAGQ = ONE/AGAMQ
3231	TESCO(2,NQ) = RAGQ
3232	IF (NQ .LT. 12) TESCO(1,NQP1) = RAGQ*RQFAC/REAL(NQP1)
3233	TESCO(3,NQM1) = RAGQ
3234	140 CONTINUE
3235	RETURN
3236	C
3237	200 PC(1) = ONE
3238	RQ1FAC = ONE
3239	DO 230 NQ = 1,5
3240	C-----------------------------------------------------------------------
3241	C THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE POLYNOMIAL
3242	C P(X) = (X+1)(X+2)...*(X+NQ).
3243	C INITIALLY, P(X) = 1.
3244	C-----------------------------------------------------------------------
3245	FNQ = REAL(NQ)
3246	NQP1 = NQ + 1
3247	C FORM COEFFICIENTS OF P(X)*(X+NQ). ------------------------------------
3248	PC(NQP1) = ZERO
3249	DO 210 IB = 1,NQ
3250	I = NQ + 2 - IB
3251	210 PC(I) = PC(I-1) + FNQ*PC(I)
3252	PC(1) = FNQ*PC(1)
3253	C STORE COEFFICIENTS IN ELCO AND TESCO. --------------------------------
3254	DO 220 I = 1,NQP1
3255	220 ELCO(I,NQ) = PC(I)/PC(2)
3256	ELCO(2,NQ) = ONE
3257	TESCO(1,NQ) = RQ1FAC
3258	TESCO(2,NQ) = REAL(NQP1)/ELCO(1,NQ)
3259	TESCO(3,NQ) = REAL(NQ+2)/ELCO(1,NQ)
3260	RQ1FAC = RQ1FAC/FNQ
3261	230 CONTINUE
3262	RETURN
3263	C----------------------- END OF SUBROUTINE CFODE -----------------------
3264	END
3265	SUBROUTINE SOLSY (WM, IWM, X, TEM)
3266	C IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3267	INCLUDE 'KPP_ROOT_Parameters.h'
3268	INCLUDE 'KPP_ROOT_Sparse.h'
3269	DIMENSION WM(), IWM(), X(), TEM()
3270	PARAMETER (ZERO=0.0D0,ONE=1.0D0)
3271	COMMON /ODE001/ ROWND, ROWNS(173),
3272	2 RDUM1(37), EL0, H, RDUM2(6),
3273	3 IOWND(14), IOWNS(4),
3274	4 IDUM1(4), IERSL, IDUM2(5),
3275	5 MITER, IDUM3(4), N, IDUM4(5)
3276	C-----------------------------------------------------------------------
3277	C THIS ROUTINE MANAGES THE SOLUTION OF THE LINEAR SYSTEM ARISING FROM
3278	C A CHORD ITERATION. IT IS CALLED IF MITER .NE. 0.
3279	C IF MITER IS 1 OR 2, IT CALLS DGESL TO ACCOMPLISH THIS.
3280	C IF MITER = 3 IT UPDATES THE COEFFICIENT H*EL0 IN THE DIAGONAL
3281	C MATRIX, AND THEN COMPUTES THE SOLUTION.
3282	C IF MITER IS 4 OR 5, IT CALLS DGBSL.
3283	C COMMUNICATION WITH SOLSY USES THE FOLLOWING VARIABLES..
3284	C WM = REAL WORK SPACE CONTAINING THE INVERSE DIAGONAL MATRIX IF
3285	C MITER = 3 AND THE LU DECOMPOSITION OF THE MATRIX OTHERWISE.
3286	C STORAGE OF MATRIX ELEMENTS STARTS AT WM(3).
3287	C WM ALSO CONTAINS THE FOLLOWING MATRIX-RELATED DATA..
3288	C WM(1) = SQRT(UROUND) (NOT USED HERE),
3289	C WM(2) = HL0, THE PREVIOUS VALUE OF H*EL0, USED IF MITER = 3.
3290	C IWM = INTEGER WORK SPACE CONTAINING PIVOT INFORMATION, STARTING AT
3291	C IWM(21), IF MITER IS 1, 2, 4, OR 5. IWM ALSO CONTAINS BAND
3292	C PARAMETERS ML = IWM(1) AND MU = IWM(2) IF MITER IS 4 OR 5.
3293	C X = THE RIGHT-HAND SIDE VECTOR ON INPUT, AND THE SOLUTION VECTOR
3294	C ON OUTPUT, OF LENGTH N.
3295	C TEM = VECTOR OF WORK SPACE OF LENGTH N, NOT USED IN THIS VERSION.
3296	C IERSL = OUTPUT FLAG (IN COMMON). IERSL = 0 IF NO TROUBLE OCCURRED.
3297	C IERSL = 1 IF A SINGULAR MATRIX AROSE WITH MITER = 3.
3298	C THIS ROUTINE ALSO USES THE COMMON VARIABLES EL0, H, MITER, AND N.
3299	C-----------------------------------------------------------------------
3300	IERSL = 0
3301	GO TO (100, 100, 300, 400, 400), MITER
3302	C 100 CALL DGESL (WM(3), N, N, IWM(21), X, 0)
3303	100 CALL KppSolve (WM(3), X)
3304	RETURN
3305	C
3306	300 PHL0 = WM(2)
3307	HL0 = H*EL0
3308	WM(2) = HL0
3309	IF (HL0 .EQ. PHL0) GO TO 330
3310	R = HL0/PHL0
3311	DO 320 I = 1,N
3312	DI = ONE - R*(ONE - ONE/WM(I+2))
3313	IF (ABS(DI) .EQ. ZERO) GO TO 390
3314	320 WM(I+2) = ONE/DI
3315	330 DO 340 I = 1,N
3316	340 X(I) = WM(I+2)*X(I)
3317	RETURN
3318	390 IERSL = 1
3319	RETURN
3320	C
3321	400 ML = IWM(1)
3322	MU = IWM(2)
3323	MEBAND = 2*ML + MU + 1
3324	CALL DGBSL (WM(3), MEBAND, N, ML, MU, IWM(21), X, 0)
3325	RETURN
3326	C----------------------- END OF SUBROUTINE SOLSY -----------------------
3327	END
3328	SUBROUTINE EWSET (N, ITOL, RTOL, ATOL, YCUR, EWT)
3329	C-----------------------------------------------------------------------
3330	C THIS SUBROUTINE SETS THE ERROR WEIGHT VECTOR EWT ACCORDING TO
3331	C EWT(I) = RTOL(I)*ABS(YCUR(I)) + ATOL(I), I = 1,...,N,
3332	C WITH THE SUBSCRIPT ON RTOL AND/OR ATOL POSSIBLY REPLACED BY 1 ABOVE,
3333	C DEPENDING ON THE VALUE OF ITOL.
3334	C-----------------------------------------------------------------------
3335	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3336	DIMENSION RTOL(), ATOL(), YCUR(N), EWT(N)
3337	RTOLI = RTOL(1)
3338	ATOLI = ATOL(1)
3339	DO 10 I = 1,N
3340	IF (ITOL .GE. 3) RTOLI = RTOL(I)
3341	IF (ITOL .EQ. 2 .OR. ITOL .EQ. 4) ATOLI = ATOL(I)
3342	EWT(I) = RTOLI*ABS(YCUR(I)) + ATOLI
3343	10 CONTINUE
3344	RETURN
3345	C----------------------- END OF SUBROUTINE EWSET -----------------------
3346	END
3347	DOUBLE PRECISION FUNCTION VNORM (N, V, W)
3348	C-----------------------------------------------------------------------
3349	C THIS FUNCTION ROUTINE COMPUTES THE WEIGHTED ROOT-MEAN-SQUARE NORM
3350	C OF THE VECTOR OF LENGTH N CONTAINED IN THE ARRAY V, WITH WEIGHTS
3351	C CONTAINED IN THE ARRAY W OF LENGTH N..
3352	C VNORM = SQRT( (1/N) * SUM( V(I)W(I) )*2 )
3353	C PROTECTION FOR UNDERFLOW/OVERFLOW IS ACCOMPLISHED USING TWO
3354	C CONSTANTS WHICH ARE HOPEFULLY APPLICABLE FOR ALL MACHINES.
3355	C THESE ARE:
3356	C CUTLO = maximum of SQRT(U/EPS) over all known machines
3357	C CUTHI = minimum of SQRT(Z) over all known machines
3358	C WHERE
3359	C EPS = smallest number s.t. EPS + 1 .GT. 1
3360	C U = smallest positive number (underflow limit)
3361	C Z = largest number (overflow limit)
3362	C
3363	C DETAILS OF THE ALGORITHM AND OF VALUES OF CUTLO AND CUTHI ARE
3364	C FOUND IN THE BLAS ROUTINE SNRM2 (SEE ALSO ALGORITHM 539, TRANS.
3365	C MATH. SOFTWARE, VOL. 5 NO. 3, 1979, 308-323.
3366	C FOR SINGLE PRECISION, THE FOLLOWING VALUES SHOULD BE UNIVERSAL:
3367	C DATA CUTLO,CUTHI /4.441E-16,1.304E19/
3368	C FOR DOUBLE PRECISION, USE
3369	C DATA CUTLO,CUTHI /8.232D-11,1.304D19/
3370	C
3371	C-----------------------------------------------------------------------
3372	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3373	INTEGER NEXT,I,J,N
3374	DIMENSION V(N),W(N)
3375	DATA CUTLO,CUTHI /8.232D-11,1.304D19/
3376	DATA ZERO,ONE/0.0D0,1.0D0/
3377	C BLAS ALGORITHM
3378	NEXT = 1
3379	SUM = ZERO
3380	I = 1
3381	20 SX = V(I)*W(I)
3382	GO TO (30,40,70,80),NEXT
3383	30 IF (ABS(SX).GT.CUTLO) GO TO 110
3384	NEXT = 2
3385	XMAX = ZERO
3386	40 IF (SX.EQ.ZERO) GO TO 130
3387	IF (ABS(SX).GT.CUTLO) GO TO 110
3388	NEXT = 3
3389	GO TO 60
3390	50 I=J
3391	NEXT = 4
3392	SUM = (SUM/SX)/SX
3393	60 XMAX = ABS(SX)
3394	GO TO 90
3395	70 IF(ABS(SX).GT.CUTLO) GO TO 100
3396	80 IF(ABS(SX).LE.XMAX) GO TO 90
3397	SUM = ONE + SUM * (XMAX/SX)**2
3398	XMAX = ABS(SX)
3399	GO TO 130
3400	90 SUM = SUM + (SX/XMAX)**2
3401	GO TO 130
3402	100 SUM = (SUMXMAX)XMAX
3403	110 HITEST = CUTHI/REAL(N)
3404	DO 120 J = I,N
3405	SX = V(J)*W(J)
3406	IF(ABS(SX).GE.HITEST) GO TO 50
3407	SUM = SUM + SX**2
3408	120 CONTINUE
3409	VNORM = SQRT(SUM)
3410	GO TO 140
3411	130 CONTINUE
3412	I = I + 1
3413	IF (I.LE.N) GO TO 20
3414	VNORM = XMAX * SQRT(SUM)
3415	140 CONTINUE
3416	RETURN
3417	C----------------------- END OF FUNCTION VNORM -------------------------
3418	END
3419	SUBROUTINE SVCOM (RSAV, ISAV)
3420	C-----------------------------------------------------------------------
3421	C THIS ROUTINE STORES IN RSAV AND ISAV THE CONTENTS OF COMMON BLOCKS
3422	C ODE001, ODE002 AND EH0001, WHICH ARE USED INTERNALLY IN THE ODESSA
3423	C PACKAGE.
3424	C RSAV = REAL ARRAY OF LENGTH 222 OR MORE.
3425	C ISAV = INTEGER ARRAY OF LENGTH 52 OR MORE.
3426	C-----------------------------------------------------------------------
3427	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3428	DIMENSION RSAV(), ISAV()
3429	COMMON /ODE001/ RODE1(219), IODE1(39)
3430	COMMON /ODE002/ RODE2(3), IODE2(11)
3431	COMMON /EH0001/ IEH(2)
3432	DATA LRODE1/219/, LIODE1/39/, LRODE2/3/, LIODE2/11/
3433	C
3434	DO 10 I = 1,LRODE1
3435	10 RSAV(I) = RODE1(I)
3436	DO 20 I = 1,LRODE2
3437	J = LRODE1 + I
3438	20 RSAV(J) = RODE2(I)
3439	DO 30 I = 1,LIODE1
3440	30 ISAV(I) = IODE1(I)
3441	DO 40 I = 1,LIODE2
3442	J = LIODE1 + I
3443	40 ISAV(J) = IODE2(I)
3444	ISAV(LIODE1+LIODE2+1) = IEH(1)
3445	ISAV(LIODE1+LIODE2+2) = IEH(2)
3446	RETURN
3447	C----------------------- END OF SUBROUTINE SVCOM -----------------------
3448	END
3449	SUBROUTINE RSCOM (RSAV, ISAV)
3450	C-----------------------------------------------------------------------
3451	C THIS ROUTINE RESTORES FROM RSAV AND ISAV THE CONTENTS OF COMMON BLOCKS
3452	C ODE001, ODE002 AND EH0001, WHICH ARE USED INTERNALLY IN THE ODESSSA
3453	C PACKAGE. THIS PRESUMES THAT RSAV AND ISAV WERE LOADED BY MEANS
3454	C OF SUBROUTINE SVCOM OR THE EQUIVALENT.
3455	C-----------------------------------------------------------------------
3456	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3457	DIMENSION RSAV(), ISAV()
3458	COMMON /ODE001/ RODE1(219), IODE1(39)
3459	COMMON /ODE002/ RODE2(3), IODE2(11)
3460	COMMON /EH0001/ IEH(2)
3461	DATA LRODE1/219/, LIODE1/39/, LRODE2/3/, LIODE2/11/
3462	C
3463	DO 10 I = 1,LRODE1
3464	10 RODE1(I) = RSAV(I)
3465	DO 20 I = 1,LRODE2
3466	J = LRODE1 + I
3467	20 RODE2(I) = RSAV(J)
3468	DO 30 I = 1,LIODE1
3469	30 IODE1(I) = ISAV(I)
3470	DO 40 I = 1,LODE2
3471	J = LIODE1 + I
3472	40 IODE2(I) = ISAV(J)
3473	IEH(1) = ISAV(LIODE1+LIODE2+1)
3474	IEH(2) = ISAV(LIODE1+LIODE2+2)
3475	RETURN
3476	C----------------------- END OF SUBROUTINE RSCOM -----------------------
3477	END
3478	SUBROUTINE DGEFA(A,LDA,N,IPVT,INFO)
3479	INTEGER LDA,N,IPVT(*),INFO
3480	DOUBLE PRECISION A(LDA,*)
3481	C
3482	C DGEFA FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION.
3483	C
3484	C DGEFA IS USUALLY CALLED BY DGECO, BUT IT CAN BE CALLED
3485	C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED.
3486	C (TIME FOR DGECO) = (1 + 9/N)*(TIME FOR DGEFA) .
3487	C
3488	C ON ENTRY
3489	C
3490	C A DOUBLE PRECISION(LDA, N)
3491	C THE MATRIX TO BE FACTORED.
3492	C
3493	C LDA INTEGER
3494	C THE LEADING DIMENSION OF THE ARRAY A .
3495	C
3496	C N INTEGER
3497	C THE ORDER OF THE MATRIX A .
3498	C
3499	C ON RETURN
3500	C
3501	C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS
3502	C WHICH WERE USED TO OBTAIN IT.
3503	C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE
3504	C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER
3505	C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR.
3506	C
3507	C IPVT INTEGER(N)
3508	C AN INTEGER VECTOR OF PIVOT INDICES.
3509	C
3510	C INFO INTEGER
3511	C = 0 NORMAL VALUE.
3512	C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR
3513	C CONDITION FOR THIS SUBROUTINE, BUT IT DOES
3514	C INDICATE THAT DGESL OR DGEDI WILL DIVIDE BY ZERO
3515	C IF CALLED. USE RCOND IN DGECO FOR A RELIABLE
3516	C INDICATION OF SINGULARITY.
3517	C
3518	C LINPACK. THIS VERSION DATED 08/14/78 .
3519	C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
3520	C
3521	C SUBROUTINES AND FUNCTIONS
3522	C
3523	C BLAS DAXPY,DSCAL,IDAMAX
3524	C
3525	C INTERNAL VARIABLES
3526	C
3527	DOUBLE PRECISION T
3528	INTEGER IDAMAX,J,K,KP1,L,NM1
3529	C
3530	C
3531	C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
3532	C
3533	INFO = 0
3534	NM1 = N - 1
3535	IF (NM1 .LT. 1) GO TO 70
3536	DO 60 K = 1, NM1
3537	KP1 = K + 1
3538	C
3539	C FIND L = PIVOT INDEX
3540	C
3541	L = IDAMAX(N-K+1,A(K,K),1) + K - 1
3542	IPVT(K) = L
3543	C
3544	C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
3545	C
3546	IF (A(L,K) .EQ. 0.0D0) GO TO 40
3547	C
3548	C INTERCHANGE IF NECESSARY
3549	C
3550	IF (L .EQ. K) GO TO 10
3551	T = A(L,K)
3552	A(L,K) = A(K,K)
3553	A(K,K) = T
3554	10 CONTINUE
3555	C
3556	C COMPUTE MULTIPLIERS
3557	C
3558	T = -1.0D0/A(K,K)
3559	CALL DSCAL(N-K,T,A(K+1,K),1)
3560	C
3561	C ROW ELIMINATION WITH COLUMN INDEXING
3562	C
3563	DO 30 J = KP1, N
3564	T = A(L,J)
3565	IF (L .EQ. K) GO TO 20
3566	A(L,J) = A(K,J)
3567	A(K,J) = T
3568	20 CONTINUE
3569	CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1)
3570	30 CONTINUE
3571	GO TO 50
3572	40 CONTINUE
3573	INFO = K
3574	50 CONTINUE
3575	60 CONTINUE
3576	70 CONTINUE
3577	IPVT(N) = N
3578	IF (A(N,N) .EQ. 0.0D0) INFO = N
3579	RETURN
3580	END
3581	SUBROUTINE DGESL(A,LDA,N,IPVT,B,JOB)
3582	INTEGER LDA,N,IPVT(*),JOB
3583	DOUBLE PRECISION A(LDA,),B()
3584	C
3585	C DGESL KppSolveS THE DOUBLE PRECISION SYSTEM
3586	C A * X = B OR TRANS(A) * X = B
3587	C USING THE FACTORS COMPUTED BY DGECO OR DGEFA.
3588	C
3589	C ON ENTRY
3590	C
3591	C A DOUBLE PRECISION(LDA, N)
3592	C THE OUTPUT FROM DGECO OR DGEFA.
3593	C
3594	C LDA INTEGER
3595	C THE LEADING DIMENSION OF THE ARRAY A .
3596	C
3597	C N INTEGER
3598	C THE ORDER OF THE MATRIX A .
3599	C
3600	C IPVT INTEGER(N)
3601	C THE PIVOT VECTOR FROM DGECO OR DGEFA.
3602	C
3603	C B DOUBLE PRECISION(N)
3604	C THE RIGHT HAND SIDE VECTOR.
3605	C
3606	C JOB INTEGER
3607	C = 0 TO KppSolve A*X = B ,
3608	C = NONZERO TO KppSolve TRANS(A)*X = B WHERE
3609	C TRANS(A) IS THE TRANSPOSE.
3610	C
3611	C ON RETURN
3612	C
3613	C B THE SOLUTION VECTOR X .
3614	C
3615	C ERROR CONDITION
3616	C
3617	C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A
3618	C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY
3619	C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER
3620	C SETTING OF LDA . IT WILL NOT OCCUR IF THE SUBROUTINES ARE
3621	C CALLED CORRECTLY AND IF DGECO HAS SET RCOND .GT. 0.0
3622	C OR DGEFA HAS SET INFO .EQ. 0 .
3623	C
3624	C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX
3625	C WITH P COLUMNS
3626	C CALL DGECO(A,LDA,N,IPVT,RCOND,Z)
3627	C IF (RCOND IS TOO SMALL) GO TO ...
3628	C DO 10 J = 1, P
3629	C CALL DGESL(A,LDA,N,IPVT,C(1,J),0)
3630	C 10 CONTINUE
3631	C
3632	C LINPACK. THIS VERSION DATED 08/14/78 .
3633	C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
3634	C
3635	C SUBROUTINES AND FUNCTIONS
3636	C
3637	C BLAS DAXPY,DDOT
3638	C
3639	C INTERNAL VARIABLES
3640	C
3641	DOUBLE PRECISION DDOT,T
3642	INTEGER K,KB,L,NM1
3643	C
3644	NM1 = N - 1
3645	IF (JOB .NE. 0) GO TO 50
3646	C
3647	C JOB = 0 , KppSolve A * X = B
3648	C FIRST KppSolve L*Y = B
3649	C
3650	IF (NM1 .LT. 1) GO TO 30
3651	DO 20 K = 1, NM1
3652	L = IPVT(K)
3653	T = B(L)
3654	IF (L .EQ. K) GO TO 10
3655	B(L) = B(K)
3656	B(K) = T
3657	10 CONTINUE
3658	CALL DAXPY(N-K,T,A(K+1,K),1,B(K+1),1)
3659	20 CONTINUE
3660	30 CONTINUE
3661	C
3662	C NOW KppSolve U*X = Y
3663	C
3664	DO 40 KB = 1, N
3665	K = N + 1 - KB
3666	B(K) = B(K)/A(K,K)
3667	T = -B(K)
3668	CALL DAXPY(K-1,T,A(1,K),1,B(1),1)
3669	40 CONTINUE
3670	GO TO 100
3671	50 CONTINUE
3672	C
3673	C JOB = NONZERO, KppSolve TRANS(A) * X = B
3674	C FIRST KppSolve TRANS(U)*Y = B
3675	C
3676	DO 60 K = 1, N
3677	T = DDOT(K-1,A(1,K),1,B(1),1)
3678	B(K) = (B(K) - T)/A(K,K)
3679	60 CONTINUE
3680	C
3681	C NOW KppSolve TRANS(L)*X = Y
3682	C
3683	IF (NM1 .LT. 1) GO TO 90
3684	DO 80 KB = 1, NM1
3685	K = N - KB
3686	B(K) = B(K) + DDOT(N-K,A(K+1,K),1,B(K+1),1)
3687	L = IPVT(K)
3688	IF (L .EQ. K) GO TO 70
3689	T = B(L)
3690	B(L) = B(K)
3691	B(K) = T
3692	70 CONTINUE
3693	80 CONTINUE
3694	90 CONTINUE
3695	100 CONTINUE
3696	RETURN
3697	END
3698	SUBROUTINE DGBFA(ABD,LDA,N,ML,MU,IPVT,INFO)
3699	INTEGER LDA,N,ML,MU,IPVT(*),INFO
3700	DOUBLE PRECISION ABD(LDA,*)
3701	C
3702	C DGBFA FACTORS A DOUBLE PRECISION BAND MATRIX BY ELIMINATION.
3703	C
3704	C DGBFA IS USUALLY CALLED BY DGBCO, BUT IT CAN BE CALLED
3705	C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED.
3706	C
3707	C ON ENTRY
3708	C
3709	C ABD DOUBLE PRECISION(LDA, N)
3710	C CONTAINS THE MATRIX IN BAND STORAGE. THE COLUMNS
3711	C OF THE MATRIX ARE STORED IN THE COLUMNS OF ABD AND
3712	C THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS
3713	C ML+1 THROUGH 2*ML+MU+1 OF ABD .
3714	C SEE THE COMMENTS BELOW FOR DETAILS.
3715	C
3716	C LDA INTEGER
3717	C THE LEADING DIMENSION OF THE ARRAY ABD .
3718	C LDA MUST BE .GE. 2*ML + MU + 1 .
3719	C
3720	C N INTEGER
3721	C THE ORDER OF THE ORIGINAL MATRIX.
3722	C
3723	C ML INTEGER
3724	C NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL.
3725	C 0 .LE. ML .LT. N .
3726	C
3727	C MU INTEGER
3728	C NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL.
3729	C 0 .LE. MU .LT. N .
3730	C MORE EFFICIENT IF ML .LE. MU .
3731	C ON RETURN
3732	C
3733	C ABD AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND
3734	C THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT.
3735	C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE
3736	C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER
3737	C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR.
3738	C
3739	C IPVT INTEGER(N)
3740	C AN INTEGER VECTOR OF PIVOT INDICES.
3741	C
3742	C INFO INTEGER
3743	C = 0 NORMAL VALUE.
3744	C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR
3745	C CONDITION FOR THIS SUBROUTINE, BUT IT DOES
3746	C INDICATE THAT DGBSL WILL DIVIDE BY ZERO IF
3747	C CALLED. USE RCOND IN DGBCO FOR A RELIABLE
3748	C INDICATION OF SINGULARITY.
3749	C
3750	C BAND STORAGE
3751	C
3752	C IF A IS A BAND MATRIX, THE FOLLOWING PROGRAM SEGMENT
3753	C WILL SET UP THE INPUT.
3754	C
3755	C ML = (BAND WIDTH BELOW THE DIAGONAL)
3756	C MU = (BAND WIDTH ABOVE THE DIAGONAL)
3757	C M = ML + MU + 1
3758	C DO 20 J = 1, N
3759	C I1 = MAX0(1, J-MU)
3760	C I2 = MIN0(N, J+ML)
3761	C DO 10 I = I1, I2
3762	C K = I - J + M
3763	C ABD(K,J) = A(I,J)
3764	C 10 CONTINUE
3765	C 20 CONTINUE
3766	C
3767	C THIS USES ROWS ML+1 THROUGH 2*ML+MU+1 OF ABD .
3768	C IN ADDITION, THE FIRST ML ROWS IN ABD ARE USED FOR
3769	C ELEMENTS GENERATED DURING THE TRIANGULARIZATION.
3770	C THE TOTAL NUMBER OF ROWS NEEDED IN ABD IS 2*ML+MU+1 .
3771	C THE ML+MU BY ML+MU UPPER LEFT TRIANGLE AND THE
3772	C ML BY ML LOWER RIGHT TRIANGLE ARE NOT REFERENCED.
3773	C
3774	C LINPACK. THIS VERSION DATED 08/14/78 .
3775	C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
3776	C
3777	C SUBROUTINES AND FUNCTIONS
3778	C
3779	C BLAS DAXPY,DSCAL,IDAMAX
3780	C FORTRAN MAX0,MIN0
3781	C
3782	C INTERNAL VARIABLES
3783	C
3784	DOUBLE PRECISION T
3785	INTEGER I,IDAMAX,I0,J,JU,JZ,J0,J1,K,KP1,L,LM,M,MM,NM1
3786	C
3787	C
3788	M = ML + MU + 1
3789	INFO = 0
3790	C
3791	C ZERO INITIAL FILL-IN COLUMNS
3792	C
3793	J0 = MU + 2
3794	J1 = MIN0(N,M) - 1
3795	IF (J1 .LT. J0) GO TO 30
3796	DO 20 JZ = J0, J1
3797	I0 = M + 1 - JZ
3798	DO 10 I = I0, ML
3799	ABD(I,JZ) = 0.0D0
3800	10 CONTINUE
3801	20 CONTINUE
3802	30 CONTINUE
3803	JZ = J1
3804	JU = 0
3805	C
3806	C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
3807	C
3808	NM1 = N - 1
3809	IF (NM1 .LT. 1) GO TO 130
3810	DO 120 K = 1, NM1
3811	KP1 = K + 1
3812	C
3813	C ZERO NEXT FILL-IN COLUMN
3814	C
3815	JZ = JZ + 1
3816	IF (JZ .GT. N) GO TO 50
3817	IF (ML .LT. 1) GO TO 50
3818	DO 40 I = 1, ML
3819	ABD(I,JZ) = 0.0D0
3820	40 CONTINUE
3821	50 CONTINUE
3822	C
3823	C FIND L = PIVOT INDEX
3824	C
3825	LM = MIN0(ML,N-K)
3826	L = IDAMAX(LM+1,ABD(M,K),1) + M - 1
3827	IPVT(K) = L + K - M
3828	C
3829	C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
3830	C
3831	IF (ABD(L,K) .EQ. 0.0D0) GO TO 100
3832	C
3833	C INTERCHANGE IF NECESSARY
3834	C
3835	IF (L .EQ. M) GO TO 60
3836	T = ABD(L,K)
3837	ABD(L,K) = ABD(M,K)
3838	ABD(M,K) = T
3839	60 CONTINUE
3840	C
3841	C COMPUTE MULTIPLIERS
3842	C
3843	T = -1.0D0/ABD(M,K)
3844	CALL DSCAL(LM,T,ABD(M+1,K),1)
3845	C
3846	C ROW ELIMINATION WITH COLUMN INDEXING
3847	C
3848	JU = MIN0(MAX0(JU,MU+IPVT(K)),N)
3849	MM = M
3850	IF (JU .LT. KP1) GO TO 90
3851	DO 80 J = KP1, JU
3852	L = L - 1
3853	MM = MM - 1
3854	T = ABD(L,J)
3855	IF (L .EQ. MM) GO TO 70
3856	ABD(L,J) = ABD(MM,J)
3857	ABD(MM,J) = T
3858	70 CONTINUE
3859	CALL DAXPY(LM,T,ABD(M+1,K),1,ABD(MM+1,J),1)
3860	80 CONTINUE
3861	90 CONTINUE
3862	GO TO 110
3863	100 CONTINUE
3864	INFO = K
3865	110 CONTINUE
3866	120 CONTINUE
3867	130 CONTINUE
3868	IPVT(N) = N
3869	IF (ABD(M,N) .EQ. 0.0D0) INFO = N
3870	RETURN
3871	END
3872	SUBROUTINE DGBSL(ABD,LDA,N,ML,MU,IPVT,B,JOB)
3873	INTEGER LDA,N,ML,MU,IPVT(*),JOB
3874	DOUBLE PRECISION ABD(LDA,),B()
3875	C
3876	C DGBSL KppSolveS THE DOUBLE PRECISION BAND SYSTEM
3877	C A * X = B OR TRANS(A) * X = B
3878	C USING THE FACTORS COMPUTED BY DGBCO OR DGBFA.
3879	C
3880	C ON ENTRY
3881	C
3882	C ABD DOUBLE PRECISION(LDA, N)
3883	C THE OUTPUT FROM DGBCO OR DGBFA.
3884	C
3885	C LDA INTEGER
3886	C THE LEADING DIMENSION OF THE ARRAY ABD .
3887	C
3888	C N INTEGER
3889	C THE ORDER OF THE ORIGINAL MATRIX.
3890	C
3891	C ML INTEGER
3892	C NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL.
3893	C
3894	C MU INTEGER
3895	C NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL.
3896	C
3897	C IPVT INTEGER(N)
3898	C THE PIVOT VECTOR FROM DGBCO OR DGBFA.
3899	C
3900	C B DOUBLE PRECISION(N)
3901	C THE RIGHT HAND SIDE VECTOR.
3902	C
3903	C JOB INTEGER
3904	C = 0 TO KppSolve A*X = B ,
3905	C = NONZERO TO KppSolve TRANS(A)*X = B , WHERE
3906	C TRANS(A) IS THE TRANSPOSE.
3907	C
3908	C ON RETURN
3909	C
3910	C B THE SOLUTION VECTOR X .
3911	C
3912	C ERROR CONDITION
3913	C
3914	C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A
3915	C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY
3916	C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER
3917	C SETTING OF LDA . IT WILL NOT OCCUR IF THE SUBROUTINES ARE
3918	C CALLED CORRECTLY AND IF DGBCO HAS SET RCOND .GT. 0.0
3919	C OR DGBFA HAS SET INFO .EQ. 0 .
3920	C
3921	C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX
3922	C WITH P COLUMNS
3923	C CALL DGBCO(ABD,LDA,N,ML,MU,IPVT,RCOND,Z)
3924	C IF (RCOND IS TOO SMALL) GO TO ...
3925	C DO 10 J = 1, P
3926	C CALL DGBSL(ABD,LDA,N,ML,MU,IPVT,C(1,J),0)
3927	C 10 CONTINUE
3928	C
3929	C LINPACK. THIS VERSION DATED 08/14/78 .
3930	C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
3931	C
3932	C SUBROUTINES AND FUNCTIONS
3933	C
3934	C BLAS DAXPY,DDOT
3935	C FORTRAN MIN0
3936	C
3937	C INTERNAL VARIABLES
3938	C
3939	DOUBLE PRECISION DDOT,T
3940	INTEGER K,KB,L,LA,LB,LM,M,NM1
3941	C
3942	M = MU + ML + 1
3943	NM1 = N - 1
3944	IF (JOB .NE. 0) GO TO 50
3945	C
3946	C JOB = 0 , KppSolve A * X = B
3947	C FIRST KppSolve L*Y = B
3948	C
3949	IF (ML .EQ. 0) GO TO 30
3950	IF (NM1 .LT. 1) GO TO 30
3951	DO 20 K = 1, NM1
3952	LM = MIN0(ML,N-K)
3953	L = IPVT(K)
3954	T = B(L)
3955	IF (L .EQ. K) GO TO 10
3956	B(L) = B(K)
3957	B(K) = T
3958	10 CONTINUE
3959	CALL DAXPY(LM,T,ABD(M+1,K),1,B(K+1),1)
3960	20 CONTINUE
3961	30 CONTINUE
3962	C
3963	C NOW KppSolve U*X = Y
3964	C
3965	DO 40 KB = 1, N
3966	K = N + 1 - KB
3967	B(K) = B(K)/ABD(M,K)
3968	LM = MIN0(K,M) - 1
3969	LA = M - LM
3970	LB = K - LM
3971	T = -B(K)
3972	CALL DAXPY(LM,T,ABD(LA,K),1,B(LB),1)
3973	40 CONTINUE
3974	GO TO 100
3975	50 CONTINUE
3976	C
3977	C JOB = NONZERO, KppSolve TRANS(A) * X = B
3978	C FIRST KppSolve TRANS(U)*Y = B
3979	C
3980	DO 60 K = 1, N
3981	LM = MIN0(K,M) - 1
3982	LA = M - LM
3983	LB = K - LM
3984	T = DDOT(LM,ABD(LA,K),1,B(LB),1)
3985	B(K) = (B(K) - T)/ABD(M,K)
3986	60 CONTINUE
3987	C
3988	C NOW KppSolve TRANS(L)*X = Y
3989	C
3990	IF (ML .EQ. 0) GO TO 90
3991	IF (NM1 .LT. 1) GO TO 90
3992	DO 80 KB = 1, NM1
3993	K = N - KB
3994	LM = MIN0(ML,N-K)
3995	B(K) = B(K) + DDOT(LM,ABD(M+1,K),1,B(K+1),1)
3996	L = IPVT(K)
3997	IF (L .EQ. K) GO TO 70
3998	T = B(L)
3999	B(L) = B(K)
4000	B(K) = T
4001	70 CONTINUE
4002	80 CONTINUE
4003	90 CONTINUE
4004	100 CONTINUE
4005	RETURN
4006	END
4007	SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)
4008	C
4009	C CONSTANT TIMES A VECTOR PLUS A VECTOR.
4010	C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
4011	C JACK DONGARRA, LINPACK, 3/11/78.
4012	C
4013	DOUBLE PRECISION DX(),DY(),DA
4014	INTEGER I,INCX,INCY,IX,IY,M,MP1,N
4015	C
4016	IF(N.LE.0)RETURN
4017	IF (DA .EQ. 0.0D0) RETURN
4018	IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
4019	C
4020	C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
4021	C NOT EQUAL TO 1
4022	C
4023	IX = 1
4024	IY = 1
4025	IF(INCX.LT.0)IX = (-N+1)*INCX + 1
4026	IF(INCY.LT.0)IY = (-N+1)*INCY + 1
4027	DO 10 I = 1,N
4028	DY(IY) = DY(IY) + DA*DX(IX)
4029	IX = IX + INCX
4030	IY = IY + INCY
4031	10 CONTINUE
4032	RETURN
4033	C
4034	C CODE FOR BOTH INCREMENTS EQUAL TO 1
4035	C
4036	C
4037	C CLEAN-UP LOOP
4038	C
4039	20 M = MOD(N,4)
4040	IF( M .EQ. 0 ) GO TO 40
4041	DO 30 I = 1,M
4042	DY(I) = DY(I) + DA*DX(I)
4043	30 CONTINUE
4044	IF( N .LT. 4 ) RETURN
4045	40 MP1 = M + 1
4046	DO 50 I = MP1,N,4
4047	DY(I) = DY(I) + DA*DX(I)
4048	DY(I + 1) = DY(I + 1) + DA*DX(I + 1)
4049	DY(I + 2) = DY(I + 2) + DA*DX(I + 2)
4050	DY(I + 3) = DY(I + 3) + DA*DX(I + 3)
4051	50 CONTINUE
4052	RETURN
4053	END
4054	SUBROUTINE DSCAL(N,DA,DX,INCX)
4055	C
4056	C SCALES A VECTOR BY A CONSTANT.
4057	C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE.
4058	C JACK DONGARRA, LINPACK, 3/11/78.
4059	C
4060	DOUBLE PRECISION DA,DX(*)
4061	INTEGER I,INCX,M,MP1,N,NINCX
4062	C
4063	IF(N.LE.0)RETURN
4064	IF(INCX.EQ.1)GO TO 20
4065	C
4066	C CODE FOR INCREMENT NOT EQUAL TO 1
4067	*
4068	C
4069	NINCX = N*INCX
4070	DO 10 I = 1,NINCX,INCX
4071	DX(I) = DA*DX(I)
4072	10 CONTINUE
4073	RETURN
4074	C
4075	C CODE FOR INCREMENT EQUAL TO 1
4076	C
4077	C
4078	C CLEAN-UP LOOP
4079	C
4080	20 M = MOD(N,5)
4081	IF( M .EQ. 0 ) GO TO 40
4082	DO 30 I = 1,M
4083	DX(I) = DA*DX(I)
4084	30 CONTINUE
4085	IF( N .LT. 5 ) RETURN
4086	40 MP1 = M + 1
4087	DO 50 I = MP1,N,5
4088	DX(I) = DA*DX(I)
4089	DX(I + 1) = DA*DX(I + 1)
4090	DX(I + 2) = DA*DX(I + 2)
4091	DX(I + 3) = DA*DX(I + 3)
4092	DX(I + 4) = DA*DX(I + 4)
4093	50 CONTINUE
4094	RETURN
4095	END
4096	DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
4097	C
4098	C FORMS THE DOT PRODUCT OF TWO VECTORS.
4099	C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
4100	C JACK DONGARRA, LINPACK, 3/11/78.
4101	C
4102	DOUBLE PRECISION DX(),DY(),DTEMP
4103	INTEGER I,INCX,INCY,IX,IY,M,MP1,N
4104	C
4105	DDOT = 0.0D0
4106	DTEMP = 0.0D0
4107	IF(N.LE.0)RETURN
4108	IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
4109	C
4110	C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
4111	C NOT EQUAL TO 1
4112	C
4113	IX = 1
4114	IY = 1
4115	IF(INCX.LT.0)IX = (-N+1)*INCX + 1
4116	IF(INCY.LT.0)IY = (-N+1)*INCY + 1
4117	DO 10 I = 1,N
4118	DTEMP = DTEMP + DX(IX)*DY(IY)
4119	IX = IX + INCX
4120	IY = IY + INCY
4121	10 CONTINUE
4122	DDOT = DTEMP
4123	RETURN
4124	C
4125	C CODE FOR BOTH INCREMENTS EQUAL TO 1
4126	C
4127	C
4128	C CLEAN-UP LOOP
4129	C
4130	20 M = MOD(N,5)
4131	IF( M .EQ. 0 ) GO TO 40
4132	DO 30 I = 1,M
4133	DTEMP = DTEMP + DX(I)*DY(I)
4134	30 CONTINUE
4135	IF( N .LT. 5 ) GO TO 60
4136	40 MP1 = M + 1
4137	DO 50 I = MP1,N,5
4138	DTEMP = DTEMP + DX(I)DY(I) + DX(I + 1)DY(I + 1) +
4139	* DX(I + 2)DY(I + 2) + DX(I + 3)DY(I + 3) + DX(I + 4)*DY(I + 4)
4140	50 CONTINUE
4141	60 DDOT = DTEMP
4142	RETURN
4143	END
4144	INTEGER FUNCTION IDAMAX(N,DX,INCX)
4145	C
4146	C FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE.
4147	C JACK DONGARRA, LINPACK, 3/11/78.
4148	C
4149	DOUBLE PRECISION DX(*),DMAX
4150	INTEGER I,INCX,IX,N
4151	C
4152	IDAMAX = 0
4153	IF( N .LT. 1 ) RETURN
4154	IDAMAX = 1
4155	IF(N.EQ.1)RETURN
4156	IF(INCX.EQ.1)GO TO 20
4157	C
4158	C CODE FOR INCREMENT NOT EQUAL TO 1
4159	C
4160	IX = 1
4161	DMAX = DABS(DX(1))
4162	IX = IX + INCX
4163	DO 10 I = 2,N
4164	IF(DABS(DX(IX)).LE.DMAX) GO TO 5
4165	IDAMAX = I
4166	DMAX = DABS(DX(IX))
4167	5 IX = IX + INCX
4168	10 CONTINUE
4169	RETURN
4170	C
4171	C CODE FOR INCREMENT EQUAL TO 1
4172	C
4173	20 DMAX = DABS(DX(1))
4174	DO 30 I = 2,N
4175	IF(DABS(DX(I)).LE.DMAX) GO TO 30
4176	IDAMAX = I
4177	DMAX = DABS(DX(I))
4178	30 CONTINUE
4179	RETURN
4180	END
4181	DOUBLE PRECISION FUNCTION D1MACH (IDUM)
4182	INTEGER IDUM
4183	C-----------------------------------------------------------------------
4184	C THIS ROUTINE COMPUTES THE UNIT ROUNDOFF OF THE MACHINE IN DOUBLE
4185	C PRECISION. THIS IS DEFINED AS THE SMALLEST POSITIVE MACHINE NUMBER
4186	C U SUCH THAT 1.0D0 + U .NE. 1.0D0 (IN DOUBLE PRECISION).
4187	C-----------------------------------------------------------------------
4188	DOUBLE PRECISION U, COMP
4189	U = 1.0D0
4190	10 U = U*0.5D0
4191	COMP = 1.0D0 + U
4192	IF (COMP .NE. 1.0D0) GO TO 10
4193	D1MACH = U*2.0D0
4194	RETURN
4195	C----------------------- END OF FUNCTION D1MACH ------------------------
4196	END
4197	SUBROUTINE XERR (MSG, NERR, IERT, NI, I1, I2, NR, R1, R2)
4198	INTEGER NERR, IERT, NI, I1, I2, NR,
4199	1 LUN, LUNIT, MESFLG
4200	DOUBLE PRECISION R1, R2
4201	CHARACTER() MSG
4202	C-------------------------------------------------------------------
4203	C
4204	C ALL ARGUMENTS ARE INPUT ARGUMENTS.
4205	C
4206	C MSG = THE MESSAGE (CHARACTER VARIABLE)
4207	C NERR = THE ERROR NUMBER (NOT USED).
4208	C IERT = THE ERROR TYPE..
4209	C 1 MEANS RECOVERABLE (CONTROL RETURNS TO CALLER).
4210	C 2 MEANS FATAL (RUN IS ABORTED--SEE NOTE BELOW).
4211	C NI = NUMBER OF INTEGERS (0, 1, OR 2) TO BE PRINTED WITH MESSAGE.
4212	C I1,I2 = INTEGERS TO BE PRINTED, DEPENDING ON NI.
4213	C NR = NUMBER OF REALS (0, 1, OR 2) TO BE PRINTED WITH MESSAGE.
4214	C R1,R2 = REALS TO BE PRINTED, DEPENDING ON NR.
4215	C
4216	C NOTES:
4217	C 1. THE DIMENSION OF MSG IS ASSUMED TO BE AT MOST 60.
4218	C (MULTI-LINE MESSAGES ARE GENERATED BY REPEATED CALLS.)
4219	C 2. IF IERT = 2, CONTROL PASSES TO THE STATEMENT STOP
4220	C TO ABORT THE RUN. THIS STATEMENT MAY BE MACHINE-DEPENDENT.
4221	C 3. R1 AND R2 ARE ASSUMED TO BE IN DOUBLE PRECISION AND ARE PRINTED
4222	C IN D21.13 FORMAT.
4223	C 4. THE COMMON BLOCK /EH0001/ BELOW IS DATA-LOADED (A MACHINE-
4224	C DEPENDENT FEATURE) WITH DEFAULT VALUES.
4225	C THIS BLOCK IS NEEDED FOR PROPER RETENTION OF PARAMETERS USED BY
4226	C THIS ROUTINE WHICH THE USER CAN RESET BY CALLING XSETF OR XSETUN.
4227	C THE VARIABLES IN THIS BLOCK ARE AS FOLLOWS..
4228	C MESFLG = PRINT CONTROL FLAG..
4229	C 1 MEANS PRINT ALL MESSAGES (THE DEFAULT).
4230	C 0 MEANS NO PRINTING.
4231	C LUNIT = LOGICAL UNIT NUMBER FOR MESSAGES.
4232	C THE DEFAULT IS 6 (MACHINE-DEPENDENT).
4233	C 5. TO CHANGE THE DEFAULT OUTPUT UNIT, CHANGE THE DATA STATEMENT
4234	C IN THE BLOCK DATA SUBPROGRAM BELOW.
4235	C
4236	C FOR A DIFFERENT RUN-ABORT COMMAND, CHANGE THE STATEMENT FOLLOWING
4237	C STATEMENT 100 AT THE END.
4238	C-----------------------------------------------------------------------
4239	COMMON /EH0001/ MESFLG, LUNIT
4240	IF (MESFLG .EQ. 0) GO TO 100
4241	C GET LOGICAL UNIT NUMBER. ---------------------------------------------
4242	LUN = LUNIT
4243	C WRITE THE MESSAGE. ---------------------------------------------------
4244	WRITE (LUN, 10) MSG
4245	10 FORMAT(1X,A)
4246	C-----------------------------------------------------------------------
4247	IF (NI .EQ. 1) WRITE (LUN, 20) I1
4248	20 FORMAT(6X,'IN ABOVE MESSAGE, I1 = ',I10)
4249	IF (NI .EQ. 2) WRITE (LUN, 30) I1,I2
4250	30 FORMAT(6X,'IN ABOVE MESSAGE, I1 = ',I10,3X,'I2 = ',I10)
4251	IF (NR .EQ. 1) WRITE (LUN, 40) R1
4252	40 FORMAT(6X,'IN ABOVE MESSAGE, R1 = ',D21.13)
4253	IF (NR .EQ. 2) WRITE (LUN, 50) R1,R2
4254	50 FORMAT(6X,'IN ABOVE, R1 = ',D21.13,3X,'R2 = ',D21.13)
4255	C ABORT THE RUN IF IERT = 2. -------------------------------------------
4256	100 IF (IERT .NE. 2) RETURN
4257	STOP
4258	C----------------------- END OF SUBROUTINE XERR ----------------------
4259	END
4260	SUBROUTINE XSETF (MFLAG)
4261	C
4262	C THIS ROUTINE RESETS THE PRINT CONTROL FLAG MFLAG.
4263	C
4264	INTEGER MFLAG, MESFLG, LUNIT
4265	COMMON /EH0001/ MESFLG, LUNIT
4266	C
4267	IF (MFLAG .EQ. 0 .OR. MFLAG .EQ. 1) MESFLG = MFLAG
4268	RETURN
4269	C----------------------- END OF SUBROUTINE XSETF -----------------------
4270	END
4271	SUBROUTINE XSETUN (LUN)
4272	C
4273	C THIS ROUTINE RESETS THE LOGICAL UNIT NUMBER FOR MESSAGES.
4274	C
4275	INTEGER LUN, MESFLG, LUNIT
4276	COMMON /EH0001/ MESFLG, LUNIT
4277	C
4278	IF (LUN .GT. 0) LUNIT = LUN
4279	RETURN
4280	C----------------------- END OF SUBROUTINE XSETUN ----------------------
4281	END
4282	BLOCK DATA
4283	C-----------------------------------------------------------------------
4284	C THIS DATA SUBPROGRAM LOADS VARIABLES INTO THE INTERNAL COMMON
4285	C BLOCKS USED BY ODESSA AND ITS VARIANTS. THE VARIABLES ARE
4286	C DEFINED AS FOLLOWS..
4287	C ILLIN = COUNTER FOR THE NUMBER OF CONSECUTIVE TIMES THE PACKAGE
4288	C WAS CALLED WITH ILLEGAL INPUT. THE RUN IS STOPPED WHEN
4289	C ILLIN REACHES 5.
4290	C NTREP = COUNTER FOR THE NUMBER OF CONSECUTIVE TIMES THE PACKAGE
4291	C WAS CALLED WITH ISTATE = 1 AND TOUT = T. THE RUN IS
4292	C STOPPED WHEN NTREP REACHES 5.
4293	C MESFLG = FLAG TO CONTROL PRINTING OF ERROR MESSAGES. 1 MEANS PRINT,
4294	C 0 MEANS NO PRINTING.
4295	C LUNIT = DEFAULT VALUE OF LOGICAL UNIT NUMBER FOR PRINTING OF ERROR
4296	C MESSAGES.
4297	C-----------------------------------------------------------------------
4298	INTEGER ILLIN, IDUMA, NTREP, IDUMB, IOWNS, ICOMM, MESFLG, LUNIT
4299	DOUBLE PRECISION ROWND, ROWNS, RCOMM
4300	COMMON /ODE001/ ROWND, ROWNS(173), RCOMM(45),
4301	1 ILLIN, IDUMA(10), NTREP, IDUMB(2), IOWNS(4), ICOMM(21)
4302	COMMON /EH0001/ MESFLG, LUNIT
4303	DATA ILLIN/0/, NTREP/0/
4304	DATA MESFLG/1/, LUNIT/6/
4305	C
4306	C------------------------ END OF BLOCK DATA ----------------------------
4307	END
4308	C-----------------------------------------------------------------------
4309	C INSTRUCTIONS FOR INSTALLING THE ODESSA PACKAGE. (see @ below.)
4310	C
4311	C ODESSA is an enhanced version of the widely disseminated ODE solver
4312	C LSODE, and as such retains the same properties regarding portability.
4313	C The notes below, adapted from the installation instructions for LSODE,
4314	C are intended to facilitate the installation of the ODESSA package in
4315	C the user's software library.
4316	C
4317	C 1. Both a single and a double precision version of ODESSA are
4318	C provided in this release. It is expected that most users will
4319	C utilize the double precision version, except in the case of
4320	C extended word-length computers. Most routines used by ODESSA
4321	C are named the same regardless of whether they are single or
4322	C double precision. The exceptions are the LINPAK and BLAS
4323	C routines that follow the LINPAK/BLAS naming conventions, i.e.
4324	C D--- for a double precision routine, and S--- for a single
4325	C precision routine. Thus, care should be taken if both single
4326	C and double precision versions are stored in the same library.
4327	C
4328	C 2. Several routines in ODESSA have the same names as the LSODE
4329	C routines from which they were derived, although they contain
4330	C different code. These are: INTDY, STODE, PREPJ, SVCOM, and
4331	C RSCOM. If ODESSA is added to a subroutine library of which
4332	C LSODE is already a member, these routine names must be changed
4333	C in one of the two programs. Also see the note regarding BLOCK
4334	C DATA subroutines below.
4335	C
4336	C 3. In many cases, ODESSA uses unaltered LSODE routines and
4337	C common library routines that may already reside on your system.
4338	C The installation of ODESSA should be done so that identical routines
4339	C are shared rather than kept as duplicate copies.
4340	C a. Normally, the user calls only subroutine ODESSA, but for optional
4341	C capabilities the user may also CALL XSETUN, XSETF, SVCOM, RSCOM,
4342	C or INTDY, as described in Part II of the Full Description in the
4343	C User Documentation (ODESSA.DOC, see below). Except for INTDY,
4344	C none of these are called from within the package.
4345	C b. Two routines, EWSET and VNORM, are optionally replaceable by the
4346	C user if the package version is unsuitable. Hence, the install-
4347	C ation of the package should be done so that the user's version
4348	C for either routine overrides the package version.
4349	C c. The function routine D1MACH is provided to compute the unit
4350	C roundoff of the machine and precision in use, in a manner com-
4351	C patible with machine parameter routines developed at Bell Lab-
4352	C oratories. If such a routine has already been installed on
4353	C your system, the version supplied here may be discarded.
4354	C d. Linear algebraic systems are solved with routines from the
4355	C LINPACK collection, in conjunction with routines from the Basic
4356	C Linear Algebra module collection (BLAS). In double precision,
4357	C the names are DGEFA, DGESL, DGBFA, and DGBSL (from LINPACK), and
4358	C DAXPY, DSCAL, IDAMAX, and DDOT (from BLAS). If these routines
4359	C have already been installed on your system, copies supplied with
4360	C ODESSA may be discarded. The single precision versions of these
4361	C routines are used in the single precision version.
4362	C
4363	C 4. There are four integer variables, in the two labeled COMMON
4364	C blocks /ODE001/ and /EH0001/, which need to be loaded with DATA
4365	C statements. They can vary during execution, and are in common to
4366	C assure their retention between calls. This is legal in ANSI Fortran
4367	C only if done in a BLOCK DATA subprogram, and this package has a
4368	C BLOCK DATA for this purpose. However, BLOCK DATA subprograms can be
4369	C difficult to install in libraries, and many compilers allow such DATA
4370	C statements in subroutines. If your system allows this, the location
4371	C of the DATA statements are just after the initial type and common
4372	C declarations in subroutines ODESSA and XERR. In ODESSA, ILLIN and
4373	C NTREP are DATA-loaded as 0. In XERR, MESFLG is loaded as 1 and
4374	C LUNIT is loaded as the appropriate default logical unit number.
4375	C
4376	C 5. The ODESSA package contains subscript expressions which may not
4377	C be accepted by some compilers. Subscripts of the form I + J, I - J,
4378	C etc., occur in various routines. If any of these forms are
4379	C unacceptable to your compiler, an extra line of code setting the
4380	C subscript to a dummy integer value should be added for each subscipt.
4381	C
4382	C 6. User documentation is provided in a two-level structure
4383	C to accommmodate both the casual and serious user. The novice or
4384	C casual user should need to read only the Summary of Usage and the
4385	C Example Problem located at the beginning of the documentation. More
4386	C experienced users, requiring the full set of available options,
4387	C should read the Full Description which follows the Example Problem.
4388	C
4389	C 7. The user documentation may need corrections in the following ways:
4390	C a. If subroutine names have been changed to avoid conflicts between
4391	C the LSODE and ODESSA packages, the corresponding name changes
4392	C should be made in the documentation.
4393	C b. In the Summary of Usage, and in the description of XSETUN under
4394	C Part II of the Full Description, the default logical unit number
4395	C should be corrected if it is not 6.
4396	C c. In the Summary of Usage, users should be instructed to execute
4397	C CALL XSETF(1) before the first CALL to ODESSA, if this is neces-
4398	C sary for proper error message handling. (see note 2(e) above.)
4399	C d. In the description of the subroutines DF and JAC in the Summary
4400	C of Usage and in Part I of the Full Description, it is stated
4401	C that dummy names may be passed if these two routines are not user
4402	C supplied. Your system may require the user to supply a dummy
4403	C subroutine instead.
4404	C e. The ODESSA package treats the arguments NEQ, RTOL, and ATOL as
4405	C arrays (possibly of length 1), while the usage documentation
4406	C states that these arguments may be either arrays or scalars.
4407	C If your system does not allow such a mismatch, then the
4408	C documentation should be changed to reflect this.
4409	C 8. A demonstration program is provided with the package for
4410	C verification.
4411	C
4412	C
4413	C Jorge R. Leis and Mark A. Kramer
4414	C Department of Chemical Engineering
4415	C Massachusetts Institute of Technology
4416	C Cambridge, Massachusetts 02139
4417	C U.S.A.
4418	C
4419	C Current address of J.R. Leis (Jan. 1988):
4420	C
4421	C Shell Development Company
4422	C Westhollow Research Center
4423	C Houston, TX
4424	C
4425	C @ Adapted from 'Instructions for Installing LSODE', written by
4426	C Alan C. Hindmarsh, Mathematics & Statistics Division, L-316,
4427	C Lawrence Livermore National Laboratory, Livermore, CA. 94550

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

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int.modified_WCOPY/kpp_odessa_ddm.f @ 3534

Download in other formats: