Home

Context Navigation

atm_odessa_ddm.f

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

Download in other formats: