1 | SUBROUTINE INTEGRATE( TIN, TOUT ) |
---|
2 | |
---|
3 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
4 | INCLUDE 'KPP_ROOT_params.h' |
---|
5 | INCLUDE 'KPP_ROOT_global.h' |
---|
6 | |
---|
7 | C TIN - Start Time |
---|
8 | KPP_REAL TIN |
---|
9 | C TOUT - End Time |
---|
10 | KPP_REAL TOUT |
---|
11 | INTEGER i |
---|
12 | |
---|
13 | PARAMETER (LWORK=2*NVAR*NVAR+14*NVAR+20,LIWORK=NVAR+20) |
---|
14 | KPP_REAL WORK(LWORK) |
---|
15 | INTEGER IWORK(LIWORK) |
---|
16 | EXTERNAL FUNC_CHEM,JAC_CHEM,SOLOUT |
---|
17 | |
---|
18 | ITOL=1 ! --- VECTOR TOLERANCES |
---|
19 | IJAC=1 ! --- COMPUTE THE JACOBIAN ANALYTICALLY |
---|
20 | MLJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX |
---|
21 | MUJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX |
---|
22 | IMAS=0 ! --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM |
---|
23 | IOUT=0 ! --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION |
---|
24 | IDFX=0 ! --- INTERNAL TIME DERIVATIVE |
---|
25 | |
---|
26 | DO i=1,20 |
---|
27 | IWORK(i) = 0 |
---|
28 | WORK(i) = 0.D0 |
---|
29 | ENDDO |
---|
30 | IWORK(2) = 6 |
---|
31 | |
---|
32 | CALL KPP_ROS4(NVAR,FUNC_CHEM,Autonomous,TIN,VAR,TOUT, |
---|
33 | & STEPMIN,RTOL,ATOL,ITOL, |
---|
34 | & JAC_CHEM,IJAC,MLJAC,MUJAC,FUNC_CHEM,IDFX, |
---|
35 | & FUNC_CHEM,IMAS,MLJAC,MUJAC, |
---|
36 | & SOLOUT,IOUT, |
---|
37 | & WORK,LWORK,IWORK,LIWORK,IDID) |
---|
38 | |
---|
39 | IF (IDID.LT.0) THEN |
---|
40 | print *,'KPP_ROS4: Unsucessful step at T=', |
---|
41 | & TIN,' (IDID=',IDID,')' |
---|
42 | ENDIF |
---|
43 | |
---|
44 | RETURN |
---|
45 | END |
---|
46 | |
---|
47 | |
---|
48 | SUBROUTINE KPP_ROS4(N,FCN,IFCN,X,Y,XEND,H, |
---|
49 | & RelTol,AbsTol,ITOL, |
---|
50 | & JAC ,IJAC,MLJAC,MUJAC,DFX,IDFX, |
---|
51 | & MAS ,IMAS,MLMAS,MUMAS, |
---|
52 | & SOLOUT,IOUT, |
---|
53 | & WORK,LWORK,IWORK,LIWORK,IDID) |
---|
54 | C ---------------------------------------------------------- |
---|
55 | C NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC) |
---|
56 | C SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS MY'=F(X,Y). |
---|
57 | C THIS IS AN EMBEDDED ROSENBROCK METHOD OF ORDER (3)4 |
---|
58 | C (WITH STEP SIZE CONTROL). |
---|
59 | C C.F. SECTION IV.7 |
---|
60 | C |
---|
61 | C AUTHORS: E. HAIRER AND G. WANNER |
---|
62 | C UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES |
---|
63 | C CH-1211 GENEVE 24, SWITZERLAND |
---|
64 | C E-MAIL: HAIRER@CGEUGE51.BITNET, WANNER@CGEUGE51.BITNET |
---|
65 | C |
---|
66 | C THIS CODE IS PART OF THE BOOK: |
---|
67 | C E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL |
---|
68 | C EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. |
---|
69 | C SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS, |
---|
70 | C SPRINGER-VERLAG (1990) |
---|
71 | C |
---|
72 | C VERSION OF NOVEMBER 17, 1992 |
---|
73 | C |
---|
74 | C INPUT PARAMETERS |
---|
75 | C ---------------- |
---|
76 | C N DIMENSION OF THE SYSTEM |
---|
77 | C |
---|
78 | C FCN NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE |
---|
79 | C VALUE OF F(X,Y): |
---|
80 | C SUBROUTINE FCN(N,X,Y,F) |
---|
81 | C KPP_REAL X,Y(N),F(N) |
---|
82 | C F(1)=... ETC. |
---|
83 | C |
---|
84 | C IFCN GIVES INFORMATION ON FCN: |
---|
85 | C IFCN=0: F(X,Y) INDEPENDENT OF X (AUTONOMOUS) |
---|
86 | C IFCN=1: F(X,Y) MAY DEPEND ON X (NON-AUTONOMOUS) |
---|
87 | C |
---|
88 | C X INITIAL X-VALUE |
---|
89 | C |
---|
90 | C Y(N) INITIAL VALUES FOR Y |
---|
91 | C |
---|
92 | C XEND FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE) |
---|
93 | C |
---|
94 | C H INITIAL STEP SIZE GUESS; |
---|
95 | C FOR STIFF EQUATIONS WITH INITIAL TRANSIENT, |
---|
96 | C H=1.D0/(NORM OF F'), USUALLY 1.D-2 OR 1.D-3, IS GOOD. |
---|
97 | C THIS CHOICE IS NOT VERY IMPORTANT, THE CODE QUICKLY |
---|
98 | C ADAPTS ITS STEP SIZE. STUDY THE CHOSEN VALUES FOR A FEW |
---|
99 | C STEPS IN SUBROUTINE "SOLOUT", WHEN YOU ARE NOT SURE. |
---|
100 | C (IF H=0.D0, THE CODE PUTS H=1.D-6). |
---|
101 | C |
---|
102 | C RelTol,AbsTol RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY |
---|
103 | C CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N. |
---|
104 | C |
---|
105 | C ITOL SWITCH FOR RelTol AND AbsTol: |
---|
106 | C ITOL=0: BOTH RelTol AND AbsTol ARE SCALARS. |
---|
107 | C THE CODE KEEPS, ROUGHLY, THE LOCAL ERROR OF |
---|
108 | C Y(I) BELOW RelTol*ABS(Y(I))+AbsTol |
---|
109 | C ITOL=1: BOTH RelTol AND AbsTol ARE VECTORS. |
---|
110 | C THE CODE KEEPS THE LOCAL ERROR OF Y(I) BELOW |
---|
111 | C RelTol(I)*ABS(Y(I))+AbsTol(I). |
---|
112 | C |
---|
113 | C JAC NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES |
---|
114 | C THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO Y |
---|
115 | C (THIS ROUTINE IS ONLY CALLED IF IJAC=1; SUPPLY |
---|
116 | C A DUMMY SUBROUTINE IN THE CASE IJAC=0). |
---|
117 | C FOR IJAC=1, THIS SUBROUTINE MUST HAVE THE FORM |
---|
118 | C SUBROUTINE JAC(N,X,Y,DFY,LDFY) |
---|
119 | C KPP_REAL X,Y(N),DFY(LDFY,N) |
---|
120 | C DFY(1,1)= ... |
---|
121 | C LDFY, THE COLOMN-LENGTH OF THE ARRAY, IS |
---|
122 | C FURNISHED BY THE CALLING PROGRAM. |
---|
123 | C IF (MLJAC.EQ.N) THE JACOBIAN IS SUPPOSED TO |
---|
124 | C BE FULL AND THE PARTIAL DERIVATIVES ARE |
---|
125 | C STORED IN DFY AS |
---|
126 | C DFY(I,J) = PARTIAL F(I) / PARTIAL Y(J) |
---|
127 | C ELSE, THE JACOBIAN IS TAKEN AS BANDED AND |
---|
128 | C THE PARTIAL DERIVATIVES ARE STORED |
---|
129 | C DIAGONAL-WISE AS |
---|
130 | C DFY(I-J+MUJAC+1,J) = PARTIAL F(I) / PARTIAL Y(J). |
---|
131 | C |
---|
132 | C IJAC SWITCH FOR THE COMPUTATION OF THE JACOBIAN: |
---|
133 | C IJAC=0: JACOBIAN IS COMPUTED INTERNALLY BY FINITE |
---|
134 | C DIFFERENCES, SUBROUTINE "JAC" IS NEVER CALLED. |
---|
135 | C IJAC=1: JACOBIAN IS SUPPLIED BY SUBROUTINE JAC. |
---|
136 | C |
---|
137 | C MLJAC SWITCH FOR THE BANDED STRUCTURE OF THE JACOBIAN: |
---|
138 | C MLJAC=N: JACOBIAN IS A FULL MATRIX. THE LINEAR |
---|
139 | C ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION. |
---|
140 | C 0<=MLJAC<N: MLJAC IS THE LOWER BANDWITH OF JACOBIAN |
---|
141 | C MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW |
---|
142 | C THE MAIN DIAGONAL). |
---|
143 | C |
---|
144 | C MUJAC UPPER BANDWITH OF JACOBIAN MATRIX (>= NUMBER OF NON- |
---|
145 | C ZERO DIAGONALS ABOVE THE MAIN DIAGONAL). |
---|
146 | C NEED NOT BE DEFINED IF MLJAC=N. |
---|
147 | C |
---|
148 | C DFX NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES |
---|
149 | C THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO X |
---|
150 | C (THIS ROUTINE IS ONLY CALLED IF IDFX=1 AND IFCN=1; |
---|
151 | C SUPPLY A DUMMY SUBROUTINE IN THE CASE IDFX=0 OR IFCN=0). |
---|
152 | C OTHERWISE, THIS SUBROUTINE MUST HAVE THE FORM |
---|
153 | C SUBROUTINE DFX(N,X,Y,FX) |
---|
154 | C KPP_REAL X,Y(N),FX(N) |
---|
155 | C FX(1)= ... |
---|
156 | C |
---|
157 | C IDFX SWITCH FOR THE COMPUTATION OF THE DF/DX: |
---|
158 | C IDFX=0: DF/DX IS COMPUTED INTERNALLY BY FINITE |
---|
159 | C DIFFERENCES, SUBROUTINE "DFX" IS NEVER CALLED. |
---|
160 | C IDFX=1: DF/DX IS SUPPLIED BY SUBROUTINE DFX. |
---|
161 | C |
---|
162 | C ---- MAS,IMAS,MLMAS, AND MUMAS HAVE ANALOG MEANINGS ----- |
---|
163 | C ---- FOR THE "MASS MATRIX" (THE MATRIX "M" OF SECTION IV.8): - |
---|
164 | C |
---|
165 | C MAS NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE MASS- |
---|
166 | C MATRIX M. |
---|
167 | C IF IMAS=0, THIS MATRIX IS ASSUMED TO BE THE IDENTITY |
---|
168 | C MATRIX AND NEEDS NOT TO BE DEFINED; |
---|
169 | C SUPPLY A DUMMY SUBROUTINE IN THIS CASE. |
---|
170 | C IF IMAS=1, THE SUBROUTINE MAS IS OF THE FORM |
---|
171 | C SUBROUTINE MAS(N,AM,LMAS) |
---|
172 | C KPP_REAL AM(LMAS,N) |
---|
173 | C AM(1,1)= .... |
---|
174 | C IF (MLMAS.EQ.N) THE MASS-MATRIX IS STORED |
---|
175 | C AS FULL MATRIX LIKE |
---|
176 | C AM(I,J) = M(I,J) |
---|
177 | C ELSE, THE MATRIX IS TAKEN AS BANDED AND STORED |
---|
178 | C DIAGONAL-WISE AS |
---|
179 | C AM(I-J+MUMAS+1,J) = M(I,J). |
---|
180 | C |
---|
181 | C IMAS GIVES INFORMATION ON THE MASS-MATRIX: |
---|
182 | C IMAS=0: M IS SUPPOSED TO BE THE IDENTITY |
---|
183 | C MATRIX, MAS IS NEVER CALLED. |
---|
184 | C IMAS=1: MASS-MATRIX IS SUPPLIED. |
---|
185 | C |
---|
186 | C MLMAS SWITCH FOR THE BANDED STRUCTURE OF THE MASS-MATRIX: |
---|
187 | C MLMAS=N: THE FULL MATRIX CASE. THE LINEAR |
---|
188 | C ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION. |
---|
189 | C 0<=MLMAS<N: MLMAS IS THE LOWER BANDWITH OF THE |
---|
190 | C MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW |
---|
191 | C THE MAIN DIAGONAL). |
---|
192 | C MLMAS IS SUPPOSED TO BE .LE. MLJAC. |
---|
193 | C |
---|
194 | C MUMAS UPPER BANDWITH OF MASS-MATRIX (>= NUMBER OF NON- |
---|
195 | C ZERO DIAGONALS ABOVE THE MAIN DIAGONAL). |
---|
196 | C NEED NOT BE DEFINED IF MLMAS=N. |
---|
197 | C MUMAS IS SUPPOSED TO BE .LE. MUJAC. |
---|
198 | C |
---|
199 | C SOLOUT NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE |
---|
200 | C NUMERICAL SOLUTION DURING INTEGRATION. |
---|
201 | C IF IOUT=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP. |
---|
202 | C SUPPLY A DUMMY SUBROUTINE IF IOUT=0. |
---|
203 | C IT MUST HAVE THE FORM |
---|
204 | C SUBROUTINE SOLOUT (NR,XOLD,X,Y,N,IRTRN) |
---|
205 | C KPP_REAL X,Y(N) |
---|
206 | C .... |
---|
207 | C SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH |
---|
208 | C GRID-POINT "X" (THEREBY THE INITIAL VALUE IS |
---|
209 | C THE FIRST GRID-POINT). |
---|
210 | C "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN |
---|
211 | C IS SET <0, ROS4 RETURNS TO THE CALLING PROGRAM. |
---|
212 | C |
---|
213 | C IOUT GIVES INFORMATION ON THE SUBROUTINE SOLOUT: |
---|
214 | C IOUT=0: SUBROUTINE IS NEVER CALLED |
---|
215 | C IOUT=1: SUBROUTINE IS USED FOR OUTPUT |
---|
216 | C |
---|
217 | C WORK ARRAY OF WORKING SPACE OF LENGTH "LWORK". |
---|
218 | C SERVES AS WORKING SPACE FOR ALL VECTORS AND MATRICES. |
---|
219 | C "LWORK" MUST BE AT LEAST |
---|
220 | C N*(LJAC+LMAS+LE1+8)+5 |
---|
221 | C WHERE |
---|
222 | C LJAC=N IF MLJAC=N (FULL JACOBIAN) |
---|
223 | C LJAC=MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.) |
---|
224 | C AND |
---|
225 | C LMAS=0 IF IMAS=0 |
---|
226 | C LMAS=N IF IMAS=1 AND MLMAS=N (FULL) |
---|
227 | C LMAS=MLMAS+MUMAS+1 IF MLMAS<N (BANDED MASS-M.) |
---|
228 | C AND |
---|
229 | C LE1=N IF MLJAC=N (FULL JACOBIAN) |
---|
230 | C LE1=2*MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.). |
---|
231 | c |
---|
232 | C IN THE USUAL CASE WHERE THE JACOBIAN IS FULL AND THE |
---|
233 | C MASS-MATRIX IS THE INDENTITY (IMAS=0), THE MINIMUM |
---|
234 | C STORAGE REQUIREMENT IS |
---|
235 | C LWORK = 2*N*N+8*N+5. |
---|
236 | C |
---|
237 | C LWORK DECLARED LENGHT OF ARRAY "WORK". |
---|
238 | C |
---|
239 | C IWORK INTEGER WORKING SPACE OF LENGTH "LIWORK". |
---|
240 | C "LIWORK" MUST BE AT LEAST N+2. |
---|
241 | C |
---|
242 | C LIWORK DECLARED LENGHT OF ARRAY "IWORK". |
---|
243 | C |
---|
244 | C ---------------------------------------------------------------------- |
---|
245 | C |
---|
246 | C SOPHISTICATED SETTING OF PARAMETERS |
---|
247 | C ----------------------------------- |
---|
248 | C SEVERAL PARAMETERS OF THE CODE ARE TUNED TO MAKE IT WORK |
---|
249 | C WELL. THEY MAY BE DEFINED BY SETTING WORK(1),..,WORK(5) |
---|
250 | C AS WELL AS IWORK(1),IWORK(2) DIFFERENT FROM ZERO. |
---|
251 | C FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES: |
---|
252 | C |
---|
253 | C IWORK(1) THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS. |
---|
254 | C THE DEFAULT VALUE (FOR IWORK(1)=0) IS 100000. |
---|
255 | C |
---|
256 | C IWORK(2) SWITCH FOR THE CHOICE OF THE COEFFICIENTS |
---|
257 | C IF IWORK(2).EQ.1 METHOD OF SHAMPINE |
---|
258 | C IF IWORK(2).EQ.2 METHOD GRK4T OF KAPS-RENTROP |
---|
259 | C IF IWORK(2).EQ.3 METHOD GRK4A OF KAPS-RENTROP |
---|
260 | C IF IWORK(2).EQ.4 METHOD OF VAN VELDHUIZEN (GAMMA=1/2) |
---|
261 | C IF IWORK(2).EQ.5 METHOD OF VAN VELDHUIZEN ("D-STABLE") |
---|
262 | C IF IWORK(2).EQ.6 AN L-STABLE METHOD |
---|
263 | C THE DEFAULT VALUE (FOR IWORK(2)=0) IS IWORK(2)=2. |
---|
264 | C |
---|
265 | C WORK(1) UROUND, THE ROUNDING UNIT, DEFAULT 1.D-16. |
---|
266 | C |
---|
267 | C WORK(2) MAXIMAL STEP SIZE, DEFAULT XEND-X. |
---|
268 | C |
---|
269 | C WORK(3), WORK(4) PARAMETERS FOR STEP SIZE SELECTION |
---|
270 | C THE NEW STEP SIZE IS CHOSEN SUBJECT TO THE RESTRICTION |
---|
271 | C WORK(3) <= HNEW/HOLD <= WORK(4) |
---|
272 | C DEFAULT VALUES: WORK(3)=0.2D0, WORK(4)=6.D0 |
---|
273 | C |
---|
274 | C WORK(5) AVOID THE HUMP: AFTER TWO CONSECUTIVE STEP REJECTIONS |
---|
275 | C THE STEP SIZE IS MULTIPLIED BY WORK(5) |
---|
276 | C DEFAULT VALUES: WORK(5)=0.1D0 |
---|
277 | C |
---|
278 | C----------------------------------------------------------------------- |
---|
279 | C |
---|
280 | C OUTPUT PARAMETERS |
---|
281 | C ----------------- |
---|
282 | C X X-VALUE WHERE THE SOLUTION IS COMPUTED |
---|
283 | C (AFTER SUCCESSFUL RETURN X=XEND) |
---|
284 | C |
---|
285 | C Y(N) SOLUTION AT X |
---|
286 | C |
---|
287 | C H PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP |
---|
288 | C |
---|
289 | C IDID REPORTS ON SUCCESSFULNESS UPON RETURN: |
---|
290 | C IDID=1 COMPUTATION SUCCESSFUL, |
---|
291 | C IDID=-1 COMPUTATION UNSUCCESSFUL. |
---|
292 | C |
---|
293 | C --------------------------------------------------------- |
---|
294 | C *** *** *** *** *** *** *** *** *** *** *** *** *** |
---|
295 | C DECLARATIONS |
---|
296 | C *** *** *** *** *** *** *** *** *** *** *** *** *** |
---|
297 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
298 | DIMENSION Y(N),AbsTol(1),RelTol(1),WORK(LWORK),IWORK(LIWORK) |
---|
299 | LOGICAL AUTNMS,IMPLCT,JBAND,ARRET |
---|
300 | EXTERNAL FCN,JAC,DFX,MAS,SOLOUT |
---|
301 | COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL |
---|
302 | C ----------------------------------------------------- |
---|
303 | C --- COMMON STAT CAN BE USED FOR STATISTICS |
---|
304 | C --- NFCN NUMBER OF FUNC_CHEMCTION EVALUATIONS (THOSE FOR NUMERICAL |
---|
305 | C EVALUATION OF THE JACOBIAN ARE NOT COUNTED) |
---|
306 | C --- NJAC NUMBER OF JACOBIAN EVALUATIONS (EITHER ANALYTICALLY |
---|
307 | C OR NUMERICALLY) |
---|
308 | C --- NSTEP NUMBER OF COMPUTED STEPS |
---|
309 | C --- NACCPT NUMBER OF ACCEPTED STEPS |
---|
310 | C --- NREJCT NUMBER OF REJECTED STEPS (AFTER AT LEAST ONE STEP |
---|
311 | C HAS BEEN ACCEPTED) |
---|
312 | C --- NDEC NUMBER OF LU-DECOMPOSITIONS (N-DIMENSIONAL MATRIX) |
---|
313 | C --- NSOL NUMBER OF FORWARD-BACKWARD SUBSTITUTIONS |
---|
314 | C ----------------------------------------------------- |
---|
315 | C *** *** *** *** *** *** *** |
---|
316 | C SETTING THE PARAMETERS |
---|
317 | C *** *** *** *** *** *** *** |
---|
318 | NFCN=0 |
---|
319 | NJAC=0 |
---|
320 | NSTEP=0 |
---|
321 | NACCPT=0 |
---|
322 | NREJCT=0 |
---|
323 | NDEC=0 |
---|
324 | NSOL=0 |
---|
325 | ARRET=.FALSE. |
---|
326 | C -------- NMAX , THE MAXIMAL NUMBER OF STEPS ----- |
---|
327 | IF(IWORK(1).EQ.0)THEN |
---|
328 | NMAX=100000 |
---|
329 | ELSE |
---|
330 | NMAX=IWORK(1) |
---|
331 | IF(NMAX.LE.0)THEN |
---|
332 | WRITE(6,*)' WRONG INPUT IWORK(1)=',IWORK(1) |
---|
333 | ARRET=.TRUE. |
---|
334 | END IF |
---|
335 | END IF |
---|
336 | C -------- METH COEFFICIENTS OF THE METHOD |
---|
337 | IF(IWORK(2).EQ.0)THEN |
---|
338 | METH=2 |
---|
339 | ELSE |
---|
340 | METH=IWORK(2) |
---|
341 | IF(METH.LE.0.OR.METH.GE.7)THEN |
---|
342 | WRITE(6,*)' CURIOUS INPUT IWORK(2)=',IWORK(2) |
---|
343 | ARRET=.TRUE. |
---|
344 | END IF |
---|
345 | END IF |
---|
346 | C -------- UROUND SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.D0 |
---|
347 | IF(WORK(1).EQ.0.D0)THEN |
---|
348 | UROUND=1.D-16 |
---|
349 | ELSE |
---|
350 | UROUND=WORK(1) |
---|
351 | IF(UROUND.LE.1.D-14.OR.UROUND.GE.1.D0)THEN |
---|
352 | WRITE(6,*)' COEFFICIENTS HAVE 16 DIGITS, UROUND=',WORK(1) |
---|
353 | ARRET=.TRUE. |
---|
354 | END IF |
---|
355 | END IF |
---|
356 | C -------- MAXIMAL STEP SIZE |
---|
357 | IF(WORK(2).EQ.0.D0)THEN |
---|
358 | HMAX=XEND-X |
---|
359 | ELSE |
---|
360 | HMAX=WORK(2) |
---|
361 | END IF |
---|
362 | C ------- FAC1,FAC2 PARAMETERS FOR STEP SIZE SELECTION |
---|
363 | IF(WORK(3).EQ.0.D0)THEN |
---|
364 | FAC1=5.D0 |
---|
365 | ELSE |
---|
366 | FAC1=1.D0/WORK(3) |
---|
367 | END IF |
---|
368 | IF(WORK(4).EQ.0.D0)THEN |
---|
369 | FAC2=1.D0/6.0D0 |
---|
370 | ELSE |
---|
371 | FAC2=1.D0/WORK(4) |
---|
372 | END IF |
---|
373 | C ------- FACREJ FOR THE HUMP |
---|
374 | IF(WORK(5).EQ.0.D0)THEN |
---|
375 | FACREJ=0.1D0 |
---|
376 | ELSE |
---|
377 | FACREJ=WORK(5) |
---|
378 | END IF |
---|
379 | C --------- CHECK IF TOLERANCES ARE O.K. |
---|
380 | IF (ITOL.EQ.0) THEN |
---|
381 | IF (AbsTol(1).LE.0.D0.OR.RelTol(1).LE.10.D0*UROUND) THEN |
---|
382 | WRITE (6,*) ' TOLERANCES ARE TOO SMALL' |
---|
383 | ARRET=.TRUE. |
---|
384 | END IF |
---|
385 | ELSE |
---|
386 | DO 15 I=1,N |
---|
387 | IF (AbsTol(I).LE.0.D0.OR.RelTol(I).LE.10.D0*UROUND) THEN |
---|
388 | WRITE (6,*) ' TOLERANCES(',I,') ARE TOO SMALL' |
---|
389 | ARRET=.TRUE. |
---|
390 | END IF |
---|
391 | 15 CONTINUE |
---|
392 | END IF |
---|
393 | C *** *** *** *** *** *** *** *** *** *** *** *** *** |
---|
394 | C COMPUTATION OF ARRAY ENTRIES |
---|
395 | C *** *** *** *** *** *** *** *** *** *** *** *** *** |
---|
396 | C ---- AUTONOMOUS, IMPLICIT, BANDED OR NOT ? |
---|
397 | AUTNMS=IFCN.EQ.0 |
---|
398 | IMPLCT=IMAS.NE.0 |
---|
399 | JBAND=MLJAC.NE.N |
---|
400 | ARRET=.FALSE. |
---|
401 | C -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS --- |
---|
402 | C -- JACOBIAN |
---|
403 | IF(JBAND)THEN |
---|
404 | LDJAC=MLJAC+MUJAC+1 |
---|
405 | ELSE |
---|
406 | LDJAC=N |
---|
407 | END IF |
---|
408 | C -- MATRIX E FOR LINEAR ALGEBRA |
---|
409 | IF(JBAND)THEN |
---|
410 | LDE=2*MLJAC+MUJAC+1 |
---|
411 | ELSE |
---|
412 | LDE=N |
---|
413 | END IF |
---|
414 | C -- MASS MATRIX |
---|
415 | IF (IMPLCT) THEN |
---|
416 | IF (MLMAS.NE.N) THEN |
---|
417 | LDMAS=MLMAS+MUMAS+1 |
---|
418 | ELSE |
---|
419 | LDMAS=N |
---|
420 | END IF |
---|
421 | C ------ BANDWITH OF "MAS" NOT LARGER THAN BANDWITH OF "JAC" |
---|
422 | IF (MLMAS.GT.MLJAC.OR.MUMAS.GT.MUJAC) THEN |
---|
423 | WRITE (6,*) 'BANDWITH OF "MAS" NOT LARGER THAN BANDWITH OF |
---|
424 | & "JAC"' |
---|
425 | ARRET=.TRUE. |
---|
426 | END IF |
---|
427 | ELSE |
---|
428 | LDMAS=0 |
---|
429 | END IF |
---|
430 | LDMAS2=MAX(1,LDMAS) |
---|
431 | C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK ----- |
---|
432 | IEYNEW=6 |
---|
433 | IEDY1=IEYNEW+N |
---|
434 | IEDY=IEDY1+N |
---|
435 | IEAK1=IEDY+N |
---|
436 | IEAK2=IEAK1+N |
---|
437 | IEAK3=IEAK2+N |
---|
438 | IEAK4=IEAK3+N |
---|
439 | IEFX =IEAK4+N |
---|
440 | IEJAC=IEFX +N |
---|
441 | IEMAS=IEJAC+N*LDJAC |
---|
442 | IEE =IEMAS+N*LDMAS |
---|
443 | C ------ TOTAL STORAGE REQUIREMENT ----------- |
---|
444 | ISTORE=IEE+N*LDE-1 |
---|
445 | IF(ISTORE.GT.LWORK)THEN |
---|
446 | WRITE(6,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE |
---|
447 | ARRET=.TRUE. |
---|
448 | END IF |
---|
449 | C ------- ENTRY POINTS FOR INTEGER WORKSPACE ----- |
---|
450 | IEIP=3 |
---|
451 | C --------- TOTAL REQUIREMENT --------------- |
---|
452 | ISTORE=IEIP+N-1 |
---|
453 | IF(ISTORE.GT.LIWORK)THEN |
---|
454 | WRITE(6,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE |
---|
455 | ARRET=.TRUE. |
---|
456 | END IF |
---|
457 | C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1 |
---|
458 | IF (ARRET) THEN |
---|
459 | IDID=-1 |
---|
460 | RETURN |
---|
461 | END IF |
---|
462 | C -------- CALL TO CORE INTEGRATOR ------------ |
---|
463 | CALL RO4COR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,ITOL,JAC,IJAC, |
---|
464 | & MLJAC,MUJAC,DFX,IDFX,MAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID, |
---|
465 | & NMAX,UROUND,METH,FAC1,FAC2,FACREJ,AUTNMS,IMPLCT,JBAND, |
---|
466 | & LDJAC,LDE,LDMAS2,WORK(IEYNEW),WORK(IEDY1),WORK(IEDY), |
---|
467 | & WORK(IEAK1),WORK(IEAK2),WORK(IEAK3),WORK(IEAK4), |
---|
468 | & WORK(IEFX),WORK(IEJAC),WORK(IEE),WORK(IEMAS),IWORK(IEIP)) |
---|
469 | C ----------- RETURN ----------- |
---|
470 | RETURN |
---|
471 | END |
---|
472 | C |
---|
473 | C |
---|
474 | C |
---|
475 | C ----- ... AND HERE IS THE CORE INTEGRATOR ---------- |
---|
476 | C |
---|
477 | SUBROUTINE RO4COR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,ITOL,JAC, |
---|
478 | & IJAC,MLJAC,MUJAC,DFX,IDFX,MAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID, |
---|
479 | & NMAX,UROUND,METH,FAC1,FAC2,FACREJ,AUTNMS,IMPLCT,BANDED, |
---|
480 | & LFJAC,LE,LDMAS,YNEW,DY1,DY,AK1,AK2,AK3,AK4,FX,FJAC,E,FMAS,IP) |
---|
481 | C ---------------------------------------------------------- |
---|
482 | C CORE INTEGRATOR FOR ROS4 |
---|
483 | C PARAMETERS SAME AS IN ROS4 WITH WORKSPACE ADDED |
---|
484 | C ---------------------------------------------------------- |
---|
485 | C DECLARATIONS |
---|
486 | C ---------------------------------------------------------- |
---|
487 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
488 | INCLUDE 'KPP_ROOT_params.h' |
---|
489 | INCLUDE 'KPP_ROOT_sparse.h' |
---|
490 | KPP_REAL Y(N),YNEW(N),DY1(N),DY(N),AK1(N), |
---|
491 | * AK2(N),AK3(N),AK4(N),FX(N), |
---|
492 | * FJAC(LU_NONZERO),E(LU_NONZERO), |
---|
493 | * FMAS(LDMAS,N),AbsTol(1),RelTol(1) |
---|
494 | INTEGER IP(N) |
---|
495 | LOGICAL REJECT,RJECT2,AUTNMS,IMPLCT,BANDED |
---|
496 | COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL |
---|
497 | EXTERNAL FCN,JAC,MAS,SOLOUT,DFX |
---|
498 | |
---|
499 | |
---|
500 | C ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ---------- |
---|
501 | IF (IMPLCT) CALL MAS(N,FMAS,LDMAS) |
---|
502 | C ---- PREPARE BANDWIDTHS ----- |
---|
503 | IF (BANDED) THEN |
---|
504 | MLE=MLJAC |
---|
505 | MUE=MUJAC |
---|
506 | MBJAC=MLJAC+MUJAC+1 |
---|
507 | MBB=MLMAS+MUMAS+1 |
---|
508 | MDIAG=MLE+MUE+1 |
---|
509 | MBDIAG=MUMAS+1 |
---|
510 | MDIFF=MLE+MUE-MUMAS |
---|
511 | END IF |
---|
512 | C *** *** *** *** *** *** *** |
---|
513 | C INITIALISATIONS |
---|
514 | C *** *** *** *** *** *** *** |
---|
515 | IF (METH.EQ.1) CALL SHAMP (A21,A31,A32,C21,C31,C32,C41,C42,C43, |
---|
516 | & B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) |
---|
517 | IF (METH.EQ.2) CALL GRK4T (A21,A31,A32,C21,C31,C32,C41,C42,C43, |
---|
518 | & B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) |
---|
519 | IF (METH.EQ.3) CALL GRK4A (A21,A31,A32,C21,C31,C32,C41,C42,C43, |
---|
520 | & B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) |
---|
521 | IF (METH.EQ.4) CALL VELDS (A21,A31,A32,C21,C31,C32,C41,C42,C43, |
---|
522 | & B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) |
---|
523 | IF (METH.EQ.5) CALL VELDD (A21,A31,A32,C21,C31,C32,C41,C42,C43, |
---|
524 | & B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) |
---|
525 | IF (METH.EQ.6) CALL LSTAB (A21,A31,A32,C21,C31,C32,C41,C42,C43, |
---|
526 | & B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) |
---|
527 | C --- INITIAL PREPARATIONS |
---|
528 | POSNEG=SIGN(1.D0,XEND-X) |
---|
529 | HMAXN=MIN(ABS(HMAX),ABS(XEND-X)) |
---|
530 | IF (ABS(H).LE.10.D0*UROUND) H=1.0D-6 |
---|
531 | H=MIN(ABS(H),HMAXN) |
---|
532 | H=SIGN(H,POSNEG) |
---|
533 | REJECT=.FALSE. |
---|
534 | NSING=0 |
---|
535 | IRTRN=1 |
---|
536 | XXOLD=X |
---|
537 | IF (IRTRN.LT.0) GOTO 79 |
---|
538 | C --- BASIC INTEGRATION STEP |
---|
539 | 1 IF (NSTEP.GT.NMAX.OR.X+.1D0*H.EQ.X.OR.ABS(H).LE.UROUND) GOTO 79 |
---|
540 | IF ((X-XEND)*POSNEG+UROUND.GT.0.D0) THEN |
---|
541 | H=HOPT |
---|
542 | IDID=1 |
---|
543 | RETURN |
---|
544 | END IF |
---|
545 | HOPT=H |
---|
546 | IF ((X+H-XEND)*POSNEG.GT.0.D0) H=XEND-X |
---|
547 | CALL FCN(N,X,Y,DY1) |
---|
548 | NFCN=NFCN+1 |
---|
549 | |
---|
550 | C *** *** *** *** *** *** *** |
---|
551 | C COMPUTATION OF THE JACOBIAN |
---|
552 | C *** *** *** *** *** *** *** |
---|
553 | NJAC=NJAC+1 |
---|
554 | CALL JAC(N,X,Y,FJAC) |
---|
555 | |
---|
556 | IF (.NOT.AUTNMS) THEN |
---|
557 | IF (IDFX.EQ.0) THEN |
---|
558 | C --- COMPUTE NUMERICALLY THE DERIVATIVE WITH RESPECT TO X |
---|
559 | DELT=DSQRT(UROUND*MAX(1.D-5,ABS(X))) |
---|
560 | XDELT=X+DELT |
---|
561 | CALL FCN(N,XDELT,Y,AK1) |
---|
562 | DO 19 J=1,N |
---|
563 | 19 FX(J)=(AK1(J)-DY1(J))/DELT |
---|
564 | ELSE |
---|
565 | C --- COMPUTE ANALYTICALLY THE DERIVATIVE WITH RESPECT TO X |
---|
566 | CALL DFX(N,X,Y,FX) |
---|
567 | END IF |
---|
568 | END IF |
---|
569 | 2 CONTINUE |
---|
570 | |
---|
571 | C *** *** *** *** *** *** *** |
---|
572 | C COMPUTE THE STAGES |
---|
573 | C *** *** *** *** *** *** *** |
---|
574 | NDEC=NDEC+1 |
---|
575 | HC21=C21/H |
---|
576 | HC31=C31/H |
---|
577 | HC32=C32/H |
---|
578 | HC41=C41/H |
---|
579 | HC42=C42/H |
---|
580 | HC43=C43/H |
---|
581 | FAC=1.D0/(H*GAMMA) |
---|
582 | C --- THE MATRIX E (B=IDENTITY, JACOBIAN A FULL MATRIX) |
---|
583 | DO 800 I=1,LU_NONZERO |
---|
584 | 800 E(I)=-FJAC(I) |
---|
585 | DO 801 J=1,N |
---|
586 | 801 E(LU_DIAG(J))=E(LU_DIAG(J))+FAC |
---|
587 | CALL KppDecomp (E,IER) |
---|
588 | IF (IER.NE.0) GOTO 80 |
---|
589 | IF (AUTNMS) THEN |
---|
590 | C --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE |
---|
591 | C --- 1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM |
---|
592 | C --- 2) THE JACOBIAN OF THE PROBLEM IS A FULL MATRIX |
---|
593 | C --- 3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS |
---|
594 | DO 803 I=1,N |
---|
595 | 803 AK1(I)=DY1(I) |
---|
596 | CALL KppSolve(E,AK1) |
---|
597 | DO 810 I=1,N |
---|
598 | 810 YNEW(I)=Y(I)+A21*AK1(I) |
---|
599 | CALL FCN(N,X,YNEW,DY) |
---|
600 | DO 811 I=1,N |
---|
601 | 811 AK2(I)=DY(I)+HC21*AK1(I) |
---|
602 | CALL KppSolve(E,AK2) |
---|
603 | DO 820 I=1,N |
---|
604 | 820 YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I) |
---|
605 | CALL FCN(N,X,YNEW,DY) |
---|
606 | DO 821 I=1,N |
---|
607 | 821 AK3(I)=DY(I)+HC31*AK1(I)+HC32*AK2(I) |
---|
608 | CALL KppSolve(E,AK3) |
---|
609 | DO 831 I=1,N |
---|
610 | 831 AK4(I)=DY(I)+HC41*AK1(I)+HC42*AK2(I)+HC43*AK3(I) |
---|
611 | CALL KppSolve(E,AK4) |
---|
612 | ELSE |
---|
613 | C --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE |
---|
614 | C --- 1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM |
---|
615 | C --- 2) THE JACOBIAN OF THE PROBLEM IS A FULL MATRIX |
---|
616 | C --- 3) THE DIFFERENTIAL EQUATION IS NON-AUTONOMOUS |
---|
617 | HD1=H*D1 |
---|
618 | HD2=H*D2 |
---|
619 | HD3=H*D3 |
---|
620 | HD4=H*D4 |
---|
621 | DO 903 I=1,N |
---|
622 | 903 AK1(I)=DY1(I)+HD1*FX(I) |
---|
623 | CALL KppSolve(E,AK1) |
---|
624 | DO 910 I=1,N |
---|
625 | 910 YNEW(I)=Y(I)+A21*AK1(I) |
---|
626 | CALL FCN(N,X+C2*H,YNEW,DY) |
---|
627 | DO 911 I=1,N |
---|
628 | 911 AK2(I)=DY(I)+HD2*FX(I)+HC21*AK1(I) |
---|
629 | CALL KppSolve(E,AK2) |
---|
630 | DO 920 I=1,N |
---|
631 | 920 YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I) |
---|
632 | CALL FCN(N,X+C3*H,YNEW,DY) |
---|
633 | DO 921 I=1,N |
---|
634 | 921 AK3(I)=DY(I)+HD3*FX(I)+HC31*AK1(I)+HC32*AK2(I) |
---|
635 | CALL KppSolve(E,AK3) |
---|
636 | DO 931 I=1,N |
---|
637 | 931 AK4(I)=DY(I)+HD4*FX(I)+HC41*AK1(I)+HC42*AK2(I) |
---|
638 | & +HC43*AK3(I) |
---|
639 | CALL KppSolve(E,AK4) |
---|
640 | END IF |
---|
641 | NSOL=NSOL+4 |
---|
642 | NFCN=NFCN+2 |
---|
643 | C *** *** *** *** *** *** *** |
---|
644 | C ERROR ESTIMATION |
---|
645 | C *** *** *** *** *** *** *** |
---|
646 | NSTEP=NSTEP+1 |
---|
647 | C ------------ NEW SOLUTION --------------- |
---|
648 | DO 240 I=1,N |
---|
649 | 240 YNEW(I)=Y(I)+B1*AK1(I)+B2*AK2(I)+B3*AK3(I)+B4*AK4(I) |
---|
650 | C ------------ COMPUTE ERROR ESTIMATION ---------------- |
---|
651 | ERR=0.D0 |
---|
652 | DO 300 I=1,N |
---|
653 | S=E1*AK1(I)+E2*AK2(I)+E3*AK3(I)+E4*AK4(I) |
---|
654 | IF (ITOL.EQ.0) THEN |
---|
655 | SK=AbsTol(1)+RelTol(1)*DMAX1(DABS(Y(I)),DABS(YNEW(I))) |
---|
656 | ELSE |
---|
657 | SK=AbsTol(I)+RelTol(I)*DMAX1(DABS(Y(I)),DABS(YNEW(I))) |
---|
658 | END IF |
---|
659 | 300 ERR=ERR+(S/SK)**2 |
---|
660 | ERR=DSQRT(ERR/N) |
---|
661 | C --- COMPUTATION OF HNEW |
---|
662 | C --- WE REQUIRE .2<=HNEW/H<=6. |
---|
663 | FAC=DMAX1(FAC2,DMIN1(FAC1,(ERR)**.25D0/.9D0)) |
---|
664 | HNEW=H/FAC |
---|
665 | C *** *** *** *** *** *** *** |
---|
666 | C IS THE ERROR SMALL ENOUGH ? |
---|
667 | C *** *** *** *** *** *** *** |
---|
668 | IF (ERR.LE.1.D0) THEN |
---|
669 | C --- STEP IS ACCEPTED |
---|
670 | NACCPT=NACCPT+1 |
---|
671 | DO 44 I=1,N |
---|
672 | 44 Y(I)=YNEW(I) |
---|
673 | XXOLD=X |
---|
674 | X=X+H |
---|
675 | IF (IRTRN.LT.0) GOTO 79 |
---|
676 | IF (DABS(HNEW).GT.HMAXN) HNEW=POSNEG*HMAXN |
---|
677 | IF (REJECT) HNEW=POSNEG*DMIN1(DABS(HNEW),DABS(H)) |
---|
678 | REJECT=.FALSE. |
---|
679 | RJECT2=.FALSE. |
---|
680 | H=HNEW |
---|
681 | GOTO 1 |
---|
682 | ELSE |
---|
683 | C --- STEP IS REJECTED |
---|
684 | IF (RJECT2) HNEW=H*FACREJ |
---|
685 | IF (REJECT) RJECT2=.TRUE. |
---|
686 | REJECT=.TRUE. |
---|
687 | H=HNEW |
---|
688 | IF (NACCPT.GE.1) NREJCT=NREJCT+1 |
---|
689 | GOTO 2 |
---|
690 | END IF |
---|
691 | C --- EXIT |
---|
692 | 80 WRITE (6,*) ' MATRIX E IS SINGULAR, INFO = ',INFO |
---|
693 | NSING=NSING+1 |
---|
694 | IF (NSING.GE.5) GOTO 79 |
---|
695 | H=H*0.5D0 |
---|
696 | GOTO 2 |
---|
697 | 79 WRITE(6,979)X,H |
---|
698 | 979 FORMAT(' EXIT OF ROS4 AT X=',D16.7,' H=',D16.7) |
---|
699 | IDID=-1 |
---|
700 | RETURN |
---|
701 | END |
---|
702 | C |
---|
703 | SUBROUTINE SHAMP (A21,A31,A32,C21,C31,C32,C41,C42,C43, |
---|
704 | & B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) |
---|
705 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
706 | A21=2.D0 |
---|
707 | A31=48.D0/25.D0 |
---|
708 | A32=6.D0/25.D0 |
---|
709 | C21=-8.D0 |
---|
710 | C31=372.D0/25.D0 |
---|
711 | C32=12.D0/5.D0 |
---|
712 | C41=-112.D0/125.D0 |
---|
713 | C42=-54.D0/125.D0 |
---|
714 | C43=-2.D0/5.D0 |
---|
715 | B1=19.D0/9.D0 |
---|
716 | B2=1.D0/2.D0 |
---|
717 | B3=25.D0/108.D0 |
---|
718 | B4=125.D0/108.D0 |
---|
719 | E1=17.D0/54.D0 |
---|
720 | E2=7.D0/36.D0 |
---|
721 | E3=0.D0 |
---|
722 | E4=125.D0/108.D0 |
---|
723 | GAMMA=.5D0 |
---|
724 | C2= 0.1000000000000000D+01 |
---|
725 | C3= 0.6000000000000000D+00 |
---|
726 | D1= 0.5000000000000000D+00 |
---|
727 | D2=-0.1500000000000000D+01 |
---|
728 | D3= 0.2420000000000000D+01 |
---|
729 | D4= 0.1160000000000000D+00 |
---|
730 | RETURN |
---|
731 | END |
---|
732 | C |
---|
733 | SUBROUTINE GRK4A (A21,A31,A32,C21,C31,C32,C41,C42,C43, |
---|
734 | & B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) |
---|
735 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
736 | A21= 0.1108860759493671D+01 |
---|
737 | A31= 0.2377085261983360D+01 |
---|
738 | A32= 0.1850114988899692D+00 |
---|
739 | C21=-0.4920188402397641D+01 |
---|
740 | C31= 0.1055588686048583D+01 |
---|
741 | C32= 0.3351817267668938D+01 |
---|
742 | C41= 0.3846869007049313D+01 |
---|
743 | C42= 0.3427109241268180D+01 |
---|
744 | C43=-0.2162408848753263D+01 |
---|
745 | B1= 0.1845683240405840D+01 |
---|
746 | B2= 0.1369796894360503D+00 |
---|
747 | B3= 0.7129097783291559D+00 |
---|
748 | B4= 0.6329113924050632D+00 |
---|
749 | E1= 0.4831870177201765D-01 |
---|
750 | E2=-0.6471108651049505D+00 |
---|
751 | E3= 0.2186876660500240D+00 |
---|
752 | E4=-0.6329113924050632D+00 |
---|
753 | GAMMA= 0.3950000000000000D+00 |
---|
754 | C2= 0.4380000000000000D+00 |
---|
755 | C3= 0.8700000000000000D+00 |
---|
756 | D1= 0.3950000000000000D+00 |
---|
757 | D2=-0.3726723954840920D+00 |
---|
758 | D3= 0.6629196544571492D-01 |
---|
759 | D4= 0.4340946962568634D+00 |
---|
760 | RETURN |
---|
761 | END |
---|
762 | C |
---|
763 | SUBROUTINE GRK4T (A21,A31,A32,C21,C31,C32,C41,C42,C43, |
---|
764 | & B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) |
---|
765 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
766 | A21= 0.2000000000000000D+01 |
---|
767 | A31= 0.4524708207373116D+01 |
---|
768 | A32= 0.4163528788597648D+01 |
---|
769 | C21=-0.5071675338776316D+01 |
---|
770 | C31= 0.6020152728650786D+01 |
---|
771 | C32= 0.1597506846727117D+00 |
---|
772 | C41=-0.1856343618686113D+01 |
---|
773 | C42=-0.8505380858179826D+01 |
---|
774 | C43=-0.2084075136023187D+01 |
---|
775 | B1= 0.3957503746640777D+01 |
---|
776 | B2= 0.4624892388363313D+01 |
---|
777 | B3= 0.6174772638750108D+00 |
---|
778 | B4= 0.1282612945269037D+01 |
---|
779 | E1= 0.2302155402932996D+01 |
---|
780 | E2= 0.3073634485392623D+01 |
---|
781 | E3=-0.8732808018045032D+00 |
---|
782 | E4=-0.1282612945269037D+01 |
---|
783 | GAMMA= 0.2310000000000000D+00 |
---|
784 | C2= 0.4620000000000000D+00 |
---|
785 | C3= 0.8802083333333334D+00 |
---|
786 | D1= 0.2310000000000000D+00 |
---|
787 | D2=-0.3962966775244303D-01 |
---|
788 | D3= 0.5507789395789127D+00 |
---|
789 | D4=-0.5535098457052764D-01 |
---|
790 | RETURN |
---|
791 | END |
---|
792 | C |
---|
793 | SUBROUTINE VELDS (A21,A31,A32,C21,C31,C32,C41,C42,C43, |
---|
794 | & B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) |
---|
795 | C --- METHOD GIVEN BY VAN VELDHUIZEN |
---|
796 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
797 | A21= 0.2000000000000000D+01 |
---|
798 | A31= 0.1750000000000000D+01 |
---|
799 | A32= 0.2500000000000000D+00 |
---|
800 | C21=-0.8000000000000000D+01 |
---|
801 | C31=-0.8000000000000000D+01 |
---|
802 | C32=-0.1000000000000000D+01 |
---|
803 | C41= 0.5000000000000000D+00 |
---|
804 | C42=-0.5000000000000000D+00 |
---|
805 | C43= 0.2000000000000000D+01 |
---|
806 | B1= 0.1333333333333333D+01 |
---|
807 | B2= 0.6666666666666667D+00 |
---|
808 | B3=-0.1333333333333333D+01 |
---|
809 | B4= 0.1333333333333333D+01 |
---|
810 | E1=-0.3333333333333333D+00 |
---|
811 | E2=-0.3333333333333333D+00 |
---|
812 | E3=-0.0000000000000000D+00 |
---|
813 | E4=-0.1333333333333333D+01 |
---|
814 | GAMMA= 0.5000000000000000D+00 |
---|
815 | C2= 0.1000000000000000D+01 |
---|
816 | C3= 0.5000000000000000D+00 |
---|
817 | D1= 0.5000000000000000D+00 |
---|
818 | D2=-0.1500000000000000D+01 |
---|
819 | D3=-0.7500000000000000D+00 |
---|
820 | D4= 0.2500000000000000D+00 |
---|
821 | RETURN |
---|
822 | END |
---|
823 | C |
---|
824 | SUBROUTINE VELDD (A21,A31,A32,C21,C31,C32,C41,C42,C43, |
---|
825 | & B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) |
---|
826 | C --- METHOD GIVEN BY VAN VELDHUIZEN |
---|
827 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
828 | A21= 0.2000000000000000D+01 |
---|
829 | A31= 0.4812234362695436D+01 |
---|
830 | A32= 0.4578146956747842D+01 |
---|
831 | C21=-0.5333333333333331D+01 |
---|
832 | C31= 0.6100529678848254D+01 |
---|
833 | C32= 0.1804736797378427D+01 |
---|
834 | C41=-0.2540515456634749D+01 |
---|
835 | C42=-0.9443746328915205D+01 |
---|
836 | C43=-0.1988471753215993D+01 |
---|
837 | B1= 0.4289339254654537D+01 |
---|
838 | B2= 0.5036098482851414D+01 |
---|
839 | B3= 0.6085736420673917D+00 |
---|
840 | B4= 0.1355958941201148D+01 |
---|
841 | E1= 0.2175672787531755D+01 |
---|
842 | E2= 0.2950911222575741D+01 |
---|
843 | E3=-0.7859744544887430D+00 |
---|
844 | E4=-0.1355958941201148D+01 |
---|
845 | GAMMA= 0.2257081148225682D+00 |
---|
846 | C2= 0.4514162296451364D+00 |
---|
847 | C3= 0.8755928946018455D+00 |
---|
848 | D1= 0.2257081148225682D+00 |
---|
849 | D2=-0.4599403502680582D-01 |
---|
850 | D3= 0.5177590504944076D+00 |
---|
851 | D4=-0.3805623938054428D-01 |
---|
852 | RETURN |
---|
853 | END |
---|
854 | C |
---|
855 | SUBROUTINE LSTAB (A21,A31,A32,C21,C31,C32,C41,C42,C43, |
---|
856 | & B1,B2,B3,B4,E1,E2,E3,E4,GAMMA,C2,C3,D1,D2,D3,D4) |
---|
857 | C --- AN L-STABLE METHOD |
---|
858 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
859 | A21= 0.2000000000000000D+01 |
---|
860 | A31= 0.1867943637803922D+01 |
---|
861 | A32= 0.2344449711399156D+00 |
---|
862 | C21=-0.7137615036412310D+01 |
---|
863 | C31= 0.2580708087951457D+01 |
---|
864 | C32= 0.6515950076447975D+00 |
---|
865 | C41=-0.2137148994382534D+01 |
---|
866 | C42=-0.3214669691237626D+00 |
---|
867 | C43=-0.6949742501781779D+00 |
---|
868 | C B_i = M_i |
---|
869 | B1= 0.2255570073418735D+01 |
---|
870 | B2= 0.2870493262186792D+00 |
---|
871 | B3= 0.4353179431840180D+00 |
---|
872 | B4= 0.1093502252409163D+01 |
---|
873 | C E_i = error estimator |
---|
874 | E1=-0.2815431932141155D+00 |
---|
875 | E2=-0.7276199124938920D-01 |
---|
876 | E3=-0.1082196201495311D+00 |
---|
877 | E4=-0.1093502252409163D+01 |
---|
878 | C gamma = gamma |
---|
879 | GAMMA= 0.5728200000000000D+00 |
---|
880 | C C_i = alpha_i |
---|
881 | C2= 0.1145640000000000D+01 |
---|
882 | C3= 0.6552168638155900D+00 |
---|
883 | C D_i = gamma_i |
---|
884 | D1= 0.5728200000000000D+00 |
---|
885 | D2=-0.1769193891319233D+01 |
---|
886 | D3= 0.7592633437920482D+00 |
---|
887 | D4=-0.1049021087100450D+00 |
---|
888 | RETURN |
---|
889 | END |
---|
890 | |
---|
891 | SUBROUTINE SOLOUT (NR,XOLD,X,Y,N,IRTRN) |
---|
892 | C --- PRINTS SOLUTION |
---|
893 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
894 | DIMENSION Y(N) |
---|
895 | COMMON /INTERN/XOUT |
---|
896 | IF (NR.EQ.1) THEN |
---|
897 | WRITE (6,99) X,Y(1),Y(2),NR-1 |
---|
898 | XOUT=0.1D0 |
---|
899 | ELSE |
---|
900 | IF (X.GE.XOUT) THEN |
---|
901 | WRITE (6,99) X,Y(1),Y(2),NR-1 |
---|
902 | XOUT=DMAX1(XOUT+0.1D0,X) |
---|
903 | END IF |
---|
904 | END IF |
---|
905 | 99 FORMAT(1X,'X =',F5.2,' Y =',2E18.10,' NSTEP =',I4) |
---|
906 | RETURN |
---|
907 | END |
---|
908 | |
---|
909 | SUBROUTINE DEC (N, NDIM, A, IP, IER) |
---|
910 | C VERSION REAL KPP_REAL |
---|
911 | INTEGER N,NDIM,IP,IER,NM1,K,KP1,M,I,J |
---|
912 | KPP_REAL A,T |
---|
913 | DIMENSION A(NDIM,N), IP(N) |
---|
914 | C----------------------------------------------------------------------- |
---|
915 | C MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION. |
---|
916 | C INPUT.. |
---|
917 | C N = ORDER OF MATRIX. |
---|
918 | C NDIM = DECLARED DIMENSION OF ARRAY A . |
---|
919 | C A = MATRIX TO BE TRIANGULARIZED. |
---|
920 | C OUTPUT.. |
---|
921 | C A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U . |
---|
922 | C A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L. |
---|
923 | C IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW. |
---|
924 | C IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O . |
---|
925 | C IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE |
---|
926 | C SINGULAR AT STAGE K. |
---|
927 | C USE SOL TO OBTAIN SOLUTION OF LINEAR SYSTEM. |
---|
928 | C DETERM(A) = IP(N)*A(1,1)*A(2,2)*...*A(N,N). |
---|
929 | C IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO. |
---|
930 | C |
---|
931 | C REFERENCE.. |
---|
932 | C C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR, |
---|
933 | C C.A.C.M. 15 (1972), P. 274. |
---|
934 | C----------------------------------------------------------------------- |
---|
935 | IER = 0 |
---|
936 | IP(N) = 1 |
---|
937 | IF (N .EQ. 1) GO TO 70 |
---|
938 | NM1 = N - 1 |
---|
939 | DO 60 K = 1,NM1 |
---|
940 | KP1 = K + 1 |
---|
941 | M = K |
---|
942 | DO 10 I = KP1,N |
---|
943 | IF (DABS(A(I,K)) .GT. DABS(A(M,K))) M = I |
---|
944 | 10 CONTINUE |
---|
945 | IP(K) = M |
---|
946 | T = A(M,K) |
---|
947 | IF (M .EQ. K) GO TO 20 |
---|
948 | IP(N) = -IP(N) |
---|
949 | A(M,K) = A(K,K) |
---|
950 | A(K,K) = T |
---|
951 | 20 CONTINUE |
---|
952 | IF (T .EQ. 0.D0) GO TO 80 |
---|
953 | T = 1.D0/T |
---|
954 | DO 30 I = KP1,N |
---|
955 | 30 A(I,K) = -A(I,K)*T |
---|
956 | DO 50 J = KP1,N |
---|
957 | T = A(M,J) |
---|
958 | A(M,J) = A(K,J) |
---|
959 | A(K,J) = T |
---|
960 | IF (T .EQ. 0.D0) GO TO 45 |
---|
961 | DO 40 I = KP1,N |
---|
962 | 40 A(I,J) = A(I,J) + A(I,K)*T |
---|
963 | 45 CONTINUE |
---|
964 | 50 CONTINUE |
---|
965 | 60 CONTINUE |
---|
966 | 70 K = N |
---|
967 | IF (A(N,N) .EQ. 0.D0) GO TO 80 |
---|
968 | RETURN |
---|
969 | 80 IER = K |
---|
970 | IP(N) = 0 |
---|
971 | RETURN |
---|
972 | C----------------------- END OF SUBROUTINE DEC ------------------------- |
---|
973 | END |
---|
974 | C |
---|
975 | C |
---|
976 | SUBROUTINE SOL (N, NDIM, A, B, IP) |
---|
977 | C VERSION REAL KPP_REAL |
---|
978 | INTEGER N,NDIM,IP,NM1,K,KP1,M,I,KB,KM1 |
---|
979 | KPP_REAL A,B,T |
---|
980 | DIMENSION A(NDIM,N), B(N), IP(N) |
---|
981 | C----------------------------------------------------------------------- |
---|
982 | C SOLUTION OF LINEAR SYSTEM, A*X = B . |
---|
983 | C INPUT.. |
---|
984 | C N = ORDER OF MATRIX. |
---|
985 | C NDIM = DECLARED DIMENSION OF ARRAY A . |
---|
986 | C A = TRIANGULARIZED MATRIX OBTAINED FROM DEC. |
---|
987 | C B = RIGHT HAND SIDE VECTOR. |
---|
988 | C IP = PIVOT VECTOR OBTAINED FROM DEC. |
---|
989 | C DO NOT USE IF DEC HAS SET IER .NE. 0. |
---|
990 | C OUTPUT.. |
---|
991 | C B = SOLUTION VECTOR, X . |
---|
992 | C----------------------------------------------------------------------- |
---|
993 | IF (N .EQ. 1) GO TO 50 |
---|
994 | NM1 = N - 1 |
---|
995 | DO 20 K = 1,NM1 |
---|
996 | KP1 = K + 1 |
---|
997 | M = IP(K) |
---|
998 | T = B(M) |
---|
999 | B(M) = B(K) |
---|
1000 | B(K) = T |
---|
1001 | DO 10 I = KP1,N |
---|
1002 | 10 B(I) = B(I) + A(I,K)*T |
---|
1003 | 20 CONTINUE |
---|
1004 | DO 40 KB = 1,NM1 |
---|
1005 | KM1 = N - KB |
---|
1006 | K = KM1 + 1 |
---|
1007 | B(K) = B(K)/A(K,K) |
---|
1008 | T = -B(K) |
---|
1009 | DO 30 I = 1,KM1 |
---|
1010 | 30 B(I) = B(I) + A(I,K)*T |
---|
1011 | 40 CONTINUE |
---|
1012 | 50 B(1) = B(1)/A(1,1) |
---|
1013 | RETURN |
---|
1014 | C----------------------- END OF SUBROUTINE SOL ------------------------- |
---|
1015 | END |
---|
1016 | |
---|
1017 | |
---|
1018 | |
---|
1019 | |
---|
1020 | SUBROUTINE FUNC_CHEM(N, T, Y, P) |
---|
1021 | INCLUDE 'KPP_ROOT_params.h' |
---|
1022 | INCLUDE 'KPP_ROOT_global.h' |
---|
1023 | INTEGER N |
---|
1024 | KPP_REAL T, Told |
---|
1025 | KPP_REAL Y(NVAR), P(NVAR) |
---|
1026 | Told = TIME |
---|
1027 | TIME = T |
---|
1028 | CALL Update_SUN() |
---|
1029 | CALL Update_RCONST() |
---|
1030 | CALL Fun( Y, FIX, RCONST, P ) |
---|
1031 | TIME = Told |
---|
1032 | RETURN |
---|
1033 | END |
---|
1034 | |
---|
1035 | |
---|
1036 | SUBROUTINE JAC_CHEM(N, T, Y, J) |
---|
1037 | INCLUDE 'KPP_ROOT_params.h' |
---|
1038 | INCLUDE 'KPP_ROOT_global.h' |
---|
1039 | INTEGER N |
---|
1040 | KPP_REAL Told, T |
---|
1041 | KPP_REAL Y(NVAR), J(LU_NONZERO) |
---|
1042 | Told = TIME |
---|
1043 | TIME = T |
---|
1044 | CALL Update_SUN() |
---|
1045 | CALL Update_RCONST() |
---|
1046 | CALL Jac_SP( Y, FIX, RCONST, J ) |
---|
1047 | TIME = Told |
---|
1048 | RETURN |
---|
1049 | END |
---|
1050 | |
---|
1051 | |
---|
1052 | |
---|