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 |
---|
17 | |
---|
18 | |
---|
19 | ITOL=1 ! --- VECTOR TOLERANCES |
---|
20 | IJAC=1 ! --- COMPUTE THE JACOBIAN ANALYTICALLY |
---|
21 | MLJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX |
---|
22 | MUJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX |
---|
23 | IMAS=0 ! --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM |
---|
24 | IOUT=0 ! --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION |
---|
25 | IDFX=0 ! --- INTERNAL TIME DERIVATIVE |
---|
26 | |
---|
27 | DO i=1,20 |
---|
28 | IWORK(i) = 0 |
---|
29 | WORK(i) = 0.D0 |
---|
30 | ENDDO |
---|
31 | |
---|
32 | IWORK(3) = 1 |
---|
33 | |
---|
34 | CALL ATMRODAS(NVAR,FUNC_CHEM,Autonomous,TIN,VAR,TOUT, |
---|
35 | & STEPMIN,RTOL,ATOL,ITOL, |
---|
36 | & JAC_CHEM,IJAC,MLJAC,MUJAC,FUNC_CHEM,IDFX, |
---|
37 | & FUNC_CHEM,IMAS, |
---|
38 | & WORK,LWORK,IWORK,LIWORK,IDID) |
---|
39 | |
---|
40 | IF (IDID.LT.0) THEN |
---|
41 | print *,'ATMRODAS: Unsucessfull exit at T=', |
---|
42 | & TIN,' (IDID=',IDID,')' |
---|
43 | ENDIF |
---|
44 | |
---|
45 | |
---|
46 | RETURN |
---|
47 | END |
---|
48 | |
---|
49 | |
---|
50 | SUBROUTINE ATMRODAS(N,FCN,IFCN,X,Y,XEND,H, |
---|
51 | & RelTol,AbsTol,ITOL, |
---|
52 | & JAC ,IJAC,MLJAC,MUJAC,DFX,IDFX, |
---|
53 | & MAS ,IMAS, |
---|
54 | & WORK,LWORK,IWORK,LIWORK,IDID) |
---|
55 | C ---------------------------------------------------------- |
---|
56 | C NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC) |
---|
57 | C SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS MY'=F(X,Y). |
---|
58 | C THIS IS AN EMBEDDED ROSENBROCK METHOD OF ORDER (3)4 |
---|
59 | C (WITH STEP SIZE CONTROL). |
---|
60 | C C.F. SECTIONS IV.7 AND VI.3 |
---|
61 | C |
---|
62 | C AUTHORS: E. HAIRER AND G. WANNER |
---|
63 | C UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES |
---|
64 | C CH-1211 GENEVE 24, SWITZERLAND |
---|
65 | C E-MAIL: HAIRER@DIVSUN.UNIGE.CH, WANNER@DIVSUN.UNIGE.CH |
---|
66 | C --------------------------------------------------------- |
---|
67 | C *** *** *** *** *** *** *** *** *** *** *** *** *** |
---|
68 | C DECLARATIONS |
---|
69 | C *** *** *** *** *** *** *** *** *** *** *** *** *** |
---|
70 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
71 | DIMENSION Y(N),AbsTol(*),RelTol(*),WORK(LWORK),IWORK(LIWORK) |
---|
72 | LOGICAL AUTNMS,IMPLCT,JBAND,ARRET,PRED |
---|
73 | EXTERNAL FCN,JAC,DFX,MAS |
---|
74 | COMMON /STATISTICS/ NFCN,NACCPT,NREJCT,NSTEP,NJAC,NDEC,NSOL |
---|
75 | C *** *** *** *** *** *** *** |
---|
76 | C SETTING THE PARAMETERS |
---|
77 | C *** *** *** *** *** *** *** |
---|
78 | ARRET=.FALSE. |
---|
79 | METH = 1 |
---|
80 | NMAX=100000 |
---|
81 | C -------- PRED STEP SIZE CONTROL |
---|
82 | IF(IWORK(3).LE.1)THEN |
---|
83 | PRED=.TRUE. |
---|
84 | ELSE |
---|
85 | PRED=.FALSE. |
---|
86 | END IF |
---|
87 | UROUND=1.D-16 |
---|
88 | NM1 = N |
---|
89 | M1 = N |
---|
90 | M2 = N |
---|
91 | C -------- MAXIMAL STEP SIZE |
---|
92 | IF(WORK(2).EQ.0.D0)THEN |
---|
93 | HMAX=XEND-X |
---|
94 | ELSE |
---|
95 | HMAX=WORK(2) |
---|
96 | END IF |
---|
97 | C ------- FAC1,FAC2 PARAMETERS FOR STEP SIZE SELECTION |
---|
98 | IF(WORK(3).EQ.0.D0)THEN |
---|
99 | FAC1=5.D0 |
---|
100 | ELSE |
---|
101 | FAC1=1.D0/WORK(3) |
---|
102 | END IF |
---|
103 | IF(WORK(4).EQ.0.D0)THEN |
---|
104 | FAC2=1.D0/6.0D0 |
---|
105 | ELSE |
---|
106 | FAC2=1.D0/WORK(4) |
---|
107 | END IF |
---|
108 | IF (FAC1.LT.1.0D0.OR.FAC2.GT.1.0D0) THEN |
---|
109 | WRITE(6,*)' CURIOUS INPUT WORK(3,4)=',WORK(3),WORK(4) |
---|
110 | ARRET=.TRUE. |
---|
111 | END IF |
---|
112 | C --------- SAFE SAFETY FACTOR IN STEP SIZE PREDICTION |
---|
113 | SAFE=0.9D0 |
---|
114 | C --------- CHECK IF TOLERANCES ARE O.K. |
---|
115 | IF (ITOL.EQ.0) THEN |
---|
116 | IF (AbsTol(1).LE.0.D0.OR.RelTol(1).LE.10.D0*UROUND) THEN |
---|
117 | WRITE (6,*) ' TOLERANCES ARE TOO SMALL' |
---|
118 | ARRET=.TRUE. |
---|
119 | END IF |
---|
120 | ELSE |
---|
121 | DO I=1,N |
---|
122 | IF (AbsTol(I).LE.0.D0.OR.RelTol(I).LE.10.D0*UROUND) THEN |
---|
123 | WRITE (6,*) ' TOLERANCES(',I,') ARE TOO SMALL' |
---|
124 | ARRET=.TRUE. |
---|
125 | END IF |
---|
126 | END DO |
---|
127 | END IF |
---|
128 | |
---|
129 | IF (ARRET) STOP |
---|
130 | NM1 = N |
---|
131 | C *** *** *** *** *** *** *** *** *** *** *** *** *** |
---|
132 | C COMPUTATION OF ARRAY ENTRIES |
---|
133 | C *** *** *** *** *** *** *** *** *** *** *** *** *** |
---|
134 | C ---- AUTONOMOUS, IMPLICIT, BANDED OR NOT ? |
---|
135 | AUTNMS=IFCN.EQ.0 |
---|
136 | IMPLCT=IMAS.NE.0 |
---|
137 | JBAND=MLJAC.LT.NM1 |
---|
138 | C -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS --- |
---|
139 | C -- JACOBIAN AND MATRIX E |
---|
140 | MLJAC=NM1 |
---|
141 | MUJAC=NM1 |
---|
142 | LDJAC=NM1 |
---|
143 | LDE=NM1 |
---|
144 | C -- MASS MATRIX |
---|
145 | IF (IMPLCT) THEN |
---|
146 | print *, 'Implicit 1' |
---|
147 | ELSE |
---|
148 | LDMAS=0 |
---|
149 | IJOB=1 |
---|
150 | END IF |
---|
151 | LDMAS2=MAX(1,LDMAS) |
---|
152 | C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK ----- |
---|
153 | IEYNEW=21 |
---|
154 | IEDY1=IEYNEW+N |
---|
155 | IEDY=IEDY1+N |
---|
156 | IEAK1=IEDY+N |
---|
157 | IEAK2=IEAK1+N |
---|
158 | IEAK3=IEAK2+N |
---|
159 | IEAK4=IEAK3+N |
---|
160 | IEAK5=IEAK4+N |
---|
161 | IEAK6=IEAK5+N |
---|
162 | IEFX =IEAK6+N |
---|
163 | IECON=IEFX+N |
---|
164 | IEJAC=IECON+4*N |
---|
165 | IEMAS=IEJAC+N*LDJAC |
---|
166 | IEE =IEMAS+NM1*LDMAS |
---|
167 | C ------ TOTAL STORAGE REQUIREMENT ----------- |
---|
168 | ISTORE=IEE+NM1*LDE-1 |
---|
169 | IF(ISTORE.GT.LWORK)THEN |
---|
170 | WRITE(6,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE |
---|
171 | ARRET=.TRUE. |
---|
172 | END IF |
---|
173 | C ------- ENTRY POINTS FOR INTEGER WORKSPACE ----- |
---|
174 | IEIP=21 |
---|
175 | ISTORE=IEIP+NM1-1 |
---|
176 | IF(ISTORE.GT.LIWORK)THEN |
---|
177 | WRITE(6,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE |
---|
178 | ARRET=.TRUE. |
---|
179 | END IF |
---|
180 | C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1 |
---|
181 | IF (ARRET) THEN |
---|
182 | IDID=-1 |
---|
183 | RETURN |
---|
184 | END IF |
---|
185 | C -------- CALL TO CORE INTEGRATOR ------------ |
---|
186 | CALL ROSCOR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,ITOL,JAC,IJAC, |
---|
187 | & MLJAC,MUJAC,DFX,IDFX,MAS,MLMAS,MUMAS,IOUT,IDID,NMAX, |
---|
188 | & UROUND,METH,IJOB,FAC1,FAC2,SAFE,AUTNMS,IMPLCT,JBAND,PRED,LDJAC, |
---|
189 | & LDE,LDMAS2,WORK(IEYNEW),WORK(IEDY1),WORK(IEDY),WORK(IEAK1), |
---|
190 | & WORK(IEAK2),WORK(IEAK3),WORK(IEAK4),WORK(IEAK5),WORK(IEAK6), |
---|
191 | & WORK(IEFX),WORK(IEJAC),WORK(IEE),WORK(IEMAS),IWORK(IEIP), |
---|
192 | & WORK(IECON), |
---|
193 | & M1,M2,NM1,NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL) |
---|
194 | IWORK(14)=NFCN |
---|
195 | IWORK(15)=NJAC |
---|
196 | IWORK(16)=NSTEP |
---|
197 | IWORK(17)=NACCPT |
---|
198 | IWORK(18)=NREJCT |
---|
199 | IWORK(19)=NDEC |
---|
200 | IWORK(20)=NSOL |
---|
201 | C ----------- RETURN ----------- |
---|
202 | RETURN |
---|
203 | END |
---|
204 | C |
---|
205 | C |
---|
206 | C |
---|
207 | C ----- ... AND HERE IS THE CORE INTEGRATOR ---------- |
---|
208 | C |
---|
209 | SUBROUTINE ROSCOR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol, |
---|
210 | & ITOL,JAC,IJAC, |
---|
211 | & MLJAC,MUJAC,DFX,IDFX,MAS,MLMAS,MUMAS,IOUT,IDID,NMAX, |
---|
212 | & UROUND,METH,IJOB,FAC1,FAC2,SAFE,AUTNMS,IMPLCT,BANDED, |
---|
213 | & PRED,LDJAC, |
---|
214 | & LDE,LDMAS,YNEW,DY1,DY,AK1,AK2,AK3,AK4,AK5,AK6, |
---|
215 | & FX,FJAC,E,FMAS,IP,CONT, |
---|
216 | & M1,M2,NM1,NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL) |
---|
217 | C ---------------------------------------------------------- |
---|
218 | C CORE INTEGRATOR FOR RODAS |
---|
219 | C PARAMETERS SAME AS IN RODAS WITH WORKSPACE ADDED |
---|
220 | C ---------------------------------------------------------- |
---|
221 | C DECLARATIONS |
---|
222 | C ---------------------------------------------------------- |
---|
223 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
224 | INCLUDE 'KPP_ROOT_params.h' |
---|
225 | INCLUDE 'KPP_ROOT_sparse.h' |
---|
226 | DIMENSION Y(N),YNEW(N),DY1(N),DY(N),AK1(N), |
---|
227 | * AK2(N),AK3(N),AK4(N),AK5(N),AK6(N),FX(N), |
---|
228 | * FJAC(LU_NONZERO),E(LDE,NM1),FMAS(LDMAS,NM1), |
---|
229 | * AbsTol(*),RelTol(*) |
---|
230 | DIMENSION CONT(4*N) |
---|
231 | INTEGER IP(NM1) |
---|
232 | LOGICAL REJECT,AUTNMS,IMPLCT,BANDED |
---|
233 | LOGICAL ONE,LAST,PRED,SINGULAR |
---|
234 | EXTERNAL FCN, MAS, JAC, DFX |
---|
235 | COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG |
---|
236 | COMMON /CONROS/XOLD,HOUT,NN |
---|
237 | C *** *** *** *** *** *** *** |
---|
238 | C INITIALISATIONS |
---|
239 | C *** *** *** *** *** *** *** |
---|
240 | NN=N |
---|
241 | NN2=2*N |
---|
242 | NN3=3*N |
---|
243 | LRC=4*N |
---|
244 | C ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ---------- |
---|
245 | IF (IMPLCT) CALL MAS (NM1,FMAS,LDMAS) |
---|
246 | C ------ SET THE PARAMETERS OF THE METHOD ----- |
---|
247 | CALL ROCOE(METH,A21,A31,A32,A41,A42,A43,A51,A52,A53,A54, |
---|
248 | & C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,C61, |
---|
249 | & C62,C63,C64,C65,GAMMA,C2,C3,C4,D1,D2,D3,D4, |
---|
250 | & D21,D22,D23,D24,D25,D31,D32,D33,D34,D35) |
---|
251 | C --- INITIAL PREPARATIONS |
---|
252 | IF (M1.GT.0) IJOB=IJOB+10 |
---|
253 | POSNEG=SIGN(1.D0,XEND-X) |
---|
254 | HMAXN=DMIN1(DABS(HMAX),DABS(XEND-X)) |
---|
255 | IF (DABS(H).LE.10.D0*UROUND) H=1.0D-6 |
---|
256 | H=DMIN1(DABS(H),HMAXN) |
---|
257 | H=SIGN(H,POSNEG) |
---|
258 | HACC = H |
---|
259 | ERRACC = 1.0d0 |
---|
260 | REJECT=.FALSE. |
---|
261 | LAST=.FALSE. |
---|
262 | NSING=0 |
---|
263 | IRTRN=1 |
---|
264 | IF (AUTNMS) THEN |
---|
265 | HD1=0.0D0 |
---|
266 | HD2=0.0D0 |
---|
267 | HD3=0.0D0 |
---|
268 | HD4=0.0D0 |
---|
269 | END IF |
---|
270 | C -------- PREPARE BAND-WIDTHS -------- |
---|
271 | MBDIAG=MUMAS+1 |
---|
272 | |
---|
273 | C --- BASIC INTEGRATION STEP |
---|
274 | LAST = .FALSE. |
---|
275 | DO WHILE (.NOT.LAST) |
---|
276 | IF (.NOT. REJECT) THEN |
---|
277 | IF (NSTEP.GT.NMAX) CALL FAIL_EXIT(3,X,IDID,H,NMAX) |
---|
278 | IF ( 0.1D0*DABS(H) .LE. DABS(X)*UROUND ) |
---|
279 | * CALL FAIL_EXIT(2,X,IDID,H,NMAX) |
---|
280 | HOPT=H |
---|
281 | IF ((X+H*1.0001D0-XEND)*POSNEG.GE.0.D0) THEN |
---|
282 | H=XEND-X |
---|
283 | LAST=.TRUE. |
---|
284 | END IF |
---|
285 | C *** *** *** *** *** *** *** |
---|
286 | C COMPUTATION OF THE JACOBIAN |
---|
287 | C *** *** *** *** *** *** *** |
---|
288 | CALL FCN(N,X,Y,DY1) |
---|
289 | CALL JAC(N,X,Y,FJAC,LDJAC) |
---|
290 | NFCN=NFCN+1 |
---|
291 | NJAC=NJAC+1 |
---|
292 | |
---|
293 | IF (.NOT.AUTNMS) THEN |
---|
294 | C --- COMPUTE NUMERICALLY THE DERIVATIVE WITH RESPECT TO X |
---|
295 | DELT=DSQRT(UROUND*DMAX1(1.D-5,DABS(X))) |
---|
296 | XDELT=X+DELT |
---|
297 | CALL FCN(N,XDELT,Y,AK1) |
---|
298 | DO J=1,N |
---|
299 | FX(J)=(AK1(J)-DY1(J))/DELT |
---|
300 | END DO |
---|
301 | END IF |
---|
302 | END IF |
---|
303 | |
---|
304 | C *** *** *** *** *** *** *** |
---|
305 | C COMPUTE THE STAGES |
---|
306 | C *** *** *** *** *** *** *** |
---|
307 | SINGULAR = .TRUE. |
---|
308 | DO WHILE (SINGULAR) |
---|
309 | FAC=1.D0/(H*GAMMA) |
---|
310 | CALL DECOMR(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
311 | & M1,M2,NM1,FAC,E,LDE,IP,IER,IJOB,IMPLCT,IP) |
---|
312 | SINGULAR = IER.NE.0 |
---|
313 | IF (SINGULAR) THEN |
---|
314 | NSING=NSING+1 |
---|
315 | IF (NSING.GE.5) CALL FAIL_EXIT(1,X,IDID,H,NMAX) |
---|
316 | H=H*0.5D0 |
---|
317 | REJECT=.TRUE. |
---|
318 | LAST=.FALSE. |
---|
319 | ONE = .FALSE. |
---|
320 | END IF |
---|
321 | END DO |
---|
322 | |
---|
323 | NDEC=NDEC+1 |
---|
324 | C --- PREPARE FOR THE COMPUTATION OF THE 6 STAGES |
---|
325 | HC21=C21/H |
---|
326 | HC31=C31/H |
---|
327 | HC32=C32/H |
---|
328 | HC41=C41/H |
---|
329 | HC42=C42/H |
---|
330 | HC43=C43/H |
---|
331 | HC51=C51/H |
---|
332 | HC52=C52/H |
---|
333 | HC53=C53/H |
---|
334 | HC54=C54/H |
---|
335 | HC61=C61/H |
---|
336 | HC62=C62/H |
---|
337 | HC63=C63/H |
---|
338 | HC64=C64/H |
---|
339 | HC65=C65/H |
---|
340 | IF (.NOT.AUTNMS) THEN |
---|
341 | HD1=H*D1 |
---|
342 | HD2=H*D2 |
---|
343 | HD3=H*D3 |
---|
344 | HD4=H*D4 |
---|
345 | END IF |
---|
346 | C --- THE STAGES |
---|
347 | CALL SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
348 | & M1,M2,NM1,FAC,E,LDE,IP,DY1,AK1,FX,YNEW,HD1,IJOB,.FALSE.) |
---|
349 | DO I=1,N |
---|
350 | YNEW(I)=Y(I)+A21*AK1(I) |
---|
351 | END DO |
---|
352 | CALL FCN(N,X+C2*H,YNEW,DY) |
---|
353 | DO I=1,N |
---|
354 | YNEW(I)=HC21*AK1(I) |
---|
355 | END DO |
---|
356 | CALL SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
357 | & M1,M2,NM1,FAC,E,LDE,IP,DY,AK2,FX,YNEW,HD2,IJOB,.TRUE.) |
---|
358 | DO I=1,N |
---|
359 | YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I) |
---|
360 | END DO |
---|
361 | CALL FCN(N,X+C3*H,YNEW,DY) |
---|
362 | DO I=1,N |
---|
363 | YNEW(I)=HC31*AK1(I)+HC32*AK2(I) |
---|
364 | END DO |
---|
365 | CALL SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
366 | & M1,M2,NM1,FAC,E,LDE,IP,DY,AK3,FX,YNEW,HD3,IJOB,.TRUE.) |
---|
367 | DO I=1,N |
---|
368 | YNEW(I)=Y(I)+A41*AK1(I)+A42*AK2(I)+A43*AK3(I) |
---|
369 | END DO |
---|
370 | CALL FCN(N,X+C4*H,YNEW,DY) |
---|
371 | DO I=1,N |
---|
372 | YNEW(I)=HC41*AK1(I)+HC42*AK2(I)+HC43*AK3(I) |
---|
373 | END DO |
---|
374 | CALL SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
375 | & M1,M2,NM1,FAC,E,LDE,IP,DY,AK4,FX,YNEW,HD4,IJOB,.TRUE.) |
---|
376 | DO I=1,N |
---|
377 | YNEW(I)=Y(I)+A51*AK1(I)+A52*AK2(I)+A53*AK3(I)+A54*AK4(I) |
---|
378 | END DO |
---|
379 | CALL FCN(N,X+H,YNEW,DY) |
---|
380 | DO I=1,N |
---|
381 | AK6(I)=HC52*AK2(I)+HC54*AK4(I)+HC51*AK1(I)+HC53*AK3(I) |
---|
382 | END DO |
---|
383 | CALL SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
384 | & M1,M2,NM1,FAC,E,LDE,IP,DY,AK5,FX,AK6,0.D0,IJOB,.TRUE.) |
---|
385 | C ------------ EMBEDDED SOLUTION --------------- |
---|
386 | DO I=1,N |
---|
387 | YNEW(I)=YNEW(I)+AK5(I) |
---|
388 | END DO |
---|
389 | CALL FCN(N,X+H,YNEW,DY) |
---|
390 | DO I=1,N |
---|
391 | AK5(I)=HC61*AK1(I)+HC62*AK2(I)+HC65*AK5(I) |
---|
392 | & +HC64*AK4(I)+HC63*AK3(I) |
---|
393 | END DO |
---|
394 | CALL SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
395 | & M1,M2,NM1,FAC,E,LDE,IP,DY,AK6,FX,AK5,0.D0,IJOB,.TRUE.) |
---|
396 | C ------------ NEW SOLUTION --------------- |
---|
397 | DO I=1,N |
---|
398 | YNEW(I)=YNEW(I)+AK6(I) |
---|
399 | END DO |
---|
400 | NSOL=NSOL+6 |
---|
401 | NFCN=NFCN+5 |
---|
402 | |
---|
403 | C *** *** *** *** *** *** *** |
---|
404 | C ERROR ESTIMATION |
---|
405 | C *** *** *** *** *** *** *** |
---|
406 | NSTEP=NSTEP+1 |
---|
407 | C ------------ COMPUTE ERROR ESTIMATION ---------------- |
---|
408 | ERR=0.0D0 |
---|
409 | DO I=1,N |
---|
410 | IF (ITOL.EQ.0) THEN |
---|
411 | SK=AbsTol(1)+RelTol(1)*DMAX1(DABS(Y(I)),DABS(YNEW(I))) |
---|
412 | ELSE |
---|
413 | SK=AbsTol(I)+RelTol(I)*DMAX1(DABS(Y(I)),DABS(YNEW(I))) |
---|
414 | END IF |
---|
415 | ERR=ERR+(AK6(I)/SK)**2 |
---|
416 | c2 ERR = DMAX1(ERR, AK6(I)/SK) |
---|
417 | END DO |
---|
418 | ERR=DSQRT(ERR/N) |
---|
419 | |
---|
420 | C --- COMPUTATION OF HNEW |
---|
421 | C --- WE REQUIRE .2<=HNEW/H<=6. |
---|
422 | FAC=DMAX1(FAC2,DMIN1(FAC1,(ERR)**0.25D0/SAFE)) |
---|
423 | HNEW=DMAX1(H/FAC, STEPMIN) |
---|
424 | |
---|
425 | C *** *** *** *** *** *** *** |
---|
426 | C IS THE ERROR SMALL ENOUGH ? |
---|
427 | C *** *** *** *** *** *** *** |
---|
428 | |
---|
429 | IF ( (ERR.LE.1.D0).or.(H.LE.STEPMIN) ) THEN |
---|
430 | C --- STEP IS ACCEPTED |
---|
431 | NACCPT=NACCPT+1 |
---|
432 | IF (PRED) THEN |
---|
433 | C --- PREDICTIVE CONTROLLER OF GUSTAFSSON |
---|
434 | IF (NACCPT.GT.1) THEN |
---|
435 | FACGUS=(HACC/H)*(ERR**2/ERRACC)**0.25D0/SAFE |
---|
436 | FACGUS=DMAX1(FAC2,DMIN1(FAC1,FACGUS)) |
---|
437 | FAC=DMAX1(FAC,FACGUS) |
---|
438 | HNEW=DMAX1(H/FAC, STEPMIN) |
---|
439 | END IF |
---|
440 | HACC=H |
---|
441 | ERRACC=DMAX1(1.0D-2,ERR) |
---|
442 | END IF |
---|
443 | DO I=1,N |
---|
444 | Y(I)=YNEW(I) |
---|
445 | END DO |
---|
446 | XOLD=X |
---|
447 | X=X+H |
---|
448 | IF (DABS(HNEW).GT.HMAXN) HNEW=POSNEG*HMAXN |
---|
449 | IF (REJECT) HNEW=POSNEG*DMIN1(DABS(HNEW),DABS(H)) |
---|
450 | REJECT=.FALSE. |
---|
451 | H=HNEW |
---|
452 | ELSE |
---|
453 | C --- STEP IS REJECTED |
---|
454 | REJECT=.TRUE. |
---|
455 | LAST=.FALSE. |
---|
456 | H=HNEW |
---|
457 | IF (NACCPT.GE.1) NREJCT=NREJCT+1 |
---|
458 | END IF |
---|
459 | END DO |
---|
460 | RETURN |
---|
461 | END |
---|
462 | C |
---|
463 | SUBROUTINE FAIL_EXIT(NERR,X,IDID,H,NMAX) |
---|
464 | INTEGER NERR, NMAX |
---|
465 | KPP_REAL X, H |
---|
466 | GO TO (1,2,3,4) NERR |
---|
467 | 1 CONTINUE |
---|
468 | WRITE(6,979)X |
---|
469 | WRITE(6,*) ' MATRIX IS REPEATEDLY SINGULAR, IER=',IER |
---|
470 | IDID=-4 |
---|
471 | STOP |
---|
472 | 2 CONTINUE |
---|
473 | WRITE(6,979)X |
---|
474 | WRITE(6,*) ' STEP SIZE TOO SMALL, H=',H |
---|
475 | IDID=-3 |
---|
476 | STOP |
---|
477 | 3 CONTINUE |
---|
478 | WRITE(6,979)X |
---|
479 | WRITE(6,*) ' MORE THAN NMAX =',NMAX,'STEPS ARE NEEDED' |
---|
480 | IDID=-2 |
---|
481 | STOP |
---|
482 | C --- EXIT CAUSED BY solout |
---|
483 | 4 CONTINUE |
---|
484 | WRITE(6,979)X |
---|
485 | 979 FORMAT(' EXIT OF RODAS AT X=',E18.4) |
---|
486 | IDID=2 |
---|
487 | RETURN |
---|
488 | END |
---|
489 | |
---|
490 | SUBROUTINE ROCOE(METH,A21,A31,A32,A41,A42,A43,A51,A52,A53,A54, |
---|
491 | & C21,C31,C32,C41,C42,C43,C51,C52,C53,C54,C61, |
---|
492 | & C62,C63,C64,C65,GAMMA,C2,C3,C4,D1,D2,D3,D4, |
---|
493 | & D21,D22,D23,D24,D25,D31,D32,D33,D34,D35) |
---|
494 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
495 | |
---|
496 | if (METH.ne.1) print *, 'WRONG CHOICE OF METHOD' |
---|
497 | C2=0.386D0 |
---|
498 | C3=0.21D0 |
---|
499 | C4=0.63D0 |
---|
500 | BET2P=0.0317D0 |
---|
501 | BET3P=0.0635D0 |
---|
502 | BET4P=0.3438D0 |
---|
503 | D1= 0.2500000000000000D+00 |
---|
504 | D2=-0.1043000000000000D+00 |
---|
505 | D3= 0.1035000000000000D+00 |
---|
506 | D4=-0.3620000000000023D-01 |
---|
507 | A21= 0.1544000000000000D+01 |
---|
508 | A31= 0.9466785280815826D+00 |
---|
509 | A32= 0.2557011698983284D+00 |
---|
510 | A41= 0.3314825187068521D+01 |
---|
511 | A42= 0.2896124015972201D+01 |
---|
512 | A43= 0.9986419139977817D+00 |
---|
513 | A51= 0.1221224509226641D+01 |
---|
514 | A52= 0.6019134481288629D+01 |
---|
515 | A53= 0.1253708332932087D+02 |
---|
516 | A54=-0.6878860361058950D+00 |
---|
517 | C21=-0.5668800000000000D+01 |
---|
518 | C31=-0.2430093356833875D+01 |
---|
519 | C32=-0.2063599157091915D+00 |
---|
520 | C41=-0.1073529058151375D+00 |
---|
521 | C42=-0.9594562251023355D+01 |
---|
522 | C43=-0.2047028614809616D+02 |
---|
523 | C51= 0.7496443313967647D+01 |
---|
524 | C52=-0.1024680431464352D+02 |
---|
525 | C53=-0.3399990352819905D+02 |
---|
526 | C54= 0.1170890893206160D+02 |
---|
527 | C61= 0.8083246795921522D+01 |
---|
528 | C62=-0.7981132988064893D+01 |
---|
529 | C63=-0.3152159432874371D+02 |
---|
530 | C64= 0.1631930543123136D+02 |
---|
531 | C65=-0.6058818238834054D+01 |
---|
532 | GAMMA= 0.2500000000000000D+00 |
---|
533 | |
---|
534 | D21= 0.1012623508344586D+02 |
---|
535 | D22=-0.7487995877610167D+01 |
---|
536 | D23=-0.3480091861555747D+02 |
---|
537 | D24=-0.7992771707568823D+01 |
---|
538 | D25= 0.1025137723295662D+01 |
---|
539 | D31=-0.6762803392801253D+00 |
---|
540 | D32= 0.6087714651680015D+01 |
---|
541 | D33= 0.1643084320892478D+02 |
---|
542 | D34= 0.2476722511418386D+02 |
---|
543 | D35=-0.6594389125716872D+01 |
---|
544 | RETURN |
---|
545 | END |
---|
546 | C |
---|
547 | |
---|
548 | C ****************************************** |
---|
549 | C VERSION OF SEPTEMBER 18, 1995 |
---|
550 | C ****************************************** |
---|
551 | C |
---|
552 | SUBROUTINE DECOMR(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
553 | & M1,M2,NM1,FAC1,E1,LDE1,IP1,IER,IJOB,CALHES,IPHES) |
---|
554 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
555 | INCLUDE 'KPP_ROOT_params.h' |
---|
556 | INCLUDE 'KPP_ROOT_sparse.h' |
---|
557 | DIMENSION FJAC(LU_NONZERO),FMAS(LDMAS,NM1),E1(LU_NONZERO), |
---|
558 | & IP1(NM1),IPHES(N) |
---|
559 | LOGICAL CALHES |
---|
560 | COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG |
---|
561 | C |
---|
562 | |
---|
563 | |
---|
564 | C --- B=IDENTITY, JACOBIAN A FULL MATRIX |
---|
565 | DO J=1,LU_NONZERO |
---|
566 | E1(J) = -FJAC(J) |
---|
567 | END DO |
---|
568 | DO J=1,N |
---|
569 | E1(LU_DIAG(J)) = E1(LU_DIAG(J)) + FAC1 |
---|
570 | END DO |
---|
571 | CALL KppDecomp (E1,IER) |
---|
572 | RETURN |
---|
573 | END |
---|
574 | C |
---|
575 | C END OF SUBROUTINE DECOMR |
---|
576 | C |
---|
577 | C *********************************************************** |
---|
578 | C |
---|
579 | C |
---|
580 | C |
---|
581 | C |
---|
582 | SUBROUTINE SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, |
---|
583 | & M1,M2,NM1,FAC1,E,LDE,IP,DY,AK,FX,YNEW,HD,IJOB,STAGE1) |
---|
584 | IMPLICIT KPP_REAL (A-H,O-Z) |
---|
585 | INCLUDE 'KPP_ROOT_params.h' |
---|
586 | INCLUDE 'KPP_ROOT_sparse.h' |
---|
587 | DIMENSION FJAC(LU_NONZERO),FMAS(LDMAS,NM1),E(LU_NONZERO), |
---|
588 | & IP(NM1),DY(N),AK(N),FX(N),YNEW(N) |
---|
589 | LOGICAL STAGE1 |
---|
590 | COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG |
---|
591 | C |
---|
592 | IF (HD.EQ.0.D0) THEN |
---|
593 | DO I=1,N |
---|
594 | AK(I)=DY(I) |
---|
595 | END DO |
---|
596 | ELSE |
---|
597 | DO I=1,N |
---|
598 | AK(I)=DY(I)+HD*FX(I) |
---|
599 | END DO |
---|
600 | END IF |
---|
601 | |
---|
602 | C --- B=IDENTITY, JACOBIAN A FULL MATRIX |
---|
603 | IF (STAGE1) THEN |
---|
604 | DO I=1,N |
---|
605 | AK(I)=AK(I)+YNEW(I) |
---|
606 | END DO |
---|
607 | END IF |
---|
608 | CALL KppSolve (E,AK) |
---|
609 | RETURN |
---|
610 | END |
---|
611 | C |
---|
612 | C END OF SUBROUTINE SLVROD |
---|
613 | C |
---|
614 | C |
---|
615 | C *********************************************************** |
---|
616 | |
---|
617 | |
---|
618 | SUBROUTINE FUNC_CHEM(N, T, Y, P) |
---|
619 | INCLUDE 'KPP_ROOT_params.h' |
---|
620 | INCLUDE 'KPP_ROOT_global.h' |
---|
621 | INTEGER N |
---|
622 | KPP_REAL T, Told |
---|
623 | KPP_REAL Y(NVAR), P(NVAR) |
---|
624 | Told = TIME |
---|
625 | TIME = T |
---|
626 | CALL Update_SUN() |
---|
627 | CALL Update_RCONST() |
---|
628 | CALL Fun( Y, FIX, RCONST, P ) |
---|
629 | TIME = Told |
---|
630 | RETURN |
---|
631 | END |
---|
632 | |
---|
633 | |
---|
634 | SUBROUTINE JAC_CHEM(N, T, Y, J) |
---|
635 | INCLUDE 'KPP_ROOT_params.h' |
---|
636 | INCLUDE 'KPP_ROOT_global.h' |
---|
637 | INTEGER N |
---|
638 | KPP_REAL Told, T |
---|
639 | KPP_REAL Y(NVAR), J(LU_NONZERO) |
---|
640 | Told = TIME |
---|
641 | TIME = T |
---|
642 | CALL Update_SUN() |
---|
643 | CALL Update_RCONST() |
---|
644 | CALL Jac_SP( Y, FIX, RCONST, J ) |
---|
645 | TIME = Told |
---|
646 | RETURN |
---|
647 | END |
---|