[2696] | 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 | |
---|