| 1 | SUBROUTINE INTEGRATE( TIN, TOUT ) |
|---|
| 2 | |
|---|
| 3 | IMPLICIT KPP_REAL (A-H,O-Z) |
|---|
| 4 | INCLUDE 'KPP_ROOT_Parameters.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 (KM=12,KM2=2+KM*(KM+3)/2,NRDENS=NVAR) |
|---|
| 14 | PARAMETER (LWORK=2*NVAR*NVAR+(KM+8)*NVAR+4*KM+20+KM2*NRDENS) |
|---|
| 15 | PARAMETER (LIWORK=2*NVAR+KM+20+NRDENS) |
|---|
| 16 | |
|---|
| 17 | KPP_REAL WORK(LWORK) |
|---|
| 18 | INTEGER IWORK(LIWORK) |
|---|
| 19 | EXTERNAL FUNC_CHEM,JAC_CHEM |
|---|
| 20 | |
|---|
| 21 | ITOL=1 ! --- VECTOR TOLERANCES |
|---|
| 22 | IJAC=1 ! --- COMPUTE THE JACOBIAN ANALYTICALLY |
|---|
| 23 | MLJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX |
|---|
| 24 | MUJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX |
|---|
| 25 | IMAS=0 ! --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM |
|---|
| 26 | IOUT=0 ! --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION |
|---|
| 27 | IDFX=0 ! --- INTERNAL TIME DERIVATIVE |
|---|
| 28 | |
|---|
| 29 | DO i=1,20 |
|---|
| 30 | IWORK(i) = 0 |
|---|
| 31 | WORK(i) = 0.D0 |
|---|
| 32 | ENDDO |
|---|
| 33 | |
|---|
| 34 | CALL ATMSEULEX(NVAR,FUNC_CHEM,Autonomous,TIN,VAR,TOUT, |
|---|
| 35 | & STEPMIN,RTOL,ATOL,ITOL, |
|---|
| 36 | & JAC_CHEM,IJAC,MLJAC,MUJAC, |
|---|
| 37 | & FUNC_CHEM,IMAS,MLJAC,MUJAC, |
|---|
| 38 | & WORK,LWORK,IWORK,LIWORK,IDID) |
|---|
| 39 | |
|---|
| 40 | IF (IDID.LT.0) THEN |
|---|
| 41 | print *,'ATMSEULEX: Unsucessfull exit at T=', |
|---|
| 42 | & TIN,' (IDID=',IDID,')' |
|---|
| 43 | ENDIF |
|---|
| 44 | |
|---|
| 45 | RETURN |
|---|
| 46 | END |
|---|
| 47 | |
|---|
| 48 | |
|---|
| 49 | SUBROUTINE ATMSEULEX(N,FCN,IFCN,X,Y,XEND,H, |
|---|
| 50 | & RelTol,AbsTol,ITOL, |
|---|
| 51 | & JAC ,IJAC,MLJAC,MUJAC, |
|---|
| 52 | & 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 EXTRAPOLATION-ALGORITHM, BASED ON THE |
|---|
| 58 | C LINEARLY IMPLICIT EULER METHOD (WITH STEP SIZE CONTROL |
|---|
| 59 | C AND ORDER SELECTION). |
|---|
| 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@DIVSUN.UNIGE.CH, WANNER@DIVSUN.UNIGE.CH |
|---|
| 65 | C INCLUSION OF DENSE OUTPUT BY E. HAIRER AND A. OSTERMANN |
|---|
| 66 | C |
|---|
| 67 | C THIS CODE IS PART OF THE BOOK: |
|---|
| 68 | C E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL |
|---|
| 69 | C EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. |
|---|
| 70 | C SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS 14, |
|---|
| 71 | C SPRINGER-VERLAG (1991) |
|---|
| 72 | C |
|---|
| 73 | C VERSION OF SEPTEMBER 30, 1995 |
|---|
| 74 | C |
|---|
| 75 | C INPUT PARAMETERS |
|---|
| 76 | C ---------------- |
|---|
| 77 | C N DIMENSION OF THE SYSTEM |
|---|
| 78 | C |
|---|
| 79 | C FCN NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE |
|---|
| 80 | C VALUE OF F(X,Y): |
|---|
| 81 | C SUBROUTINE FCN(N,X,Y,F) |
|---|
| 82 | C KPP_REAL X,Y(N),F(N) |
|---|
| 83 | C F(1)=... ETC. |
|---|
| 84 | C RPAR, IPAR (SEE BELOW) |
|---|
| 85 | C |
|---|
| 86 | C IFCN GIVES INFORMATION ON FCN: |
|---|
| 87 | C IFCN=0: F(X,Y) INDEPENDENT OF X (AUTONOMOUS) |
|---|
| 88 | C IFCN=1: F(X,Y) MAY DEPEND ON X (NON-AUTONOMOUS) |
|---|
| 89 | C |
|---|
| 90 | C X INITIAL X-VALUE |
|---|
| 91 | C |
|---|
| 92 | C Y(N) INITIAL VALUES FOR Y |
|---|
| 93 | C |
|---|
| 94 | C XEND FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE) |
|---|
| 95 | C |
|---|
| 96 | C H INITIAL STEP SIZE GUESS; |
|---|
| 97 | C FOR STIFF EQUATIONS WITH INITIAL TRANSIENT, |
|---|
| 98 | C H=1.D0/(NORM OF F'), USUALLY 1.D-2 OR 1.D-3, IS GOOD. |
|---|
| 99 | C THIS CHOICE IS NOT VERY IMPORTANT, THE CODE QUICKLY |
|---|
| 100 | C ADAPTS ITS STEP SIZE (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 ---- MAS,IMAS,MLMAS, AND MUMAS HAVE ANALOG MEANINGS ----- |
|---|
| 149 | C ---- FOR THE "MASS MATRIX" (THE MATRIX "M" OF SECTION IV.8): - |
|---|
| 150 | C |
|---|
| 151 | C MAS NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE MASS- |
|---|
| 152 | C MATRIX M. |
|---|
| 153 | C IF IMAS=0, THIS MATRIX IS ASSUMED TO BE THE IDENTITY |
|---|
| 154 | C MATRIX AND NEEDS NOT TO BE DEFINED; |
|---|
| 155 | C SUPPLY A DUMMY SUBROUTINE IN THIS CASE. |
|---|
| 156 | C IF IMAS=1, THE SUBROUTINE MAS IS OF THE FORM |
|---|
| 157 | C SUBROUTINE MAS(N,AM,LMAS) |
|---|
| 158 | C KPP_REAL AM(LMAS,N) |
|---|
| 159 | C AM(1,1)= .... |
|---|
| 160 | C IF (MLMAS.EQ.N) THE MASS-MATRIX IS STORED |
|---|
| 161 | C AS FULL MATRIX LIKE |
|---|
| 162 | C AM(I,J) = M(I,J) |
|---|
| 163 | C ELSE, THE MATRIX IS TAKEN AS BANDED AND STORED |
|---|
| 164 | C DIAGONAL-WISE AS |
|---|
| 165 | C AM(I-J+MUMAS+1,J) = M(I,J). |
|---|
| 166 | C |
|---|
| 167 | C IMAS GIVES INFORMATION ON THE MASS-MATRIX: |
|---|
| 168 | C IMAS=0: M IS SUPPOSED TO BE THE IDENTITY |
|---|
| 169 | C MATRIX, MAS IS NEVER CALLED. |
|---|
| 170 | C IMAS=1: MASS-MATRIX IS SUPPLIED. |
|---|
| 171 | C |
|---|
| 172 | C MLMAS SWITCH FOR THE BANDED STRUCTURE OF THE MASS-MATRIX: |
|---|
| 173 | C MLMAS=N: THE FULL MATRIX CASE. THE LINEAR |
|---|
| 174 | C ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION. |
|---|
| 175 | C 0<=MLMAS<N: MLMAS IS THE LOWER BANDWITH OF THE |
|---|
| 176 | C MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW |
|---|
| 177 | C THE MAIN DIAGONAL). |
|---|
| 178 | C MLMAS IS SUPPOSED TO BE .LE. MLJAC. |
|---|
| 179 | C |
|---|
| 180 | C MUMAS UPPER BANDWITH OF MASS-MATRIX (>= NUMBER OF NON- |
|---|
| 181 | C ZERO DIAGONALS ABOVE THE MAIN DIAGONAL). |
|---|
| 182 | C NEED NOT BE DEFINED IF MLMAS=N. |
|---|
| 183 | C MUMAS IS SUPPOSED TO BE .LE. MUJAC. |
|---|
| 184 | C |
|---|
| 185 | C SOLOUT NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE |
|---|
| 186 | C NUMERICAL SOLUTION DURING INTEGRATION. |
|---|
| 187 | C IF IOUT>=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP. |
|---|
| 188 | C SUPPLY A DUMMY SUBROUTINE IF IOUT=0. |
|---|
| 189 | C IT MUST HAVE THE FORM |
|---|
| 190 | C SUBROUTINE SOLOUT (NR,XOLD,X,Y,RC,LRC,IC,LIC,N, |
|---|
| 191 | C RPAR,IPAR,IRTRN) |
|---|
| 192 | C KPP_REAL X,Y(N),RC(LRC),IC(LIC) |
|---|
| 193 | C .... |
|---|
| 194 | C SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH |
|---|
| 195 | C GRID-POINT "X" (THEREBY THE INITIAL VALUE IS |
|---|
| 196 | C THE FIRST GRID-POINT). |
|---|
| 197 | C "XOLD" IS THE PRECEEDING GRID-POINT. |
|---|
| 198 | C "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN |
|---|
| 199 | C IS SET <0, SEULEX RETURNS TO THE CALLING PROGRAM. |
|---|
| 200 | C DO NOT CHANGE THE ENTRIES OF RC(LRC),IC(LIC)! |
|---|
| 201 | C |
|---|
| 202 | C ----- CONTINUOUS OUTPUT (IF IOUT=2): ----- |
|---|
| 203 | C DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION |
|---|
| 204 | C FOR THE INTERVAL [XOLD,X] IS AVAILABLE THROUGH |
|---|
| 205 | C THE KPP_REAL FUNCTION |
|---|
| 206 | C >>> CONTEX(I,S,RC,LRC,IC,LIC) <<< |
|---|
| 207 | C WHICH PROVIDES AN APPROXIMATION TO THE I-TH |
|---|
| 208 | C COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE |
|---|
| 209 | C S SHOULD LIE IN THE INTERVAL [XOLD,X]. |
|---|
| 210 | C |
|---|
| 211 | C IOUT GIVES INFORMATION ON THE SUBROUTINE SOLOUT: |
|---|
| 212 | C IOUT=0: SUBROUTINE IS NEVER CALLED |
|---|
| 213 | C IOUT=1: SUBROUTINE IS USED FOR OUTPUT |
|---|
| 214 | C IOUT=2: DENSE OUTPUT IS PERFORMED IN SOLOUT |
|---|
| 215 | C |
|---|
| 216 | C WORK ARRAY OF WORKING SPACE OF LENGTH "LWORK". |
|---|
| 217 | C SERVES AS WORKING SPACE FOR ALL VECTORS AND MATRICES. |
|---|
| 218 | C "LWORK" MUST BE AT LEAST |
|---|
| 219 | C N*(LJAC+LMAS+LE1+KM+8)+4*KM+20+KM2*NRDENS |
|---|
| 220 | C WHERE |
|---|
| 221 | C KM2=2+KM*(KM+3)/2 AND NRDENS=IWORK(6) (SEE BELOW) |
|---|
| 222 | C AND |
|---|
| 223 | C LJAC=N IF MLJAC=N (FULL JACOBIAN) |
|---|
| 224 | C LJAC=MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.) |
|---|
| 225 | C AND |
|---|
| 226 | C LMAS=0 IF IMAS=0 |
|---|
| 227 | C LMAS=N IF IMAS=1 AND MLMAS=N (FULL) |
|---|
| 228 | C LMAS=MLMAS+MUMAS+1 IF MLMAS<N (BANDED MASS-M.) |
|---|
| 229 | C AND |
|---|
| 230 | C LE1=N IF MLJAC=N (FULL JACOBIAN) |
|---|
| 231 | C LE1=2*MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.). |
|---|
| 232 | C AND |
|---|
| 233 | C KM=12 IF IWORK(3)=0 |
|---|
| 234 | C KM=IWORK(3) IF IWORK(3).GT.0 |
|---|
| 235 | C |
|---|
| 236 | C IN THE USUAL CASE WHERE THE JACOBIAN IS FULL AND THE |
|---|
| 237 | C MASS-MATRIX IS THE INDENTITY (IMAS=0), THE MINIMUM |
|---|
| 238 | C STORAGE REQUIREMENT IS |
|---|
| 239 | C LWORK = 2*N*N+(KM+8)*N+4*KM+13+KM2*NRDENS. |
|---|
| 240 | C IF IWORK(9)=M1>0 THEN "LWORK" MUST BE AT LEAST |
|---|
| 241 | C N*(LJAC+KM+8)+(N-M1)*(LMAS+LE1)+4*KM+20+KM2*NRDENS |
|---|
| 242 | C WHERE IN THE DEFINITIONS OF LJAC, LMAS AND LE1 THE |
|---|
| 243 | C NUMBER N CAN BE REPLACED BY N-M1. |
|---|
| 244 | C |
|---|
| 245 | C LWORK DECLARED LENGTH OF ARRAY "WORK". |
|---|
| 246 | C |
|---|
| 247 | C IWORK INTEGER WORKING SPACE OF LENGTH "LIWORK". |
|---|
| 248 | C "LIWORK" MUST BE AT LEAST 2*N+KM+20+NRDENS. |
|---|
| 249 | C |
|---|
| 250 | C LIWORK DECLARED LENGTH OF ARRAY "IWORK". |
|---|
| 251 | C |
|---|
| 252 | C RPAR, IPAR REAL AND INTEGER PARAMETERS (OR PARAMETER ARRAYS) WHICH |
|---|
| 253 | C CAN BE USED FOR COMMUNICATION BETWEEN YOUR CALLING |
|---|
| 254 | C PROGRAM AND THE FCN, JAC, MAS, SOLOUT SUBROUTINES. |
|---|
| 255 | C |
|---|
| 256 | C ---------------------------------------------------------------------- |
|---|
| 257 | C |
|---|
| 258 | C SOPHISTICATED SETTING OF PARAMETERS |
|---|
| 259 | C ----------------------------------- |
|---|
| 260 | C SEVERAL PARAMETERS OF THE CODE ARE TUNED TO MAKE IT WORK |
|---|
| 261 | C WELL. THEY MAY BE DEFINED BY SETTING WORK(1),..,WORK(13) |
|---|
| 262 | C AS WELL AS IWORK(1),..,IWORK(4) DIFFERENT FROM ZERO. |
|---|
| 263 | C FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES: |
|---|
| 264 | C |
|---|
| 265 | C IWORK(1) IF IWORK(1).NE.0, THE CODE TRANSFORMS THE JACOBIAN |
|---|
| 266 | C MATRIX TO HESSENBERG FORM. THIS IS PARTICULARLY |
|---|
| 267 | C ADVANTAGEOUS FOR LARGE SYSTEMS WITH FULL JACOBIAN. |
|---|
| 268 | C IT DOES NOT WORK FOR BANDED JACOBIAN (MLJAC<N) |
|---|
| 269 | C AND NOT FOR IMPLICIT SYSTEMS (IMAS=1). |
|---|
| 270 | C |
|---|
| 271 | C IWORK(2) THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS. |
|---|
| 272 | C THE DEFAULT VALUE (FOR IWORK(2)=0) IS 100000. |
|---|
| 273 | C |
|---|
| 274 | C IWORK(3) THE MAXIMUM NUMBER OF COLUMNS IN THE EXTRAPOLATION |
|---|
| 275 | C TABLE. THE DEFAULT VALUE (FOR IWORK(3)=0) IS 12. |
|---|
| 276 | C IF IWORK(3).NE.0 THEN IWORK(3) SHOULD BE .GE.3. |
|---|
| 277 | C |
|---|
| 278 | C IWORK(4) SWITCH FOR THE STEP SIZE SEQUENCE |
|---|
| 279 | C IF IWORK(4).EQ.1 THEN 1,2,3,4,6,8,12,16,24,32,48,... |
|---|
| 280 | C IF IWORK(4).EQ.2 THEN 2,3,4,6,8,12,16,24,32,48,64,... |
|---|
| 281 | C IF IWORK(4).EQ.3 THEN 1,2,3,4,5,6,7,8,9,10,... |
|---|
| 282 | C IF IWORK(4).EQ.4 THEN 2,3,4,5,6,7,8,9,10,11,... |
|---|
| 283 | C THE DEFAULT VALUE (FOR IWORK(4)=0) IS IWORK(4)=2. |
|---|
| 284 | C |
|---|
| 285 | C IWORK(5) PARAMETER "LAMBDA" OF DENSE OUTPUT; POSSIBLE VALUES |
|---|
| 286 | C ARE 0 AND 1; DEFAULT IWORK(5)=0. |
|---|
| 287 | C |
|---|
| 288 | C IWORK(6) = NRDENS = NUMBER OF COMPONENTS, FOR WHICH DENSE OUTPUT |
|---|
| 289 | C IS REQUIRED |
|---|
| 290 | C |
|---|
| 291 | C IWORK(21),...,IWORK(NRDENS+20) INDICATE THE COMPONENTS, FOR WHICH |
|---|
| 292 | C DENSE OUTPUT IS REQUIRED |
|---|
| 293 | C |
|---|
| 294 | C IF THE DIFFERENTIAL SYSTEM HAS THE SPECIAL STRUCTURE THAT |
|---|
| 295 | C Y(I)' = Y(I+M2) FOR I=1,...,M1, |
|---|
| 296 | C WITH M1 A MULTIPLE OF M2, A SUBSTANTIAL GAIN IN COMPUTERTIME |
|---|
| 297 | C CAN BE ACHIEVED BY SETTING THE FOLLOWING TWO PARAMETERS. E.G., |
|---|
| 298 | C FOR SECOND ORDER SYSTEMS P'=V, V'=G(P,V), WHERE P AND V ARE |
|---|
| 299 | C VECTORS OF DIMENSION N/2, ONE HAS TO PUT M1=M2=N/2. |
|---|
| 300 | C FOR M1>0 SOME OF THE INPUT PARAMETERS HAVE DIFFERENT MEANINGS: |
|---|
| 301 | C - JAC: ONLY THE ELEMENTS OF THE NON-TRIVIAL PART OF THE |
|---|
| 302 | C JACOBIAN HAVE TO BE STORED |
|---|
| 303 | C IF (MLJAC.EQ.N-M1) THE JACOBIAN IS SUPPOSED TO BE FULL |
|---|
| 304 | C DFY(I,J) = PARTIAL F(I+M1) / PARTIAL Y(J) |
|---|
| 305 | C FOR I=1,N-M1 AND J=1,N. |
|---|
| 306 | C ELSE, THE JACOBIAN IS BANDED ( M1 = M2 * MM ) |
|---|
| 307 | C DFY(I-J+MUJAC+1,J+K*M2) = PARTIAL F(I+M1) / PARTIAL Y(J+K*M2) |
|---|
| 308 | C FOR I=1,MLJAC+MUJAC+1 AND J=1,M2 AND K=0,MM. |
|---|
| 309 | C - MLJAC: MLJAC=N-M1: IF THE NON-TRIVIAL PART OF THE JACOBIAN IS FULL |
|---|
| 310 | C 0<=MLJAC<N-M1: IF THE (MM+1) SUBMATRICES (FOR K=0,MM) |
|---|
| 311 | C PARTIAL F(I+M1) / PARTIAL Y(J+K*M2), I,J=1,M2 |
|---|
| 312 | C ARE BANDED, MLJAC IS THE MAXIMAL LOWER BANDWIDTH |
|---|
| 313 | C OF THESE MM+1 SUBMATRICES |
|---|
| 314 | C - MUJAC: MAXIMAL UPPER BANDWIDTH OF THESE MM+1 SUBMATRICES |
|---|
| 315 | C NEED NOT BE DEFINED IF MLJAC=N-M1 |
|---|
| 316 | C - MAS: IF IMAS=0 THIS MATRIX IS ASSUMED TO BE THE IDENTITY AND |
|---|
| 317 | C NEED NOT BE DEFINED. SUPPLY A DUMMY SUBROUTINE IN THIS CASE. |
|---|
| 318 | C IT IS ASSUMED THAT ONLY THE ELEMENTS OF RIGHT LOWER BLOCK OF |
|---|
| 319 | C DIMENSION N-M1 DIFFER FROM THAT OF THE IDENTITY MATRIX. |
|---|
| 320 | C IF (MLMAS.EQ.N-M1) THIS SUBMATRIX IS SUPPOSED TO BE FULL |
|---|
| 321 | C AM(I,J) = M(I+M1,J+M1) FOR I=1,N-M1 AND J=1,N-M1. |
|---|
| 322 | C ELSE, THE MASS MATRIX IS BANDED |
|---|
| 323 | C AM(I-J+MUMAS+1,J) = M(I+M1,J+M1) |
|---|
| 324 | C - MLMAS: MLMAS=N-M1: IF THE NON-TRIVIAL PART OF M IS FULL |
|---|
| 325 | C 0<=MLMAS<N-M1: LOWER BANDWIDTH OF THE MASS MATRIX |
|---|
| 326 | C - MUMAS: UPPER BANDWIDTH OF THE MASS MATRIX |
|---|
| 327 | C NEED NOT BE DEFINED IF MLMAS=N-M1 |
|---|
| 328 | C |
|---|
| 329 | C IWORK(9) THE VALUE OF M1. DEFAULT M1=0. |
|---|
| 330 | C |
|---|
| 331 | C IWORK(10) THE VALUE OF M2. DEFAULT M2=M1. |
|---|
| 332 | C |
|---|
| 333 | C WORK(1) UROUND, THE ROUNDING UNIT, DEFAULT 1.D-16. |
|---|
| 334 | C |
|---|
| 335 | C WORK(2) MAXIMAL STEP SIZE, DEFAULT XEND-X. |
|---|
| 336 | C |
|---|
| 337 | C WORK(3) DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED; |
|---|
| 338 | C INCREASE WORK(3), TO 0.01 SAY, WHEN JACOBIAN EVALUATIONS |
|---|
| 339 | C ARE COSTLY. FOR SMALL SYSTEMS WORK(3) SHOULD BE SMALLER. |
|---|
| 340 | C DEFAULT MIN(1.0D-4,RelTol(1)) |
|---|
| 341 | C |
|---|
| 342 | C WORK(4), WORK(5) PARAMETERS FOR STEP SIZE SELECTION |
|---|
| 343 | C THE NEW STEP SIZE FOR THE J-TH DIAGONAL ENTRY IS |
|---|
| 344 | C CHOSEN SUBJECT TO THE RESTRICTION |
|---|
| 345 | C FACMIN/WORK(5) <= HNEW(J)/HOLD <= 1/FACMIN |
|---|
| 346 | C WHERE FACMIN=WORK(4)**(1/(J-1)) |
|---|
| 347 | C DEFAULT VALUES: WORK(4)=0.1D0, WORK(5)=4.D0 |
|---|
| 348 | C |
|---|
| 349 | C WORK(6), WORK(7) PARAMETERS FOR THE ORDER SELECTION |
|---|
| 350 | C ORDER IS DECREASED IF W(K-1) <= W(K)*WORK(6) |
|---|
| 351 | C ORDER IS INCREASED IF W(K) <= W(K-1)*WORK(7) |
|---|
| 352 | C DEFAULT VALUES: WORK(6)=0.7D0, WORK(7)=0.9D0 |
|---|
| 353 | C |
|---|
| 354 | C WORK(8), WORK(9) SAFETY FACTORS FOR STEP CONTROL ALGORITHM |
|---|
| 355 | C HNEW=H*WORK(9)*(WORK(8)*TOL/ERR)**(1/(J-1)) |
|---|
| 356 | C DEFAULT VALUES: WORK(8)=0.8D0, WORK(9)=0.93D0 |
|---|
| 357 | C |
|---|
| 358 | C WORK(10), WORK(11), WORK(12), WORK(13) ESTIMATED WORKS FOR |
|---|
| 359 | C A CALL TO FCN, JAC, DEC, SOL, RESPECTIVELY. |
|---|
| 360 | C DEFAULT VALUES ARE: WORK(10)=1.D0, WORK(11)=5.D0, |
|---|
| 361 | C WORK(12)=1.D0, WORK(13)=1.D0. |
|---|
| 362 | C |
|---|
| 363 | C----------------------------------------------------------------------- |
|---|
| 364 | C |
|---|
| 365 | C OUTPUT PARAMETERS |
|---|
| 366 | C ----------------- |
|---|
| 367 | C X X-VALUE WHERE THE SOLUTION IS COMPUTED |
|---|
| 368 | C (AFTER SUCCESSFUL RETURN X=XEND) |
|---|
| 369 | C |
|---|
| 370 | C Y(N) SOLUTION AT X |
|---|
| 371 | C |
|---|
| 372 | C H PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP |
|---|
| 373 | C |
|---|
| 374 | C IDID REPORTS ON SUCCESSFULNESS UPON RETURN: |
|---|
| 375 | C IDID=1 COMPUTATION SUCCESSFUL, |
|---|
| 376 | C IDID=-1 COMPUTATION UNSUCCESSFUL. |
|---|
| 377 | C |
|---|
| 378 | C IWORK(14) NFCN NUMBER OF FUNCTION EVALUATIONS (THOSE FOR NUMERICAL |
|---|
| 379 | C EVALUATION OF THE JACOBIAN ARE NOT COUNTED) |
|---|
| 380 | C IWORK(15) NJAC NUMBER OF JACOBIAN EVALUATIONS (EITHER ANALYTICALLY |
|---|
| 381 | C OR NUMERICALLY) |
|---|
| 382 | C IWORK(16) NSTEP NUMBER OF COMPUTED STEPS |
|---|
| 383 | C IWORK(17) NACCPT NUMBER OF ACCEPTED STEPS |
|---|
| 384 | C IWORK(18) NREJCT NUMBER OF REJECTED STEPS (AFTER AT LEAST ONE STEP |
|---|
| 385 | C HAS BEEN ACCEPTED) |
|---|
| 386 | C IWORK(19) NDEC NUMBER OF LU-DECOMPOSITIONS (N-DIMENSIONAL MATRIX) |
|---|
| 387 | C IWORK(20) NSOL NUMBER OF FORWARD-BACKWARD SUBSTITUTIONS |
|---|
| 388 | C----------------------------------------------------------------------- |
|---|
| 389 | C *** *** *** *** *** *** *** *** *** *** *** *** *** |
|---|
| 390 | C DECLARATIONS |
|---|
| 391 | C *** *** *** *** *** *** *** *** *** *** *** *** *** |
|---|
| 392 | IMPLICIT KPP_REAL (A-H,O-Z) |
|---|
| 393 | DIMENSION Y(N),AbsTol(*),RelTol(*),WORK(LWORK),IWORK(LIWORK) |
|---|
| 394 | LOGICAL AUTNMS,IMPLCT,ARRET,JBAND |
|---|
| 395 | EXTERNAL FCN,JAC,MAS |
|---|
| 396 | C *** *** *** *** *** *** *** |
|---|
| 397 | C SETTING THE PARAMETERS |
|---|
| 398 | C *** *** *** *** *** *** *** |
|---|
| 399 | IOUT = 0 |
|---|
| 400 | NFCN=0 |
|---|
| 401 | NJAC=0 |
|---|
| 402 | NSTEP=0 |
|---|
| 403 | NACCPT=0 |
|---|
| 404 | NREJCT=0 |
|---|
| 405 | NDEC=0 |
|---|
| 406 | NSOL=0 |
|---|
| 407 | ARRET=.FALSE. |
|---|
| 408 | C -------- NMAX , THE MAXIMAL NUMBER OF STEPS ----- |
|---|
| 409 | IF(IWORK(2).EQ.0)THEN |
|---|
| 410 | NMAX=100000 |
|---|
| 411 | ELSE |
|---|
| 412 | NMAX=IWORK(2) |
|---|
| 413 | IF(NMAX.LE.0)THEN |
|---|
| 414 | WRITE(6,*)' WRONG INPUT IWORK(2)=',IWORK(2) |
|---|
| 415 | ARRET=.TRUE. |
|---|
| 416 | END IF |
|---|
| 417 | END IF |
|---|
| 418 | C -------- KM MAXIMUM NUMBER OF COLUMNS IN THE EXTRAPOLATION |
|---|
| 419 | IF(IWORK(3).EQ.0)THEN |
|---|
| 420 | KM=12 |
|---|
| 421 | ELSE |
|---|
| 422 | KM=IWORK(3) |
|---|
| 423 | IF(KM.LE.2)THEN |
|---|
| 424 | WRITE(6,*)' CURIOUS INPUT IWORK(3)=',IWORK(3) |
|---|
| 425 | ARRET=.TRUE. |
|---|
| 426 | END IF |
|---|
| 427 | END IF |
|---|
| 428 | C -------- NSEQU CHOICE OF STEP SIZE SEQUENCE |
|---|
| 429 | NSEQU=IWORK(4) |
|---|
| 430 | IF(IWORK(4).EQ.0) NSEQU=2 |
|---|
| 431 | IF(NSEQU.LE.0.OR.NSEQU.GE.5)THEN |
|---|
| 432 | WRITE(6,*)' CURIOUS INPUT IWORK(4)=',IWORK(4) |
|---|
| 433 | ARRET=.TRUE. |
|---|
| 434 | END IF |
|---|
| 435 | C -------- LAMBDA PARAMETER FOR DENSE OUTPUT |
|---|
| 436 | LAMBDA=IWORK(5) |
|---|
| 437 | IF(LAMBDA.LT.0.OR.LAMBDA.GE.2)THEN |
|---|
| 438 | WRITE(6,*)' CURIOUS INPUT IWORK(5)=',IWORK(5) |
|---|
| 439 | ARRET=.TRUE. |
|---|
| 440 | END IF |
|---|
| 441 | C -------- NRDENS NUMBER OF DENSE OUTPUT COMPONENTS |
|---|
| 442 | NRDENS=IWORK(6) |
|---|
| 443 | IF(NRDENS.LT.0.OR.NRDENS.GT.N)THEN |
|---|
| 444 | WRITE(6,*)' CURIOUS INPUT IWORK(6)=',IWORK(6) |
|---|
| 445 | ARRET=.TRUE. |
|---|
| 446 | END IF |
|---|
| 447 | C -------- PARAMETER FOR SECOND ORDER EQUATIONS |
|---|
| 448 | M1=IWORK(9) |
|---|
| 449 | M2=IWORK(10) |
|---|
| 450 | NM1=N-M1 |
|---|
| 451 | IF (M1.EQ.0) M2=N |
|---|
| 452 | IF (M2.EQ.0) M2=M1 |
|---|
| 453 | IF (M1.LT.0.OR.M2.LT.0.OR.M1+M2.GT.N) THEN |
|---|
| 454 | WRITE(6,*)' CURIOUS INPUT FOR IWORK(9,10)=',M1,M2 |
|---|
| 455 | ARRET=.TRUE. |
|---|
| 456 | END IF |
|---|
| 457 | C -------- UROUND SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.D0 |
|---|
| 458 | IF(WORK(1).EQ.0.D0)THEN |
|---|
| 459 | UROUND=1.D-16 |
|---|
| 460 | ELSE |
|---|
| 461 | UROUND=WORK(1) |
|---|
| 462 | IF(UROUND.LE.0.D0.OR.UROUND.GE.1.D0)THEN |
|---|
| 463 | WRITE(6,*)' UROUND=',WORK(1) |
|---|
| 464 | ARRET=.TRUE. |
|---|
| 465 | END IF |
|---|
| 466 | END IF |
|---|
| 467 | C -------- MAXIMAL STEP SIZE |
|---|
| 468 | IF(WORK(2).EQ.0.D0)THEN |
|---|
| 469 | HMAX=XEND-X |
|---|
| 470 | ELSE |
|---|
| 471 | HMAX=WORK(2) |
|---|
| 472 | END IF |
|---|
| 473 | C ------ THET DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED; |
|---|
| 474 | IF(WORK(3).EQ.0.D0)THEN |
|---|
| 475 | THET=MIN(1.0D-4,RelTol(1)) |
|---|
| 476 | ELSE |
|---|
| 477 | THET=WORK(3) |
|---|
| 478 | END IF |
|---|
| 479 | C ------- FAC1,FAC2 PARAMETERS FOR STEP SIZE SELECTION |
|---|
| 480 | IF(WORK(4).EQ.0.D0)THEN |
|---|
| 481 | FAC1=0.1D0 |
|---|
| 482 | ELSE |
|---|
| 483 | FAC1=WORK(4) |
|---|
| 484 | END IF |
|---|
| 485 | IF(WORK(5).EQ.0.D0)THEN |
|---|
| 486 | FAC2=4.0D0 |
|---|
| 487 | ELSE |
|---|
| 488 | FAC2=WORK(5) |
|---|
| 489 | END IF |
|---|
| 490 | C ------- FAC3, FAC4 PARAMETERS FOR THE ORDER SELECTION |
|---|
| 491 | IF(WORK(6).EQ.0.D0)THEN |
|---|
| 492 | FAC3=0.7D0 |
|---|
| 493 | ELSE |
|---|
| 494 | FAC3=WORK(6) |
|---|
| 495 | END IF |
|---|
| 496 | IF(WORK(7).EQ.0.D0)THEN |
|---|
| 497 | FAC4=0.9D0 |
|---|
| 498 | ELSE |
|---|
| 499 | FAC4=WORK(7) |
|---|
| 500 | END IF |
|---|
| 501 | C ------- SAFE1, SAFE2 SAFETY FACTORS FOR STEP SIZE PREDICTION |
|---|
| 502 | IF(WORK(8).EQ.0.D0)THEN |
|---|
| 503 | SAFE1=0.6D0 |
|---|
| 504 | ELSE |
|---|
| 505 | SAFE1=WORK(8) |
|---|
| 506 | END IF |
|---|
| 507 | IF(WORK(9).EQ.0.D0)THEN |
|---|
| 508 | SAFE2=0.93D0 |
|---|
| 509 | ELSE |
|---|
| 510 | SAFE2=WORK(9) |
|---|
| 511 | END IF |
|---|
| 512 | C ------- WKFCN,WKJAC,WKDEC,WKSOL ESTIMATED WORK FOR FCN,JAC,DEC,SOL |
|---|
| 513 | IF(WORK(10).EQ.0.D0)THEN |
|---|
| 514 | WKFCN=1.D0 |
|---|
| 515 | ELSE |
|---|
| 516 | WKFCN=WORK(10) |
|---|
| 517 | END IF |
|---|
| 518 | IF(WORK(11).EQ.0.D0)THEN |
|---|
| 519 | WKJAC=5.D0 |
|---|
| 520 | ELSE |
|---|
| 521 | WKJAC=WORK(11) |
|---|
| 522 | END IF |
|---|
| 523 | IF(WORK(12).EQ.0.D0)THEN |
|---|
| 524 | WKDEC=1.D0 |
|---|
| 525 | ELSE |
|---|
| 526 | WKDEC=WORK(12) |
|---|
| 527 | END IF |
|---|
| 528 | IF(WORK(13).EQ.0.D0)THEN |
|---|
| 529 | WKSOL=1.D0 |
|---|
| 530 | ELSE |
|---|
| 531 | WKSOL=WORK(13) |
|---|
| 532 | END IF |
|---|
| 533 | WKROW=WKFCN+WKSOL |
|---|
| 534 | C --------- CHECK IF TOLERANCES ARE O.K. |
|---|
| 535 | IF (ITOL.EQ.0) THEN |
|---|
| 536 | IF (AbsTol(1).LE.0.D0.OR.RelTol(1).LE.10.D0*UROUND) THEN |
|---|
| 537 | WRITE (6,*) ' TOLERANCES ARE TOO SMALL' |
|---|
| 538 | ARRET=.TRUE. |
|---|
| 539 | END IF |
|---|
| 540 | ELSE |
|---|
| 541 | DO I=1,N |
|---|
| 542 | IF (AbsTol(I).LE.0.D0.OR.RelTol(I).LE.10.D0*UROUND) THEN |
|---|
| 543 | WRITE (6,*) ' TOLERANCES(',I,') ARE TOO SMALL' |
|---|
| 544 | ARRET=.TRUE. |
|---|
| 545 | END IF |
|---|
| 546 | END DO |
|---|
| 547 | END IF |
|---|
| 548 | C *** *** *** *** *** *** *** *** *** *** *** *** *** |
|---|
| 549 | C COMPUTATION OF ARRAY ENTRIES |
|---|
| 550 | C *** *** *** *** *** *** *** *** *** *** *** *** *** |
|---|
| 551 | C ---- AUTONOMOUS, IMPLICIT, BANDED OR NOT ? |
|---|
| 552 | AUTNMS=IFCN.EQ.0 |
|---|
| 553 | IMPLCT=IMAS.NE.0 |
|---|
| 554 | JBAND=MLJAC.LT.NM1 |
|---|
| 555 | C -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS --- |
|---|
| 556 | C -- JACOBIAN AND MATRIX E |
|---|
| 557 | IF(JBAND)THEN |
|---|
| 558 | LDJAC=MLJAC+MUJAC+1 |
|---|
| 559 | LDE=MLJAC+LDJAC |
|---|
| 560 | ELSE |
|---|
| 561 | MLJAC=NM1 |
|---|
| 562 | MUJAC=NM1 |
|---|
| 563 | LDJAC=NM1 |
|---|
| 564 | LDE=NM1 |
|---|
| 565 | END IF |
|---|
| 566 | C -- MASS MATRIX |
|---|
| 567 | IF (IMPLCT) THEN |
|---|
| 568 | IF (MLMAS.NE.NM1) THEN |
|---|
| 569 | LDMAS=MLMAS+MUMAS+1 |
|---|
| 570 | IF (JBAND) THEN |
|---|
| 571 | IJOB=4 |
|---|
| 572 | ELSE |
|---|
| 573 | IJOB=3 |
|---|
| 574 | END IF |
|---|
| 575 | ELSE |
|---|
| 576 | LDMAS=NM1 |
|---|
| 577 | IJOB=5 |
|---|
| 578 | END IF |
|---|
| 579 | C ------ BANDWITH OF "MAS" NOT LARGER THAN BANDWITH OF "JAC" |
|---|
| 580 | IF (MLMAS.GT.MLJAC.OR.MUMAS.GT.MUJAC) THEN |
|---|
| 581 | WRITE (6,*) 'BANDWITH OF "MAS" NOT LARGER THAN BANDWITH OF |
|---|
| 582 | & "JAC"' |
|---|
| 583 | ARRET=.TRUE. |
|---|
| 584 | END IF |
|---|
| 585 | ELSE |
|---|
| 586 | LDMAS=0 |
|---|
| 587 | IF (JBAND) THEN |
|---|
| 588 | IJOB=2 |
|---|
| 589 | ELSE |
|---|
| 590 | IJOB=1 |
|---|
| 591 | IF (N.GT.2.AND.IWORK(1).NE.0) IJOB=7 |
|---|
| 592 | END IF |
|---|
| 593 | END IF |
|---|
| 594 | LDMAS2=MAX(1,LDMAS) |
|---|
| 595 | C ------ HESSENBERG OPTION ONLY FOR EXPLICIT EQU. WITH FULL JACOBIAN |
|---|
| 596 | IF ((IMPLCT.OR.JBAND).AND.IJOB.EQ.7) THEN |
|---|
| 597 | WRITE(6,*)' HESSENBERG OPTION ONLY FOR EXPLICIT EQUATIONS WITH |
|---|
| 598 | &FULL JACOBIAN' |
|---|
| 599 | ARRET=.TRUE. |
|---|
| 600 | END IF |
|---|
| 601 | C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK ----- |
|---|
| 602 | KM2=(KM*(KM+1))/2 |
|---|
| 603 | IEYH=21 |
|---|
| 604 | IEDY=IEYH+N |
|---|
| 605 | IEFX=IEDY+N |
|---|
| 606 | IEYHH=IEFX+N |
|---|
| 607 | IEDYH=IEYHH+N |
|---|
| 608 | IEDEL=IEDYH+N |
|---|
| 609 | IEWH =IEDEL+N |
|---|
| 610 | IESCAL=IEWH+N |
|---|
| 611 | IEHH =IESCAL+N |
|---|
| 612 | IEW =IEHH+KM |
|---|
| 613 | IEA =IEW+KM |
|---|
| 614 | IEJAC =IEA+KM |
|---|
| 615 | IEE =IEJAC+N*LDJAC |
|---|
| 616 | IEMAS=IEE+NM1*LDE |
|---|
| 617 | IET=IEMAS+NM1*LDMAS |
|---|
| 618 | IFAC=IET+N*KM |
|---|
| 619 | IEDE=IFAC+KM |
|---|
| 620 | IFSAFE=IEDE+(KM+2)*NRDENS |
|---|
| 621 | C ------ TOTAL STORAGE REQUIREMENT ----------- |
|---|
| 622 | ISTORE=IFSAFE+KM2*NRDENS-1 |
|---|
| 623 | IF(ISTORE.GT.LWORK)THEN |
|---|
| 624 | WRITE(6,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE |
|---|
| 625 | ARRET=.TRUE. |
|---|
| 626 | END IF |
|---|
| 627 | C ------- ENTRY POINTS FOR INTEGER WORKSPACE ----- |
|---|
| 628 | IECO=21 |
|---|
| 629 | IEIP=21+NRDENS |
|---|
| 630 | IENJ=IEIP+N |
|---|
| 631 | IEIPH=IENJ+KM |
|---|
| 632 | C --------- TOTAL REQUIREMENT --------------- |
|---|
| 633 | ISTORE=IECO+KM-1 |
|---|
| 634 | IF(ISTORE.GT.LIWORK)THEN |
|---|
| 635 | WRITE(6,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE |
|---|
| 636 | ARRET=.TRUE. |
|---|
| 637 | END IF |
|---|
| 638 | C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1 |
|---|
| 639 | IF (ARRET) THEN |
|---|
| 640 | IDID=-1 |
|---|
| 641 | RETURN |
|---|
| 642 | END IF |
|---|
| 643 | NRD=MAX(1,NRDENS) |
|---|
| 644 | C -------- CALL TO CORE INTEGRATOR ------------ |
|---|
| 645 | CALL SEUCOR(N,FCN,X,Y,XEND,HMAX,H,KM,RelTol,AbsTol,ITOL,JAC,IJAC, |
|---|
| 646 | & MLJAC,MUJAC,MLMAS,MUMAS,IOUT,IDID,IJOB,M1,M2,NM1, |
|---|
| 647 | & NMAX,UROUND,NSEQU,AUTNMS,IMPLCT,JBAND,LDJAC,LDE,LDMAS2, |
|---|
| 648 | & WORK(IEYH),WORK(IEDY),WORK(IEFX),WORK(IEYHH),WORK(IEDYH), |
|---|
| 649 | & WORK(IEDEL),WORK(IEWH),WORK(IESCAL),WORK(IEHH), |
|---|
| 650 | & WORK(IEW),WORK(IEA),WORK(IEJAC),WORK(IEE),WORK(IEMAS), |
|---|
| 651 | & WORK(IET),IWORK(IEIP),IWORK(IENJ),IWORK(IEIPH),FAC1,FAC2,FAC3, |
|---|
| 652 | & FAC4,THET,SAFE1,SAFE2,WKJAC,WKDEC,WKROW,KM2,NRD,WORK(IFAC), |
|---|
| 653 | & WORK(IFSAFE),LAMBDA,NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL, |
|---|
| 654 | & WORK(IEDE),IWORK(IECO)) |
|---|
| 655 | IWORK(14)=NFCN |
|---|
| 656 | IWORK(15)=NJAC |
|---|
| 657 | IWORK(16)=NSTEP |
|---|
| 658 | IWORK(17)=NACCPT |
|---|
| 659 | IWORK(18)=NREJCT |
|---|
| 660 | IWORK(19)=NDEC |
|---|
| 661 | IWORK(20)=NSOL |
|---|
| 662 | C ----------- RETURN ----------- |
|---|
| 663 | RETURN |
|---|
| 664 | END |
|---|
| 665 | C |
|---|
| 666 | C |
|---|
| 667 | C ----- ... AND HERE IS THE CORE INTEGRATOR ---------- |
|---|
| 668 | C |
|---|
| 669 | SUBROUTINE SEUCOR(N,FCN,X,Y,XEND,HMAX,H,KM,RelTol,AbsTol,ITOL, |
|---|
| 670 | & JAC,IJAC,MLJAC,MUJAC,MLB,MUB,IOUT,IDID,IJOB,M1,M2, |
|---|
| 671 | & NM1,NMAX,UROUND,NSEQU,AUTNMS,IMPLCT,BANDED,LFJAC,LE, |
|---|
| 672 | & LDMAS,YH,DY,FX,YHH,DYH,DEL,WH,SCAL,HH,W,A,FJAC,E,FMAS,T,IP, |
|---|
| 673 | & NJ,IPHES,FAC1,FAC2,FAC3,FAC4,THET,SAFE1,SAFE2,WKJAC,WKDEC,WKROW, |
|---|
| 674 | & KM2,NRD,FACUL,FSAFE,LAMBDA,NFCN,NJAC,NSTEP,NACCPT,NREJCT, |
|---|
| 675 | & NDEC,NSOL,DENS,ICOMP) |
|---|
| 676 | C ---------------------------------------------------------- |
|---|
| 677 | C CORE INTEGRATOR FOR SEULEX |
|---|
| 678 | C PARAMETERS SAME AS IN SEULEX WITH WORKSPACE ADDED |
|---|
| 679 | C ---------------------------------------------------------- |
|---|
| 680 | C DECLARATIONS |
|---|
| 681 | C ---------------------------------------------------------- |
|---|
| 682 | IMPLICIT KPP_REAL (A-H,O-Z) |
|---|
| 683 | INCLUDE 'KPP_ROOT_Parameters.h' |
|---|
| 684 | INCLUDE 'KPP_ROOT_Sparse.h' |
|---|
| 685 | INTEGER KM, KM2, K, KC, KRIGHT, KLR, KK, KRN, KOPT |
|---|
| 686 | DIMENSION Y(N),YH(N),DY(N),FX(N),YHH(N),DYH(N),DEL(N) |
|---|
| 687 | DIMENSION WH(N),SCAL(N),HH(KM),W(KM),A(KM),FJAC(LU_NONZERO) |
|---|
| 688 | DIMENSION FMAS(LDMAS,NM1),T(KM,N),IP(NM1),NJ(KM) |
|---|
| 689 | DIMENSION RelTol(*),AbsTol(*) |
|---|
| 690 | DIMENSION IPHES(N),FSAFE(KM2,NRD),FACUL(KM),E(LE,NM1) |
|---|
| 691 | DIMENSION DENS((KM+2)*NRD),ICOMP(NRD) |
|---|
| 692 | LOGICAL REJECT,LAST,ATOV,CALJAC,CALHES,AUTNMS,IMPLCT,BANDED |
|---|
| 693 | EXTERNAL FCN, JAC |
|---|
| 694 | COMMON /COSEU/XOLDD,HHH,NNRD,KRIGHT |
|---|
| 695 | COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG |
|---|
| 696 | C --- COMPUTE COEFFICIENTS FOR DENSE OUTPUT |
|---|
| 697 | IF (IOUT.EQ.2) THEN |
|---|
| 698 | NNRD=NRD |
|---|
| 699 | C --- COMPUTE THE FACTORIALS -------- |
|---|
| 700 | FACUL(1)=1.D0 |
|---|
| 701 | DO I=1,KM-1 |
|---|
| 702 | FACUL(I+1)=I*FACUL(I) |
|---|
| 703 | END DO |
|---|
| 704 | END IF |
|---|
| 705 | |
|---|
| 706 | C *** *** *** *** *** *** *** |
|---|
| 707 | C INITIALISATIONS |
|---|
| 708 | C *** *** *** *** *** *** *** |
|---|
| 709 | LRDE=(KM+2)*NRD |
|---|
| 710 | MLE=MLJAC |
|---|
| 711 | MUE=MUJAC |
|---|
| 712 | MBJAC=MLJAC+MUJAC+1 |
|---|
| 713 | MBB=MLB+MUB+1 |
|---|
| 714 | MDIAG=MLE+MUE+1 |
|---|
| 715 | MDIFF=MLE+MUE-MUB |
|---|
| 716 | MBDIAG=MUB+1 |
|---|
| 717 | IF (M1.GT.0) IJOB=IJOB+10 |
|---|
| 718 | C --- DEFINE THE STEP SIZE SEQUENCE |
|---|
| 719 | IF (NSEQU.EQ.1) THEN |
|---|
| 720 | NJ(1)=1 |
|---|
| 721 | NJ(2)=2 |
|---|
| 722 | NJ(3)=3 |
|---|
| 723 | DO I=4,KM |
|---|
| 724 | NJ(I)=2*NJ(I-2) |
|---|
| 725 | END DO |
|---|
| 726 | END IF |
|---|
| 727 | IF (NSEQU.EQ.2) THEN |
|---|
| 728 | NJ(1)=2 |
|---|
| 729 | NJ(2)=3 |
|---|
| 730 | DO I=3,KM |
|---|
| 731 | NJ(I)=2*NJ(I-2) |
|---|
| 732 | END DO |
|---|
| 733 | END IF |
|---|
| 734 | DO I=1,KM |
|---|
| 735 | IF (NSEQU.EQ.3) NJ(I)=I |
|---|
| 736 | IF (NSEQU.EQ.4) NJ(I)=I+1 |
|---|
| 737 | END DO |
|---|
| 738 | A(1)=WKJAC+NJ(1)*WKROW+WKDEC |
|---|
| 739 | DO I=2,KM |
|---|
| 740 | A(I)=A(I-1)+(NJ(I)-1)*WKROW+WKDEC |
|---|
| 741 | END DO |
|---|
| 742 | K=MAX0(3,MIN0(KM-2,INT(-DLOG10(RelTol(1)+AbsTol(1))*.6D0+1.5D0))) |
|---|
| 743 | HMAXN=MIN(ABS(HMAX),ABS(XEND-X)) |
|---|
| 744 | IF (ABS(H).LE.10.D0*UROUND) H=1.0D-6 |
|---|
| 745 | H=MIN(ABS(H),HMAXN) |
|---|
| 746 | THETA=2*ABS(THET) |
|---|
| 747 | ERR=0.D0 |
|---|
| 748 | W(1)=1.D30 |
|---|
| 749 | DO I=1,N |
|---|
| 750 | IF (ITOL.EQ.0) THEN |
|---|
| 751 | SCAL(I)=AbsTol(1)+RelTol(1)*DABS(Y(I)) |
|---|
| 752 | ELSE |
|---|
| 753 | SCAL(I)=AbsTol(I)+RelTol(I)*DABS(Y(I)) |
|---|
| 754 | END IF |
|---|
| 755 | END DO |
|---|
| 756 | CALJAC=.FALSE. |
|---|
| 757 | REJECT=.FALSE. |
|---|
| 758 | LAST=.FALSE. |
|---|
| 759 | 10 CONTINUE |
|---|
| 760 | IF (REJECT) THETA=2*ABS(THET) |
|---|
| 761 | ATOV=.FALSE. |
|---|
| 762 | C *** *** *** *** *** *** *** |
|---|
| 763 | C --- IS XEND REACHED IN THE NEXT STEP? |
|---|
| 764 | C *** *** *** *** *** *** *** |
|---|
| 765 | H1=XEND-X |
|---|
| 766 | IF (H1.LE.UROUND) GO TO 110 |
|---|
| 767 | HOPT=H |
|---|
| 768 | H=MIN(H,H1,HMAXN) |
|---|
| 769 | IF (H.GE.H1-UROUND) LAST=.TRUE. |
|---|
| 770 | IF (AUTNMS) THEN |
|---|
| 771 | CALL FCN(N,X,Y,DY) |
|---|
| 772 | NFCN=NFCN+1 |
|---|
| 773 | END IF |
|---|
| 774 | IF (THETA.GT.THET.AND..NOT.CALJAC) THEN |
|---|
| 775 | C *** *** *** *** *** *** *** |
|---|
| 776 | C COMPUTATION OF THE JACOBIAN |
|---|
| 777 | C *** *** *** *** *** *** *** |
|---|
| 778 | NJAC=NJAC+1 |
|---|
| 779 | C --- COMPUTE JACOBIAN MATRIX ANALYTICALLY |
|---|
| 780 | CALL JAC(N,X,Y,FJAC) |
|---|
| 781 | CALJAC=.TRUE. |
|---|
| 782 | CALHES=.FALSE. |
|---|
| 783 | END IF |
|---|
| 784 | C *** *** *** *** *** *** *** |
|---|
| 785 | C --- THE FIRST AND LAST STEP |
|---|
| 786 | C *** *** *** *** *** *** *** |
|---|
| 787 | IF (NSTEP.EQ.0.OR.LAST) THEN |
|---|
| 788 | IPT=0 |
|---|
| 789 | NSTEP=NSTEP+1 |
|---|
| 790 | DO J=1,K |
|---|
| 791 | KC=J |
|---|
| 792 | CALL SEUL(J,N,FCN,X,Y,DY,FX,FJAC,LFJAC,FMAS,LDMAS,E,LE,IP,H,KM, |
|---|
| 793 | & HMAXN,T,SCAL,NJ,HH,W,A,YHH,DYH,DEL,WH,ERR,SAFE1,FAC, |
|---|
| 794 | & FAC1,FAC2,SAFE2,THETA,MLJAC,MUJAC,NFCN,NDEC,NSOL,MLB, |
|---|
| 795 | & MUB,ERROLD,IPHES,ICOMP,AUTNMS,IMPLCT,BANDED,REJECT, |
|---|
| 796 | & ATOV,FSAFE,KM2,NRD,IOUT,IPT,M1,M2,NM1,IJOB,CALHES) |
|---|
| 797 | IF (ATOV) GOTO 10 |
|---|
| 798 | IF (J.GT.1.AND.ERR.LE.1.D0) GO TO 60 |
|---|
| 799 | END DO |
|---|
| 800 | GO TO 55 |
|---|
| 801 | END IF |
|---|
| 802 | C --- BASIC INTEGRATION STEP |
|---|
| 803 | 30 CONTINUE |
|---|
| 804 | IPT=0 |
|---|
| 805 | NSTEP=NSTEP+1 |
|---|
| 806 | IF (NSTEP.GE.NMAX) GO TO 120 |
|---|
| 807 | KC=K-1 |
|---|
| 808 | DO J=1,KC |
|---|
| 809 | CALL SEUL(J,N,FCN,X,Y,DY,FX,FJAC,LFJAC,FMAS,LDMAS,E,LE,IP,H,KM, |
|---|
| 810 | & HMAXN,T,SCAL,NJ,HH,W,A,YHH,DYH,DEL,WH,ERR,SAFE1, |
|---|
| 811 | & FAC,FAC1,FAC2,SAFE2,THETA,MLJAC,MUJAC,NFCN,NDEC,NSOL, |
|---|
| 812 | & MLB,MUB,ERROLD,IPHES,ICOMP,AUTNMS,IMPLCT,BANDED,REJECT, |
|---|
| 813 | & ATOV,FSAFE,KM2,NRD,IOUT,IPT,M1,M2,NM1,IJOB,CALHES) |
|---|
| 814 | IF (ATOV) GOTO 10 |
|---|
| 815 | END DO |
|---|
| 816 | C *** *** *** *** *** *** *** |
|---|
| 817 | C --- CONVERGENCE MONITOR |
|---|
| 818 | C *** *** *** *** *** *** *** |
|---|
| 819 | IF (K.EQ.2.OR.REJECT) GO TO 50 |
|---|
| 820 | IF (ERR.LE.1.D0) GO TO 60 |
|---|
| 821 | IF (ERR.GT.DBLE(NJ(K+1)*NJ(K))*4.D0) GO TO 100 |
|---|
| 822 | 50 CALL SEUL(K,N,FCN,X,Y,DY,FX,FJAC,LFJAC,FMAS,LDMAS,E,LE,IP,H,KM, |
|---|
| 823 | & HMAXN,T,SCAL,NJ,HH,W,A,YHH,DYH,DEL,WH,ERR,SAFE1, |
|---|
| 824 | & FAC,FAC1,FAC2,SAFE2,THETA,MLJAC,MUJAC,NFCN,NDEC,NSOL, |
|---|
| 825 | & MLB,MUB,ERROLD,IPHES,ICOMP,AUTNMS,IMPLCT,BANDED,REJECT, |
|---|
| 826 | & ATOV,FSAFE,KM2,NRD,IOUT,IPT,M1,M2,NM1,IJOB,CALHES) |
|---|
| 827 | IF (ATOV) GOTO 10 |
|---|
| 828 | KC=K |
|---|
| 829 | IF (ERR.LE.1.D0) GO TO 60 |
|---|
| 830 | C --- HOPE FOR CONVERGENCE IN LINE K+1 |
|---|
| 831 | 55 IF (ERR.GT.DBLE(NJ(K+1))*2.D0) GO TO 100 |
|---|
| 832 | KC=K+1 |
|---|
| 833 | CALL SEUL(KC,N,FCN,X,Y,DY,FX,FJAC,LFJAC,FMAS,LDMAS,E,LE,IP,H,KM, |
|---|
| 834 | & HMAXN,T,SCAL,NJ,HH,W,A,YHH,DYH,DEL,WH,ERR,SAFE1, |
|---|
| 835 | & FAC,FAC1,FAC2,SAFE2,THETA,MLJAC,MUJAC,NFCN,NDEC,NSOL, |
|---|
| 836 | & MLB,MUB,ERROLD,IPHES,ICOMP,AUTNMS,IMPLCT,BANDED,REJECT, |
|---|
| 837 | & ATOV,FSAFE,KM2,NRD,IOUT,IPT,M1,M2,NM1,IJOB,CALHES) |
|---|
| 838 | IF (ATOV) GOTO 10 |
|---|
| 839 | CAdi IF (ERR.GT.1.D0) GO TO 100 |
|---|
| 840 | IF ((ERR.GT.1.D0).and.(H.gt.STEPMIN)) GO TO 100 ! Adi |
|---|
| 841 | C *** *** *** *** *** *** *** |
|---|
| 842 | C --- STEP IS ACCEPTED |
|---|
| 843 | C *** *** *** *** *** *** *** |
|---|
| 844 | 60 XOLD=X |
|---|
| 845 | X=X+H |
|---|
| 846 | IF (IOUT.EQ.2) THEN |
|---|
| 847 | KRIGHT=KC |
|---|
| 848 | DO I=1,NRD |
|---|
| 849 | DENS(I)=Y(ICOMP(I)) |
|---|
| 850 | END DO |
|---|
| 851 | END IF |
|---|
| 852 | DO I=1,N |
|---|
| 853 | T1I=T(1,I) |
|---|
| 854 | IF (ITOL.EQ.0) THEN |
|---|
| 855 | SCAL(I)=AbsTol(1)+RelTol(1)*DABS(T1I) |
|---|
| 856 | ELSE |
|---|
| 857 | SCAL(I)=AbsTol(I)+RelTol(I)*DABS(T1I) |
|---|
| 858 | END IF |
|---|
| 859 | Y(I)=T1I |
|---|
| 860 | END DO |
|---|
| 861 | NACCPT=NACCPT+1 |
|---|
| 862 | CALJAC=.FALSE. |
|---|
| 863 | IF (IOUT.EQ.2) THEN |
|---|
| 864 | XOLDD=XOLD |
|---|
| 865 | HHH=H |
|---|
| 866 | DO I=1,NRD |
|---|
| 867 | DENS(NRD+I)=Y(ICOMP(I)) |
|---|
| 868 | END DO |
|---|
| 869 | DO KLR=1,KRIGHT-1 |
|---|
| 870 | C --- COMPUTE DIFFERENCES |
|---|
| 871 | IF (KLR.GE.2) THEN |
|---|
| 872 | DO KK=KLR,KC |
|---|
| 873 | LBEG=((KK+1)*KK)/2 |
|---|
| 874 | LEND=LBEG-KK+2 |
|---|
| 875 | DO L=LBEG,LEND,-1 |
|---|
| 876 | DO I=1,NRD |
|---|
| 877 | FSAFE(L,I)=FSAFE(L,I)-FSAFE(L-1,I) |
|---|
| 878 | END DO |
|---|
| 879 | END DO |
|---|
| 880 | END DO |
|---|
| 881 | END IF |
|---|
| 882 | C --- COMPUTE DERIVATIVES AT RIGHT END ---- |
|---|
| 883 | DO KK=KLR+LAMBDA,KC |
|---|
| 884 | FACNJ=NJ(KK) |
|---|
| 885 | FACNJ=FACNJ**KLR/FACUL(KLR+1) |
|---|
| 886 | IPT=((KK+1)*KK)/2 |
|---|
| 887 | DO I=1,NRD |
|---|
| 888 | KRN=(KK-LAMBDA+1)*NRD |
|---|
| 889 | DENS(KRN+I)=FSAFE(IPT,I)*FACNJ |
|---|
| 890 | END DO |
|---|
| 891 | END DO |
|---|
| 892 | DO J=KLR+LAMBDA+1,KC |
|---|
| 893 | DBLENJ=NJ(J) |
|---|
| 894 | DO L=J,KLR+LAMBDA+1,-1 |
|---|
| 895 | FACTOR=DBLENJ/NJ(L-1)-1.D0 |
|---|
| 896 | DO I=1,NRD |
|---|
| 897 | KRN=(L-LAMBDA+1)*NRD+I |
|---|
| 898 | DENS(KRN-NRD)=DENS(KRN) |
|---|
| 899 | & +(DENS(KRN)-DENS(KRN-NRD))/FACTOR |
|---|
| 900 | END DO |
|---|
| 901 | END DO |
|---|
| 902 | END DO |
|---|
| 903 | END DO |
|---|
| 904 | C --- COMPUTE THE COEFFICIENTS OF THE INTERPOLATION POLYNOMIAL |
|---|
| 905 | DO IN=1,NRD |
|---|
| 906 | DO J=1,KRIGHT |
|---|
| 907 | II=NRD*J+IN |
|---|
| 908 | DENS(II)=DENS(II)-DENS(II-NRD) |
|---|
| 909 | END DO |
|---|
| 910 | END DO |
|---|
| 911 | END IF |
|---|
| 912 | C --- COMPUTE OPTIMAL ORDER |
|---|
| 913 | IF (KC.EQ.2) THEN |
|---|
| 914 | KOPT=3 |
|---|
| 915 | IF (REJECT) KOPT=2 |
|---|
| 916 | GO TO 80 |
|---|
| 917 | END IF |
|---|
| 918 | IF (KC.LE.K) THEN |
|---|
| 919 | KOPT=KC |
|---|
| 920 | IF (W(KC-1).LT.W(KC)*FAC3) KOPT=KC-1 |
|---|
| 921 | IF (W(KC).LT.W(KC-1)*FAC4) KOPT=MIN0(KC+1,KM-1) |
|---|
| 922 | ELSE |
|---|
| 923 | KOPT=KC-1 |
|---|
| 924 | IF (KC.GT.3.AND.W(KC-2).LT.W(KC-1)*FAC3) KOPT=KC-2 |
|---|
| 925 | IF (W(KC).LT.W(KOPT)*FAC4) KOPT=MIN0(KC,KM-1) |
|---|
| 926 | END IF |
|---|
| 927 | C --- AFTER A REJECTED STEP |
|---|
| 928 | 80 IF (REJECT) THEN |
|---|
| 929 | K=MIN0(KOPT,KC) |
|---|
| 930 | H=DMIN1(H,HH(K)) |
|---|
| 931 | REJECT=.FALSE. |
|---|
| 932 | GO TO 10 |
|---|
| 933 | END IF |
|---|
| 934 | C --- COMPUTE STEP SIZE FOR NEXT STEP |
|---|
| 935 | IF (KOPT.LE.KC) THEN |
|---|
| 936 | H=HH(KOPT) |
|---|
| 937 | ELSE |
|---|
| 938 | IF (KC.LT.K.AND.W(KC).LT.W(KC-1)*FAC4) THEN |
|---|
| 939 | H=HH(KC)*A(KOPT+1)/A(KC) |
|---|
| 940 | ELSE |
|---|
| 941 | H=HH(KC)*A(KOPT)/A(KC) |
|---|
| 942 | END IF |
|---|
| 943 | END IF |
|---|
| 944 | K=KOPT |
|---|
| 945 | H = DMAX1(H, STEPMIN) ! Adi |
|---|
| 946 | GO TO 10 |
|---|
| 947 | C *** *** *** *** *** *** *** |
|---|
| 948 | C --- STEP IS REJECTED |
|---|
| 949 | C *** *** *** *** *** *** *** |
|---|
| 950 | 100 K=MIN0(K,KC) |
|---|
| 951 | IF (K.GT.2.AND.W(K-1).LT.W(K)*FAC3) K=K-1 |
|---|
| 952 | NREJCT=NREJCT+1 |
|---|
| 953 | H=HH(K) |
|---|
| 954 | LAST=.FALSE. |
|---|
| 955 | REJECT=.TRUE. |
|---|
| 956 | IF (CALJAC) GOTO 30 |
|---|
| 957 | GO TO 10 |
|---|
| 958 | C --- SOLUTION EXIT |
|---|
| 959 | 110 CONTINUE |
|---|
| 960 | H=HOPT |
|---|
| 961 | IDID=1 |
|---|
| 962 | RETURN |
|---|
| 963 | C --- FAIL EXIT |
|---|
| 964 | 120 WRITE (6,979) X,H |
|---|
| 965 | 979 FORMAT(' EXIT OF SEULEX AT X=',D14.7,' H=',D14.7) |
|---|
| 966 | IDID=-1 |
|---|
| 967 | RETURN |
|---|
| 968 | END |
|---|
| 969 | C |
|---|
| 970 | C |
|---|
| 971 | C *** *** *** *** *** *** *** |
|---|
| 972 | C S U B R O U T I N E S E U L |
|---|
| 973 | C *** *** *** *** *** *** *** |
|---|
| 974 | C |
|---|
| 975 | SUBROUTINE SEUL(JJ,N,FCN,X,Y,DY,FX,FJAC,LFJAC,FMAS,LDMAS,E,LE,IP, |
|---|
| 976 | & H,KM,HMAXN,T,SCAL,NJ,HH,W,A,YH,DYH,DEL,WH,ERR,SAFE1, |
|---|
| 977 | & FAC,FAC1,FAC2,SAFE2,THETA,MLJAC,MUJAC,NFCN,NDEC,NSOL, |
|---|
| 978 | & MLB,MUB,ERROLD,IPHES,ICOMP, |
|---|
| 979 | & AUTNMS,IMPLCT,BANDED,REJECT,ATOV,FSAFE,KM2,NRD,IOUT, |
|---|
| 980 | & IPT,M1,M2,NM1,IJOB,CALHES) |
|---|
| 981 | C --- THIS SUBROUTINE COMPUTES THE J-TH LINE OF THE |
|---|
| 982 | C --- EXTRAPOLATION TABLE AND PROVIDES AN ESTIMATE |
|---|
| 983 | C --- OF THE OPTIMAL STEP SIZE |
|---|
| 984 | IMPLICIT KPP_REAL (A-H,O-Z) |
|---|
| 985 | INCLUDE 'KPP_ROOT_Parameters.h' |
|---|
| 986 | INCLUDE 'KPP_ROOT_Sparse.h' |
|---|
| 987 | INTEGER KM, KM2 |
|---|
| 988 | DIMENSION Y(N),YH(N),DY(N),FX(N),DYH(N),DEL(N) |
|---|
| 989 | DIMENSION WH(N),SCAL(N),HH(KM),W(KM),A(KM) |
|---|
| 990 | DIMENSION FJAC(LU_NONZERO),E(LU_NONZERO) |
|---|
| 991 | DIMENSION FMAS(LDMAS,N),T(KM,N),IP(N),NJ(KM),IPHES(N) |
|---|
| 992 | DIMENSION FSAFE(KM2,NRD),ICOMP(NRD) |
|---|
| 993 | LOGICAL ATOV,REJECT,AUTNMS,IMPLCT,BANDED,CALHES |
|---|
| 994 | COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG |
|---|
| 995 | EXTERNAL FCN |
|---|
| 996 | C *** *** *** *** *** *** *** |
|---|
| 997 | C COMPUTE THE MATRIX E AND ITS DECOMPOSITION |
|---|
| 998 | C *** *** *** *** *** *** *** |
|---|
| 999 | HJ=H/NJ(JJ) |
|---|
| 1000 | HJI=1.D0/HJ |
|---|
| 1001 | DO I=1,LU_NONZERO |
|---|
| 1002 | E(I)=-FJAC(I) |
|---|
| 1003 | END DO |
|---|
| 1004 | DO J=1,N |
|---|
| 1005 | E(LU_DIAG(J))=E(LU_DIAG(J))+HJI |
|---|
| 1006 | END DO |
|---|
| 1007 | CALL KppDecomp (E,IER) |
|---|
| 1008 | IF (IER.NE.0) GOTO 79 |
|---|
| 1009 | NDEC=NDEC+1 |
|---|
| 1010 | C *** *** *** *** *** *** *** |
|---|
| 1011 | C --- STARTING PROCEDURE |
|---|
| 1012 | C *** *** *** *** *** *** *** |
|---|
| 1013 | IF (.NOT.AUTNMS) THEN |
|---|
| 1014 | CALL FCN(N,X+HJ,Y,DY) |
|---|
| 1015 | NFCN=NFCN+1 |
|---|
| 1016 | END IF |
|---|
| 1017 | DO I=1,N |
|---|
| 1018 | YH(I)=Y(I) |
|---|
| 1019 | DEL(I)=DY(I) |
|---|
| 1020 | END DO |
|---|
| 1021 | CALL KppSolve (E,DEL) |
|---|
| 1022 | NSOL=NSOL+1 |
|---|
| 1023 | M=NJ(JJ) |
|---|
| 1024 | IF (IOUT.EQ.2.AND.M.EQ.JJ) THEN |
|---|
| 1025 | IPT=IPT+1 |
|---|
| 1026 | DO I=1,NRD |
|---|
| 1027 | FSAFE(IPT,I)=DEL(ICOMP(I)) |
|---|
| 1028 | END DO |
|---|
| 1029 | END IF |
|---|
| 1030 | C *** *** *** *** *** *** *** |
|---|
| 1031 | C --- SEMI-IMPLICIT EULER METHOD |
|---|
| 1032 | C *** *** *** *** *** *** *** |
|---|
| 1033 | IF (M.GT.1) THEN |
|---|
| 1034 | DO MM=1,M-1 |
|---|
| 1035 | DO I=1,N |
|---|
| 1036 | YH(I)=YH(I)+DEL(I) |
|---|
| 1037 | END DO |
|---|
| 1038 | IF (AUTNMS) THEN |
|---|
| 1039 | CALL FCN(N,X+HJ*MM,YH,DYH) |
|---|
| 1040 | ELSE |
|---|
| 1041 | CALL FCN(N,X+HJ*(MM+1),YH,DYH) |
|---|
| 1042 | END IF |
|---|
| 1043 | NFCN=NFCN+1 |
|---|
| 1044 | IF (MM.EQ.1.AND.JJ.LE.2) THEN |
|---|
| 1045 | C --- STABILITY CHECK |
|---|
| 1046 | DEL1=0.D0 |
|---|
| 1047 | DO I=1,N |
|---|
| 1048 | DEL1=DEL1+(DEL(I)/SCAL(I))**2 |
|---|
| 1049 | END DO |
|---|
| 1050 | DEL1=DSQRT(DEL1) |
|---|
| 1051 | IF (IMPLCT) THEN |
|---|
| 1052 | DO I=1,NM1 |
|---|
| 1053 | WH(I)=DEL(I+M1) |
|---|
| 1054 | END DO |
|---|
| 1055 | IF (MLB.EQ.NM1) THEN |
|---|
| 1056 | DO I=1,NM1 |
|---|
| 1057 | SUM=0.D0 |
|---|
| 1058 | DO J=1,NM1 |
|---|
| 1059 | SUM=SUM+FMAS(I,J)*WH(J) |
|---|
| 1060 | END DO |
|---|
| 1061 | DEL(I+M1)=SUM |
|---|
| 1062 | END DO |
|---|
| 1063 | ELSE |
|---|
| 1064 | DO I=1,NM1 |
|---|
| 1065 | SUM=0.D0 |
|---|
| 1066 | DO J=MAX(1,I-MLB),MIN(NM1,I+MUB) |
|---|
| 1067 | SUM=SUM+FMAS(I-J+MBDIAG,J)*WH(J) |
|---|
| 1068 | END DO |
|---|
| 1069 | DEL(I+M1)=SUM |
|---|
| 1070 | END DO |
|---|
| 1071 | END IF |
|---|
| 1072 | END IF |
|---|
| 1073 | IF (.NOT.AUTNMS) THEN |
|---|
| 1074 | CALL FCN(N,X+HJ,YH,WH) |
|---|
| 1075 | NFCN=NFCN+1 |
|---|
| 1076 | DO I=1,N |
|---|
| 1077 | DEL(I)=WH(I)-DEL(I)*HJI |
|---|
| 1078 | END DO |
|---|
| 1079 | ELSE |
|---|
| 1080 | DO I=1,N |
|---|
| 1081 | DEL(I)=DYH(I)-DEL(I)*HJI |
|---|
| 1082 | END DO |
|---|
| 1083 | END IF |
|---|
| 1084 | CALL KppSolve (E,DEL) |
|---|
| 1085 | NSOL=NSOL+1 |
|---|
| 1086 | DEL2=0.D0 |
|---|
| 1087 | DO I=1,N |
|---|
| 1088 | DEL2=DEL2+(DEL(I)/SCAL(I))**2 |
|---|
| 1089 | END DO |
|---|
| 1090 | DEL2=DSQRT(DEL2) |
|---|
| 1091 | THETA=DEL2/MAX(1.D0,DEL1) |
|---|
| 1092 | IF (THETA.GT.1.D0) GOTO 79 |
|---|
| 1093 | END IF |
|---|
| 1094 | CALL KppSolve (E,DYH) |
|---|
| 1095 | NSOL=NSOL+1 |
|---|
| 1096 | DO I=1,N |
|---|
| 1097 | DEL(I)=DYH(I) |
|---|
| 1098 | END DO |
|---|
| 1099 | IF (IOUT.EQ.2.AND.MM.GE.M-JJ) THEN |
|---|
| 1100 | IPT=IPT+1 |
|---|
| 1101 | DO I=1,NRD |
|---|
| 1102 | FSAFE(IPT,I)=DEL(ICOMP(I)) |
|---|
| 1103 | END DO |
|---|
| 1104 | END IF |
|---|
| 1105 | END DO |
|---|
| 1106 | END IF |
|---|
| 1107 | DO I=1,N |
|---|
| 1108 | T(JJ,I)=YH(I)+DEL(I) |
|---|
| 1109 | END DO |
|---|
| 1110 | C *** *** *** *** *** *** *** |
|---|
| 1111 | C --- POLYNOMIAL EXTRAPOLATION |
|---|
| 1112 | C *** *** *** *** *** *** *** |
|---|
| 1113 | IF (JJ.EQ.1) RETURN |
|---|
| 1114 | DO L=JJ,2,-1 |
|---|
| 1115 | FAC=(DBLE(NJ(JJ))/DBLE(NJ(L-1)))-1.D0 |
|---|
| 1116 | DO I=1,N |
|---|
| 1117 | T(L-1,I)=T(L,I)+(T(L,I)-T(L-1,I))/FAC |
|---|
| 1118 | END DO |
|---|
| 1119 | END DO |
|---|
| 1120 | ERR=0.D0 |
|---|
| 1121 | DO I=1,N |
|---|
| 1122 | ERR=ERR+MIN(ABS((T(1,I)-T(2,I)))/SCAL(I),1.D15)**2 |
|---|
| 1123 | END DO |
|---|
| 1124 | IF (ERR.GE.1.D30) GOTO 79 |
|---|
| 1125 | ERR=DSQRT(ERR/DBLE(N)) |
|---|
| 1126 | IF (JJ.GT.2.AND.ERR.GE.ERROLD) GOTO 79 |
|---|
| 1127 | ERROLD=DMAX1(4*ERR,1.D0) |
|---|
| 1128 | C --- COMPUTE OPTIMAL STEP SIZES |
|---|
| 1129 | EXPO=1.D0/JJ |
|---|
| 1130 | FACMIN=FAC1**EXPO |
|---|
| 1131 | FAC=DMIN1(FAC2/FACMIN,DMAX1(FACMIN,(ERR/SAFE1)**EXPO/SAFE2)) |
|---|
| 1132 | FAC=1.D0/FAC |
|---|
| 1133 | HH(JJ)=DMIN1(H*FAC,HMAXN) |
|---|
| 1134 | W(JJ)=A(JJ)/HH(JJ) |
|---|
| 1135 | RETURN |
|---|
| 1136 | 79 ATOV=.TRUE. |
|---|
| 1137 | H=H*0.5D0 |
|---|
| 1138 | REJECT=.TRUE. |
|---|
| 1139 | RETURN |
|---|
| 1140 | END |
|---|
| 1141 | |
|---|
| 1142 | |
|---|
| 1143 | |
|---|
| 1144 | SUBROUTINE FUNC_CHEM(N, T, Y, P) |
|---|
| 1145 | INCLUDE 'KPP_ROOT_Parameters.h' |
|---|
| 1146 | INCLUDE 'KPP_ROOT_Global.h' |
|---|
| 1147 | INTEGER N |
|---|
| 1148 | KPP_REAL T, Told |
|---|
| 1149 | KPP_REAL Y(NVAR), P(NVAR) |
|---|
| 1150 | Told = TIME |
|---|
| 1151 | TIME = T |
|---|
| 1152 | CALL Update_SUN() |
|---|
| 1153 | CALL Update_RCONST() |
|---|
| 1154 | CALL Fun( Y, FIX, RCONST, P ) |
|---|
| 1155 | TIME = Told |
|---|
| 1156 | RETURN |
|---|
| 1157 | END |
|---|
| 1158 | |
|---|
| 1159 | |
|---|
| 1160 | SUBROUTINE JAC_CHEM(N, T, Y, J) |
|---|
| 1161 | INCLUDE 'KPP_ROOT_Parameters.h' |
|---|
| 1162 | INCLUDE 'KPP_ROOT_Global.h' |
|---|
| 1163 | INTEGER N |
|---|
| 1164 | KPP_REAL Told, T |
|---|
| 1165 | KPP_REAL Y(NVAR), J(LU_NONZERO) |
|---|
| 1166 | Told = TIME |
|---|
| 1167 | TIME = T |
|---|
| 1168 | CALL Update_SUN() |
|---|
| 1169 | CALL Update_RCONST() |
|---|
| 1170 | CALL Jac_SP( Y, FIX, RCONST, J ) |
|---|
| 1171 | TIME = Told |
|---|
| 1172 | RETURN |
|---|
| 1173 | END |
|---|
| 1174 | |
|---|
