1 SUBROUTINE INTEGRATE
( TIN
, TOUT
)
3 IMPLICIT KPP_REAL
(A
-H
,O
-Z
)
4 INCLUDE
'KPP_ROOT_params.h'
5 INCLUDE
'KPP_ROOT_global.h'
13 PARAMETER (LWORK
=2*NVAR*NVAR
+14*NVAR
+20,LIWORK
=NVAR
+20)
16 EXTERNAL FUNC_CHEM
,JAC_CHEM
,SOLOUT
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
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
,
37 & WORK
,LWORK
,IWORK
,LIWORK
,IDID
)
40 print
*,'KPP_ROS4: Unsucessful step at T=',
41 & TIN
,' (IDID=',IDID
,')'
48 SUBROUTINE KPP_ROS4
(N
,FCN
,IFCN
,X
,Y
,XEND
,H
,
50 & JAC
,IJAC
,MLJAC
,MUJAC
,DFX
,IDFX
,
51 & MAS
,IMAS
,MLMAS
,MUMAS
,
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).
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
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)
72 C VERSION OF NOVEMBER 17, 1992
76 C N DIMENSION OF THE SYSTEM
78 C FCN NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE
80 C SUBROUTINE FCN(N,X,Y,F)
81 C KPP_REAL X,Y(N),F(N)
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)
90 C Y(N) INITIAL VALUES FOR Y
92 C XEND FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE)
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).
102 C RelTol,AbsTol RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY
103 C CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N.
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).
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)
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
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
130 C DFY(I-J+MUJAC+1,J) = PARTIAL F(I) / PARTIAL Y(J).
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.
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).
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.
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)
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.
162 C ---- MAS,IMAS,MLMAS, AND MUMAS HAVE ANALOG MEANINGS -----
163 C ---- FOR THE "MASS MATRIX" (THE MATRIX "M" OF SECTION IV.8): -
165 C MAS NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE MASS-
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)
174 C IF (MLMAS.EQ.N) THE MASS-MATRIX IS STORED
175 C AS FULL MATRIX LIKE
177 C ELSE, THE MATRIX IS TAKEN AS BANDED AND STORED
179 C AM(I-J+MUMAS+1,J) = M(I,J).
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.
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.
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.
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)
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.
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
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
222 C LJAC=N IF MLJAC=N (FULL JACOBIAN)
223 C LJAC=MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.)
226 C LMAS=N IF IMAS=1 AND MLMAS=N (FULL)
227 C LMAS=MLMAS+MUMAS+1 IF MLMAS<N (BANDED MASS-M.)
229 C LE1=N IF MLJAC=N (FULL JACOBIAN)
230 C LE1=2*MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.).
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.
237 C LWORK DECLARED LENGHT OF ARRAY "WORK".
239 C IWORK INTEGER WORKING SPACE OF LENGTH "LIWORK".
240 C "LIWORK" MUST BE AT LEAST N+2.
242 C LIWORK DECLARED LENGHT OF ARRAY "IWORK".
244 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:
253 C IWORK(1) THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS.
254 C THE DEFAULT VALUE (FOR IWORK(1)=0) IS 100000.
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.
265 C WORK(1) UROUND, THE ROUNDING UNIT, DEFAULT 1.D-16.
267 C WORK(2) MAXIMAL STEP SIZE, DEFAULT XEND-X.
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
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
278 C-----------------------------------------------------------------------
282 C X X-VALUE WHERE THE SOLUTION IS COMPUTED
283 C (AFTER SUCCESSFUL RETURN X=XEND)
287 C H PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP
289 C IDID REPORTS ON SUCCESSFULNESS UPON RETURN:
290 C IDID=1 COMPUTATION SUCCESSFUL,
291 C IDID=-1 COMPUTATION UNSUCCESSFUL.
293 C ---------------------------------------------------------
294 C *** *** *** *** *** *** *** *** *** *** *** *** ***
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
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
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 *** *** *** *** *** *** ***
326 C -------- NMAX , THE MAXIMAL NUMBER OF STEPS -----
327 IF(IWORK
(1).EQ
.0)THEN
332 WRITE(6,*)' WRONG INPUT IWORK(1)=',IWORK
(1)
336 C -------- METH COEFFICIENTS OF THE METHOD
337 IF(IWORK
(2).EQ
.0)THEN
341 IF(METH
.LE
.0.OR
.METH
.GE
.7)THEN
342 WRITE(6,*)' CURIOUS INPUT IWORK(2)=',IWORK
(2)
346 C -------- UROUND SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.D0
347 IF(WORK
(1).EQ
.0.D0
)THEN
351 IF(UROUND
.LE
.1.D
-14.OR
.UROUND
.GE
.1.D0
)THEN
352 WRITE(6,*)' COEFFICIENTS HAVE 16 DIGITS, UROUND=',WORK
(1)
356 C -------- MAXIMAL STEP SIZE
357 IF(WORK
(2).EQ
.0.D0
)THEN
362 C ------- FAC1,FAC2 PARAMETERS FOR STEP SIZE SELECTION
363 IF(WORK
(3).EQ
.0.D0
)THEN
368 IF(WORK
(4).EQ
.0.D0
)THEN
373 C ------- FACREJ FOR THE HUMP
374 IF(WORK
(5).EQ
.0.D0
)THEN
379 C --------- CHECK IF TOLERANCES ARE O.K.
381 IF (AbsTol
(1).LE
.0.D0
.OR
.RelTol
(1).LE
.10.D0*UROUND
) THEN
382 WRITE (6,*) ' TOLERANCES ARE TOO SMALL'
387 IF (AbsTol
(I
).LE
.0.D0
.OR
.RelTol
(I
).LE
.10.D0*UROUND
) THEN
388 WRITE (6,*) ' TOLERANCES(',I
,') ARE TOO SMALL'
393 C *** *** *** *** *** *** *** *** *** *** *** *** ***
394 C COMPUTATION OF ARRAY ENTRIES
395 C *** *** *** *** *** *** *** *** *** *** *** *** ***
396 C ---- AUTONOMOUS, IMPLICIT, BANDED OR NOT ?
401 C -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS ---
408 C -- MATRIX E FOR LINEAR ALGEBRA
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
431 C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
443 C ------ TOTAL STORAGE REQUIREMENT -----------
445 IF(ISTORE
.GT
.LWORK
)THEN
446 WRITE(6,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE
449 C ------- ENTRY POINTS FOR INTEGER WORKSPACE -----
451 C --------- TOTAL REQUIREMENT ---------------
453 IF(ISTORE
.GT
.LIWORK
)THEN
454 WRITE(6,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE
457 C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1
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 -----------
475 C ----- ... AND HERE IS THE CORE INTEGRATOR ----------
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 ----------------------------------------------------------
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)
495 LOGICAL REJECT
,RJECT2
,AUTNMS
,IMPLCT
,BANDED
496 COMMON/STAT
/NFCN
,NJAC
,NSTEP
,NACCPT
,NREJCT
,NDEC
,NSOL
497 EXTERNAL FCN
,JAC
,MAS
,SOLOUT
,DFX
500 C ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ----------
501 IF (IMPLCT
) CALL MAS
(N
,FMAS
,LDMAS
)
502 C ---- PREPARE BANDWIDTHS -----
512 C *** *** *** *** *** *** ***
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
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
546 IF ((X
+H
-XEND
)*POSNEG
.GT
.0.D0
) H
=XEND
-X
550 C *** *** *** *** *** *** ***
551 C COMPUTATION OF THE JACOBIAN
552 C *** *** *** *** *** *** ***
556 IF (.NOT
.AUTNMS
) THEN
558 C --- COMPUTE NUMERICALLY THE DERIVATIVE WITH RESPECT TO X
559 DELT
=DSQRT
(UROUND*MAX
(1.D
-5,ABS
(X
)))
561 CALL FCN
(N
,XDELT
,Y
,AK1
)
563 19 FX
(J
)=(AK1
(J
)-DY1
(J
))/DELT
565 C --- COMPUTE ANALYTICALLY THE DERIVATIVE WITH RESPECT TO X
571 C *** *** *** *** *** *** ***
573 C *** *** *** *** *** *** ***
582 C --- THE MATRIX E (B=IDENTITY, JACOBIAN A FULL MATRIX)
583 DO 800 I
=1,LU_NONZERO
586 801 E
(LU_DIAG
(J
))=E
(LU_DIAG
(J
))+FAC
587 CALL KppDecomp
(E
,IER
)
588 IF (IER
.NE
.0) GOTO 80
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
598 810 YNEW
(I
)=Y
(I
)+A21*AK1
(I
)
599 CALL FCN
(N
,X
,YNEW
,DY
)
601 811 AK2
(I
)=DY
(I
)+HC21*AK1
(I
)
604 820 YNEW
(I
)=Y
(I
)+A31*AK1
(I
)+A32*AK2
(I
)
605 CALL FCN
(N
,X
,YNEW
,DY
)
607 821 AK3
(I
)=DY
(I
)+HC31*AK1
(I
)+HC32*AK2
(I
)
610 831 AK4
(I
)=DY
(I
)+HC41*AK1
(I
)+HC42*AK2
(I
)+HC43*AK3
(I
)
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
622 903 AK1
(I
)=DY1
(I
)+HD1*FX
(I
)
625 910 YNEW
(I
)=Y
(I
)+A21*AK1
(I
)
626 CALL FCN
(N
,X
+C2*H
,YNEW
,DY
)
628 911 AK2
(I
)=DY
(I
)+HD2*FX
(I
)+HC21*AK1
(I
)
631 920 YNEW
(I
)=Y
(I
)+A31*AK1
(I
)+A32*AK2
(I
)
632 CALL FCN
(N
,X
+C3*H
,YNEW
,DY
)
634 921 AK3
(I
)=DY
(I
)+HD3*FX
(I
)+HC31*AK1
(I
)+HC32*AK2
(I
)
637 931 AK4
(I
)=DY
(I
)+HD4*FX
(I
)+HC41*AK1
(I
)+HC42*AK2
(I
)
643 C *** *** *** *** *** *** ***
645 C *** *** *** *** *** *** ***
647 C ------------ NEW SOLUTION ---------------
649 240 YNEW
(I
)=Y
(I
)+B1*AK1
(I
)+B2*AK2
(I
)+B3*AK3
(I
)+B4*AK4
(I
)
650 C ------------ COMPUTE ERROR ESTIMATION ----------------
653 S
=E1*AK1
(I
)+E2*AK2
(I
)+E3*AK3
(I
)+E4*AK4
(I
)
655 SK
=AbsTol
(1)+RelTol
(1)*DMAX1
(DABS
(Y
(I
)),DABS
(YNEW
(I
)))
657 SK
=AbsTol
(I
)+RelTol
(I
)*DMAX1
(DABS
(Y
(I
)),DABS
(YNEW
(I
)))
659 300 ERR
=ERR
+(S
/SK
)**2
661 C --- COMPUTATION OF HNEW
662 C --- WE REQUIRE .2<=HNEW/H<=6.
663 FAC
=DMAX1
(FAC2
,DMIN1
(FAC1
,(ERR
)**.25D0
/.9D0
))
665 C *** *** *** *** *** *** ***
666 C IS THE ERROR SMALL ENOUGH ?
667 C *** *** *** *** *** *** ***
668 IF (ERR
.LE
.1.D0
) THEN
669 C --- STEP IS ACCEPTED
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
))
683 C --- STEP IS REJECTED
684 IF (RJECT2
) HNEW
=H*FACREJ
685 IF (REJECT
) RJECT2
=.TRUE
.
688 IF (NACCPT
.GE
.1) NREJCT
=NREJCT
+1
692 80 WRITE (6,*) ' MATRIX E IS SINGULAR, INFO = ',INFO
694 IF (NSING
.GE
.5) GOTO 79
698 979 FORMAT(' EXIT OF ROS4 AT X=',D16
.7
,' H=',D16
.7
)
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
)
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
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
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
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
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
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
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
879 GAMMA
= 0.5728200000000000D
+00
881 C2
= 0.1145640000000000D
+01
882 C3
= 0.6552168638155900D
+00
884 D1
= 0.5728200000000000D
+00
885 D2
=-0.1769193891319233D
+01
886 D3
= 0.7592633437920482D
+00
887 D4
=-0.1049021087100450D
+00
891 SUBROUTINE SOLOUT
(NR
,XOLD
,X
,Y
,N
,IRTRN
)
892 C --- PRINTS SOLUTION
893 IMPLICIT KPP_REAL
(A
-H
,O
-Z
)
897 WRITE (6,99) X
,Y
(1),Y
(2),NR
-1
901 WRITE (6,99) X
,Y
(1),Y
(2),NR
-1
902 XOUT
=DMAX1
(XOUT
+0.1D0
,X
)
905 99 FORMAT(1X
,'X =',F5
.2
,' Y =',2E18
.10
,' NSTEP =',I4
)
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
913 DIMENSION A
(NDIM
,N
), IP
(N
)
914 C-----------------------------------------------------------------------
915 C MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION.
917 C N = ORDER OF MATRIX.
918 C NDIM = DECLARED DIMENSION OF ARRAY A .
919 C A = MATRIX TO BE TRIANGULARIZED.
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.
932 C C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR,
933 C C.A.C.M. 15 (1972), P. 274.
934 C-----------------------------------------------------------------------
937 IF (N
.EQ
. 1) GO TO 70
943 IF (DABS
(A
(I
,K
)) .GT
. DABS
(A
(M
,K
))) M
= I
947 IF (M
.EQ
. K
) GO TO 20
952 IF (T
.EQ
. 0.D0
) GO TO 80
955 30 A
(I
,K
) = -A
(I
,K
)*T
960 IF (T
.EQ
. 0.D0
) GO TO 45
962 40 A
(I
,J
) = A
(I
,J
) + A
(I
,K
)*T
967 IF (A
(N
,N
) .EQ
. 0.D0
) GO TO 80
972 C----------------------- END OF SUBROUTINE DEC -------------------------
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
980 DIMENSION A
(NDIM
,N
), B
(N
), IP
(N
)
981 C-----------------------------------------------------------------------
982 C SOLUTION OF LINEAR SYSTEM, A*X = B .
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.
991 C B = SOLUTION VECTOR, X .
992 C-----------------------------------------------------------------------
993 IF (N
.EQ
. 1) GO TO 50
1002 10 B
(I
) = B
(I
) + A
(I
,K
)*T
1010 30 B
(I
) = B
(I
) + A
(I
,K
)*T
1012 50 B
(1) = B
(1)/A
(1,1)
1014 C----------------------- END OF SUBROUTINE SOL -------------------------
1020 SUBROUTINE FUNC_CHEM
(N
, T
, Y
, P
)
1021 INCLUDE
'KPP_ROOT_params.h'
1022 INCLUDE
'KPP_ROOT_global.h'
1025 KPP_REAL Y
(NVAR
), P
(NVAR
)
1029 CALL Update_RCONST
()
1030 CALL Fun
( Y
, FIX
, RCONST
, P
)
1036 SUBROUTINE JAC_CHEM
(N
, T
, Y
, J
)
1037 INCLUDE
'KPP_ROOT_params.h'
1038 INCLUDE
'KPP_ROOT_global.h'
1041 KPP_REAL Y
(NVAR
), J
(LU_NONZERO
)
1045 CALL Update_RCONST
()
1046 CALL Jac_SP
( Y
, FIX
, RCONST
, J
)