| 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 (LWORK=2*NVAR*NVAR+12*NVAR+7,LIWORK=2*NVAR+7) |
|---|
| 14 | PARAMETER (LRCONT=5*NVAR+2) |
|---|
| 15 | |
|---|
| 16 | KPP_REAL WORK(LWORK) |
|---|
| 17 | INTEGER IWORK(LIWORK) |
|---|
| 18 | COMMON /CONT/ ICONT(4),RCONT(LRCONT) |
|---|
| 19 | EXTERNAL FUNC_CHEM,JAC_CHEM |
|---|
| 20 | |
|---|
| 21 | ITOL=1 ! --- VECTOR TOLERANCES |
|---|
| 22 | IJAC=1 ! --- COMPUTE THE JACOBIAN ANALYTICALLY |
|---|
| 23 | IMAS=0 ! --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM |
|---|
| 24 | |
|---|
| 25 | DO i=1,20 |
|---|
| 26 | IWORK(i) = 0 |
|---|
| 27 | WORK(i) = 0.D0 |
|---|
| 28 | ENDDO |
|---|
| 29 | |
|---|
| 30 | IWORK(3) = 8 |
|---|
| 31 | |
|---|
| 32 | CALL ATMSDIRK(NVAR,FUNC_CHEM,TIN,VAR,TOUT,STEPMIN, |
|---|
| 33 | & RTOL,ATOL,ITOL, |
|---|
| 34 | & JAC_CHEM ,IJAC, FUNC_CHEM ,IMAS, |
|---|
| 35 | & WORK,LWORK,IWORK,LIWORK,LRCONT,IDID) |
|---|
| 36 | |
|---|
| 37 | IF (IDID.LT.0) THEN |
|---|
| 38 | print *,'ATMSDIRK: Unsucessfull exit at T=', |
|---|
| 39 | & TIN,' (IDID=',IDID,')' |
|---|
| 40 | ENDIF |
|---|
| 41 | |
|---|
| 42 | RETURN |
|---|
| 43 | END |
|---|
| 44 | |
|---|
| 45 | |
|---|
| 46 | SUBROUTINE ATMSDIRK(N,FCN,X,Y,XEND,H, |
|---|
| 47 | & RelTol,AbsTol,ITOL, |
|---|
| 48 | & JAC ,IJAC, MAS ,IMAS, |
|---|
| 49 | & WORK,LWORK,IWORK,LIWORK,LRCONT,IDID) |
|---|
| 50 | C ---------------------------------------------------------- |
|---|
| 51 | C *** *** *** *** *** *** *** *** *** *** *** *** *** |
|---|
| 52 | C DECLARATIONS |
|---|
| 53 | C *** *** *** *** *** *** *** *** *** *** *** *** *** |
|---|
| 54 | IMPLICIT KPP_REAL (A-H,O-Z) |
|---|
| 55 | DIMENSION Y(N),AbsTol(1),RelTol(1),WORK(LWORK),IWORK(LIWORK) |
|---|
| 56 | LOGICAL IMPLCT,JBAND,ARRET |
|---|
| 57 | EXTERNAL FCN,JAC,MAS |
|---|
| 58 | COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL |
|---|
| 59 | |
|---|
| 60 | C *** *** *** *** *** *** *** |
|---|
| 61 | C SETTING THE PARAMETERS |
|---|
| 62 | C *** *** *** *** *** *** *** |
|---|
| 63 | NFCN=0 |
|---|
| 64 | NJAC=0 |
|---|
| 65 | NSTEP=0 |
|---|
| 66 | NACCPT=0 |
|---|
| 67 | NREJCT=0 |
|---|
| 68 | NDEC=0 |
|---|
| 69 | NSOL=0 |
|---|
| 70 | ARRET=.FALSE. |
|---|
| 71 | C -------- SWITCH FOR TRANSFORMATION OF JACOBIAN TO HESS_CHEM FORM --- |
|---|
| 72 | NHESS1 = 0 ! ADRIAN |
|---|
| 73 | C -------- NMAX , THE MAXIMAL NUMBER OF STEPS ----- |
|---|
| 74 | IF(IWORK(2).EQ.0)THEN |
|---|
| 75 | NMAX=100000 |
|---|
| 76 | ELSE |
|---|
| 77 | NMAX=IWORK(2) |
|---|
| 78 | IF(NMAX.LE.0)THEN |
|---|
| 79 | WRITE(6,*)' WRONG INPUT IWORK(2)=',IWORK(2) |
|---|
| 80 | ARRET=.TRUE. |
|---|
| 81 | END IF |
|---|
| 82 | END IF |
|---|
| 83 | C -------- NIT MAXIMAL NUMBER OF NEWTON ITERATIONS |
|---|
| 84 | IF(IWORK(3).EQ.0)THEN |
|---|
| 85 | NIT=8 |
|---|
| 86 | ELSE |
|---|
| 87 | NIT=IWORK(3) |
|---|
| 88 | IF(NIT.LE.0)THEN |
|---|
| 89 | WRITE(6,*)' CURIOUS INPUT IWORK(3)=',IWORK(3) |
|---|
| 90 | ARRET=.TRUE. |
|---|
| 91 | END IF |
|---|
| 92 | END IF |
|---|
| 93 | C -------- METH SWITCH FOR THE COEFFICIENTS OF THE METHOD |
|---|
| 94 | METH = 2 |
|---|
| 95 | C -------- UROUND SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.D0 |
|---|
| 96 | IF(WORK(1).EQ.0.D0)THEN |
|---|
| 97 | UROUND=1.D-16 |
|---|
| 98 | ELSE |
|---|
| 99 | UROUND=WORK(1) |
|---|
| 100 | IF(UROUND.LE.1.D-19.OR.UROUND.GE.1.D0)THEN |
|---|
| 101 | WRITE(6,*)' COEFFICIENTS HAVE 20 DIGITS, UROUND=',WORK(1) |
|---|
| 102 | ARRET=.TRUE. |
|---|
| 103 | END IF |
|---|
| 104 | END IF |
|---|
| 105 | C --------- SAFE SAFETY FACTOR IN STEP SIZE PREDICTION |
|---|
| 106 | IF(WORK(2).EQ.0.D0)THEN |
|---|
| 107 | SAFE=0.9D0 |
|---|
| 108 | ELSE |
|---|
| 109 | SAFE=WORK(2) |
|---|
| 110 | IF(SAFE.LE..001D0.OR.SAFE.GE.1.D0)THEN |
|---|
| 111 | WRITE(6,*)' CURIOUS INPUT FOR WORK(2)=',WORK(2) |
|---|
| 112 | ARRET=.TRUE. |
|---|
| 113 | END IF |
|---|
| 114 | END IF |
|---|
| 115 | C ------ THET DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED; |
|---|
| 116 | IF(WORK(3).EQ.0.D0)THEN |
|---|
| 117 | THET=0.001D0 |
|---|
| 118 | ELSE |
|---|
| 119 | THET=WORK(3) |
|---|
| 120 | END IF |
|---|
| 121 | C --- FNEWT STOPPING CRIERION FOR NEWTON'S METHOD, USUALLY CHOSEN <1. |
|---|
| 122 | IF(WORK(4).EQ.0.D0)THEN |
|---|
| 123 | FNEWT=0.03D0 |
|---|
| 124 | ELSE |
|---|
| 125 | FNEWT=WORK(4) |
|---|
| 126 | END IF |
|---|
| 127 | C --- QUOT1 AND QUOT2: IF QUOT1 < HNEW/HOLD < QUOT2, STEP SIZE = CONST. |
|---|
| 128 | IF(WORK(5).EQ.0.D0)THEN |
|---|
| 129 | QUOT1=1.D0 |
|---|
| 130 | ELSE |
|---|
| 131 | QUOT1=WORK(5) |
|---|
| 132 | END IF |
|---|
| 133 | IF(WORK(6).EQ.0.D0)THEN |
|---|
| 134 | QUOT2=1.2D0 |
|---|
| 135 | ELSE |
|---|
| 136 | QUOT2=WORK(6) |
|---|
| 137 | END IF |
|---|
| 138 | C -------- MAXIMAL STEP SIZE |
|---|
| 139 | IF(WORK(7).EQ.0.D0)THEN |
|---|
| 140 | HMAX=XEND-X |
|---|
| 141 | ELSE |
|---|
| 142 | HMAX=WORK(7) |
|---|
| 143 | END IF |
|---|
| 144 | C --------- CHECK IF TOLERANCES ARE O.K. |
|---|
| 145 | IF (ITOL.EQ.0) THEN |
|---|
| 146 | IF (AbsTol(1).LE.0.D0.OR.RelTol(1).LE.10.D0*UROUND) THEN |
|---|
| 147 | WRITE (6,*) ' TOLERANCES ARE TOO SMALL' |
|---|
| 148 | ARRET=.TRUE. |
|---|
| 149 | END IF |
|---|
| 150 | ELSE |
|---|
| 151 | DO 15 I=1,N |
|---|
| 152 | IF (AbsTol(I).LE.0.D0.OR.RelTol(I).LE.10.D0*UROUND) THEN |
|---|
| 153 | WRITE (6,*) ' TOLERANCES(',I,') ARE TOO SMALL' |
|---|
| 154 | ARRET=.TRUE. |
|---|
| 155 | END IF |
|---|
| 156 | 15 CONTINUE |
|---|
| 157 | END IF |
|---|
| 158 | |
|---|
| 159 | C *** *** *** *** *** *** *** *** *** *** *** *** *** |
|---|
| 160 | C COMPUTATION OF ARRAY ENTRIES |
|---|
| 161 | C *** *** *** *** *** *** *** *** *** *** *** *** *** |
|---|
| 162 | C ---- IMPLICIT, BANDED OR NOT ? |
|---|
| 163 | IMPLCT=IMAS.NE.0 |
|---|
| 164 | ARRET=.FALSE. |
|---|
| 165 | C -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS --- |
|---|
| 166 | C -- JACOBIAN |
|---|
| 167 | LDJAC=N |
|---|
| 168 | C -- MATRIX E FOR LINEAR ALGEBRA |
|---|
| 169 | LDE=N |
|---|
| 170 | C -- MASS MATRIX |
|---|
| 171 | IF (IMPLCT) THEN |
|---|
| 172 | print *,'IMPLCT 1' |
|---|
| 173 | ELSE |
|---|
| 174 | LDMAS=0 |
|---|
| 175 | END IF |
|---|
| 176 | LDMAS2=MAX(1,LDMAS) |
|---|
| 177 | |
|---|
| 178 | C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK ----- |
|---|
| 179 | IEYHAT=8 |
|---|
| 180 | IEZ=IEYHAT+N |
|---|
| 181 | IEY0=IEZ+N |
|---|
| 182 | IEZ1=IEY0+N |
|---|
| 183 | IEZ2=IEZ1+N |
|---|
| 184 | IEZ3=IEZ2+N |
|---|
| 185 | IEZ4=IEZ3+N |
|---|
| 186 | IEZ5=IEZ4+N |
|---|
| 187 | IESCAL=IEZ5+N |
|---|
| 188 | IEF1=IESCAL+N |
|---|
| 189 | IEG1=IEF1+N |
|---|
| 190 | IEH1=IEG1+N |
|---|
| 191 | IEJAC=IEH1+N |
|---|
| 192 | IEMAS=IEJAC+N*LDJAC |
|---|
| 193 | IEE=IEMAS+N*LDMAS |
|---|
| 194 | |
|---|
| 195 | C ------ TOTAL STORAGE REQUIREMENT ----------- |
|---|
| 196 | ISTORE=IEE+N*LDE-1 |
|---|
| 197 | IF(ISTORE.GT.LWORK)THEN |
|---|
| 198 | WRITE(6,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE |
|---|
| 199 | ARRET=.TRUE. |
|---|
| 200 | END IF |
|---|
| 201 | C ------- ENTRY POINTS FOR INTEGER WORKSPACE ----- |
|---|
| 202 | IEIP=5 |
|---|
| 203 | IEHES=IEIP+N |
|---|
| 204 | C --------- TOTAL REQUIREMENT --------------- |
|---|
| 205 | ISTORE=IEHES+N-1 |
|---|
| 206 | IF(ISTORE.GT.LIWORK)THEN |
|---|
| 207 | WRITE(6,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE |
|---|
| 208 | ARRET=.TRUE. |
|---|
| 209 | END IF |
|---|
| 210 | C --------- CONTROL OF LENGTH OF COMMON BLOCK "CONT" ------- |
|---|
| 211 | IF(LRCONT.LT.(5*N+2))THEN |
|---|
| 212 | WRITE(6,*)' INSUFF. STORAGE FOR RCONT, MIN. LRCONT=',5*N+2 |
|---|
| 213 | ARRET=.TRUE. |
|---|
| 214 | END IF |
|---|
| 215 | C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1 |
|---|
| 216 | IF (ARRET) THEN |
|---|
| 217 | IDID=-1 |
|---|
| 218 | RETURN |
|---|
| 219 | END IF |
|---|
| 220 | C -------- CALL TO CORE INTEGRATOR ------------ |
|---|
| 221 | CALL SDICOR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,ITOL, |
|---|
| 222 | & JAC,IJAC,MLJAC,MUJAC,MAS,MLMAS,MUMAS,IOUT,IDID, |
|---|
| 223 | & NMAX,UROUND,SAFE,THET,FNEWT,QUOT1,QUOT2,NIT,METH,NHESS1, |
|---|
| 224 | & IMPLCT,JBAND,LDJAC,LDE,LDMAS2, |
|---|
| 225 | & WORK(IEYHAT),WORK(IEZ),WORK(IEY0),WORK(IEZ1),WORK(IEZ2), |
|---|
| 226 | & WORK(IEZ3),WORK(IEZ4),WORK(IEZ5),WORK(IESCAL),WORK(IEF1), |
|---|
| 227 | & WORK(IEG1),WORK(IEH1),WORK(IEJAC),WORK(IEE), |
|---|
| 228 | & WORK(IEMAS),IWORK(IEIP),IWORK(IEHES)) |
|---|
| 229 | C ----------- RETURN ----------- |
|---|
| 230 | RETURN |
|---|
| 231 | END |
|---|
| 232 | C |
|---|
| 233 | C |
|---|
| 234 | C |
|---|
| 235 | C ----- ... AND HERE IS THE CORE INTEGRATOR ---------- |
|---|
| 236 | C |
|---|
| 237 | SUBROUTINE SDICOR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,ITOL, |
|---|
| 238 | & JAC,IJAC,MLJAC,MUJAC,MAS,MLMAS,MUMAS,IOUT,IDID, |
|---|
| 239 | & NMAX,UROUND,SAFE,THET,FNEWT,QUOT1,QUOT2,NIT,METH,NHESS1, |
|---|
| 240 | & IMPLCT,BANDED,LDJAC,LE,LDMAS, |
|---|
| 241 | & YHAT,Z,Y0,Z1,Z2,Z3,Z4,Z5,SCAL,F1,G1,H1,FJAC,E,FMAS,IP,IPHES) |
|---|
| 242 | C ---------------------------------------------------------- |
|---|
| 243 | C CORE INTEGRATOR FOR SDIRK4 |
|---|
| 244 | C PARAMETERS SAME AS IN SDIRK4 WITH WORKSPACE ADDED |
|---|
| 245 | C ---------------------------------------------------------- |
|---|
| 246 | C DECLARATIONS |
|---|
| 247 | C ---------------------------------------------------------- |
|---|
| 248 | IMPLICIT KPP_REAL (A-H,O-Z) |
|---|
| 249 | INCLUDE 'KPP_ROOT_Parameters.h' |
|---|
| 250 | INCLUDE 'KPP_ROOT_Sparse.h' |
|---|
| 251 | KPP_REAL Y(N),YHAT(N),Z(N),Y0(N),Z1(N),Z2(N),Z3(N),Z4(N),Z5(N) |
|---|
| 252 | KPP_REAL SCAL(N),F1(N),G1(N),H1(N) |
|---|
| 253 | KPP_REAL FJAC(LU_NONZERO),E(LU_NONZERO),FMAS(LDMAS,N) |
|---|
| 254 | KPP_REAL AbsTol(1),RelTol(1) |
|---|
| 255 | INTEGER IP(N),IPHES(N) |
|---|
| 256 | LOGICAL REJECT,FIRST,IMPLCT,BANDED,CALJAC,NEWTRE |
|---|
| 257 | COMMON /CONT/NN,NN2,NN3,NN4,XOLD,HSOL,CONT(5*NVAR) |
|---|
| 258 | COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL |
|---|
| 259 | EXTERNAL MAS, FCN, JAC |
|---|
| 260 | |
|---|
| 261 | C *** *** *** *** *** *** *** |
|---|
| 262 | C INITIALISATIONS |
|---|
| 263 | C *** *** *** *** *** *** *** |
|---|
| 264 | |
|---|
| 265 | C --------- DUPLIFY N FOR COMMON BLOCK CONT ----- |
|---|
| 266 | NN=N |
|---|
| 267 | NN2=2*N |
|---|
| 268 | NN3=3*N |
|---|
| 269 | NN4=4*N |
|---|
| 270 | |
|---|
| 271 | C ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ---------- |
|---|
| 272 | IF(IMPLCT) CALL MAS(N,FMAS,LDMAS) |
|---|
| 273 | |
|---|
| 274 | C ---------- CONSTANTS --------- |
|---|
| 275 | MBDIAG=MUMAS+1 |
|---|
| 276 | IF (METH.EQ.2) THEN |
|---|
| 277 | C ---------- METHOD WITH GAMMA = 4/15 --------------- |
|---|
| 278 | GAMMA=4.0D0/15.0D0 |
|---|
| 279 | C2=23.0D0/30.0D0 |
|---|
| 280 | C3=17.0D0/30.0D0 |
|---|
| 281 | C4=2881.0D0/28965.0D0+GAMMA |
|---|
| 282 | ALPH21=15.0D0/8.0D0 |
|---|
| 283 | ALPH31=1577061.0D0/922880.0D0 |
|---|
| 284 | ALPH32=-23427.0D0/115360.0D0 |
|---|
| 285 | ALPH41=647163682356923881.0D0/2414496535205978880.0D0 |
|---|
| 286 | ALPH42=-593512117011179.0D0/3245291041943520.0D0 |
|---|
| 287 | ALPH43=559907973726451.0D0/1886325418129671.0D0 |
|---|
| 288 | ALPH51=724545451.0D0/796538880.0D0 |
|---|
| 289 | ALPH52=-830832077.0D0/267298560.0D0 |
|---|
| 290 | ALPH53=30957577.0D0/2509272.0D0 |
|---|
| 291 | ALPH54=-69863904375173.0D0/6212571137048.0D0 |
|---|
| 292 | E1=7752107607.0D0/11393456128.0D0 |
|---|
| 293 | E2=-17881415427.0D0/11470078208.0D0 |
|---|
| 294 | E3=2433277665.0D0/179459416.0D0 |
|---|
| 295 | E4=-96203066666797.0D0/6212571137048.0D0 |
|---|
| 296 | D11= 24.74416644927758D0 |
|---|
| 297 | D12= -4.325375951824688D0 |
|---|
| 298 | D13= 41.39683763286316D0 |
|---|
| 299 | D14= -61.04144619901784D0 |
|---|
| 300 | D15= -3.391332232917013D0 |
|---|
| 301 | D21= -51.98245719616925D0 |
|---|
| 302 | D22= 10.52501981094525D0 |
|---|
| 303 | D23= -154.2067922191855D0 |
|---|
| 304 | D24= 214.3082125319825D0 |
|---|
| 305 | D25= 14.71166018088679D0 |
|---|
| 306 | D31= 33.14347947522142D0 |
|---|
| 307 | D32= -19.72986789558523D0 |
|---|
| 308 | D33= 230.4878502285804D0 |
|---|
| 309 | D34= -287.6629744338197D0 |
|---|
| 310 | D35= -18.99932366302254D0 |
|---|
| 311 | D41= -5.905188728329743D0 |
|---|
| 312 | D42= 13.53022403646467D0 |
|---|
| 313 | D43= -117.6778956422581D0 |
|---|
| 314 | D44= 134.3962081008550D0 |
|---|
| 315 | D45= 8.678995715052762D0 |
|---|
| 316 | ETA1=23.D0/8.D0 |
|---|
| 317 | ANU1= 0.9838473040915402D0 |
|---|
| 318 | ANU2= 0.3969226768377252D0 |
|---|
| 319 | AMU1= 0.6563374010466914D0 |
|---|
| 320 | AMU3= 0.3372498196189311D0 |
|---|
| 321 | ELSE |
|---|
| 322 | PRINT *, 'WRONG CHOICE OF <METH>' |
|---|
| 323 | END IF |
|---|
| 324 | POSNEG=SIGN(1.D0,XEND-X) |
|---|
| 325 | HMAX1=MIN(ABS(HMAX),ABS(XEND-X)) |
|---|
| 326 | IF (ABS(H).LE.10.D0*UROUND) H=1.0D-6 |
|---|
| 327 | H=MIN(ABS(H),HMAX1) |
|---|
| 328 | H=SIGN(H,POSNEG) |
|---|
| 329 | HOLD=H |
|---|
| 330 | CFAC=SAFE*(1+2*NIT) |
|---|
| 331 | NEWTRE=.FALSE. |
|---|
| 332 | REJECT=.FALSE. |
|---|
| 333 | FIRST=.TRUE. |
|---|
| 334 | FACCO1=1.D0 |
|---|
| 335 | FACCO2=1.D0 |
|---|
| 336 | FACCO3=1.D0 |
|---|
| 337 | FACCO4=1.D0 |
|---|
| 338 | FACCO5=1.D0 |
|---|
| 339 | NSING=0 |
|---|
| 340 | XOLD=X |
|---|
| 341 | IF (ITOL.EQ.0) THEN |
|---|
| 342 | DO 8 I=1,N |
|---|
| 343 | 8 SCAL(I)=1.D0 / ( AbsTol(1)+RelTol(1)*DABS(Y(I)) ) |
|---|
| 344 | ELSE |
|---|
| 345 | DO 9 I=1,N |
|---|
| 346 | 9 SCAL(I)=1.D0 / ( AbsTol(I)+RelTol(I)*DABS(Y(I)) ) |
|---|
| 347 | END IF |
|---|
| 348 | |
|---|
| 349 | C --- BASIC INTEGRATION STEP |
|---|
| 350 | 10 CONTINUE |
|---|
| 351 | |
|---|
| 352 | C *** *** *** *** *** *** *** |
|---|
| 353 | C COMPUTATION OF THE JACOBIAN |
|---|
| 354 | C *** *** *** *** *** *** *** |
|---|
| 355 | NJAC=NJAC+1 |
|---|
| 356 | CALL JAC(N,X,Y,FJAC) |
|---|
| 357 | CALJAC=.TRUE. |
|---|
| 358 | 20 CONTINUE |
|---|
| 359 | |
|---|
| 360 | C *** *** *** *** *** *** *** |
|---|
| 361 | C COMPUTE THE MATRIX E AND ITS DECOMPOSITION |
|---|
| 362 | C *** *** *** *** *** *** *** |
|---|
| 363 | FAC1=1.D0/(H*GAMMA) |
|---|
| 364 | IF (IMPLCT) THEN |
|---|
| 365 | print *, 'IMPLCT 4' |
|---|
| 366 | ELSE ! EXPLICIT SYSTEM |
|---|
| 367 | C --- THE MATRIX E (MAS=IDENTITY, JACOBIAN A FULL MATRIX) |
|---|
| 368 | c DO 526 J=1,N |
|---|
| 369 | c DO 525 I=1,N |
|---|
| 370 | c 525 E(I,J)=-FJAC(I,J) |
|---|
| 371 | c 526 E(J,J)=E(J,J)+FAC1 |
|---|
| 372 | c CALL DEC(N,LE,E,IP,IER) |
|---|
| 373 | DO K=1,LU_NONZERO |
|---|
| 374 | E(K) = -FJAC(K) |
|---|
| 375 | END DO |
|---|
| 376 | DO I=1,N |
|---|
| 377 | IDG = LU_DIAG(I) |
|---|
| 378 | E(IDG) = E(IDG) + FAC1 |
|---|
| 379 | END DO |
|---|
| 380 | CALL KppDecomp ( E, IER) |
|---|
| 381 | |
|---|
| 382 | IF (IER.NE.0) GOTO 79 |
|---|
| 383 | END IF |
|---|
| 384 | NDEC=NDEC+1 |
|---|
| 385 | 30 CONTINUE |
|---|
| 386 | |
|---|
| 387 | IF (NSTEP.GT.NMAX.OR.X+.1D0*H.EQ.X.OR.ABS(H).LE.UROUND) GOTO 79 |
|---|
| 388 | XPH=X+H |
|---|
| 389 | C --- LOOP FOR THE 5 STAGES |
|---|
| 390 | FACCO1=DMAX1(FACCO1,UROUND)**0.8D0 |
|---|
| 391 | FACCO2=DMAX1(FACCO2,UROUND)**0.8D0 |
|---|
| 392 | FACCO3=DMAX1(FACCO3,UROUND)**0.8D0 |
|---|
| 393 | FACCO4=DMAX1(FACCO4,UROUND)**0.8D0 |
|---|
| 394 | FACCO5=DMAX1(FACCO5,UROUND)**0.8D0 |
|---|
| 395 | |
|---|
| 396 | C *** *** *** *** *** *** *** |
|---|
| 397 | C STARTING VALUES FOR NEWTON ITERATION |
|---|
| 398 | C *** *** *** *** *** *** *** |
|---|
| 399 | DO 59 ISTAGE=1,5 |
|---|
| 400 | IF (ISTAGE.EQ.1) THEN |
|---|
| 401 | XCH=X+GAMMA*H |
|---|
| 402 | IF (FIRST.OR.NEWTRE) THEN |
|---|
| 403 | DO 132 I=1,N |
|---|
| 404 | 132 Z(I)=0.D0 |
|---|
| 405 | ELSE |
|---|
| 406 | S=1.D0+GAMMA*H/HOLD |
|---|
| 407 | DO 232 I=1,N |
|---|
| 408 | c 232 Z(I) = 0.D0 |
|---|
| 409 | 232 Z(I)=S*(CONT(I+NN)+S*(CONT(I+NN2)+S*(CONT(I+NN3) |
|---|
| 410 | & +S*CONT(I+NN4))))-YHAT(I) |
|---|
| 411 | |
|---|
| 412 | END IF |
|---|
| 413 | DO 31 I=1,N |
|---|
| 414 | 31 G1(I)=0.D0 |
|---|
| 415 | FACCON=FACCO1 |
|---|
| 416 | END IF |
|---|
| 417 | IF (ISTAGE.EQ.2) THEN |
|---|
| 418 | XCH=X+C2*H |
|---|
| 419 | DO 131 I=1,N |
|---|
| 420 | Z1I=Z1(I) |
|---|
| 421 | Z(I)=ETA1*Z1I |
|---|
| 422 | 131 G1(I)=ALPH21*Z1I |
|---|
| 423 | FACCON=FACCO2 |
|---|
| 424 | END IF |
|---|
| 425 | IF (ISTAGE.EQ.3) THEN |
|---|
| 426 | XCH=X+C3*H |
|---|
| 427 | DO 231 I=1,N |
|---|
| 428 | Z1I=Z1(I) |
|---|
| 429 | Z2I=Z2(I) |
|---|
| 430 | Z(I)=ANU1*Z1I+ANU2*Z2I |
|---|
| 431 | 231 G1(I)=ALPH31*Z1I+ALPH32*Z2I |
|---|
| 432 | FACCON=FACCO3 |
|---|
| 433 | END IF |
|---|
| 434 | IF (ISTAGE.EQ.4) THEN |
|---|
| 435 | XCH=X+C4*H |
|---|
| 436 | DO 331 I=1,N |
|---|
| 437 | Z1I=Z1(I) |
|---|
| 438 | Z3I=Z3(I) |
|---|
| 439 | Z(I)=AMU1*Z1I+AMU3*Z3I |
|---|
| 440 | 331 G1(I)=ALPH41*Z1I+ALPH42*Z2(I)+ALPH43*Z3I |
|---|
| 441 | FACCON=FACCO4 |
|---|
| 442 | END IF |
|---|
| 443 | IF (ISTAGE.EQ.5) THEN |
|---|
| 444 | XCH=XPH |
|---|
| 445 | DO 431 I=1,N |
|---|
| 446 | Z1I=Z1(I) |
|---|
| 447 | Z2I=Z2(I) |
|---|
| 448 | Z3I=Z3(I) |
|---|
| 449 | Z4I=Z4(I) |
|---|
| 450 | Z(I)=E1*Z1I+E2*Z2I+E3*Z3I+E4*Z4I |
|---|
| 451 | YHAT(I)=Z(I) |
|---|
| 452 | 431 G1(I)=ALPH51*Z1I+ALPH52*Z2I+ALPH53*Z3I+ALPH54*Z4I |
|---|
| 453 | FACCON=FACCO5 |
|---|
| 454 | END IF |
|---|
| 455 | |
|---|
| 456 | |
|---|
| 457 | |
|---|
| 458 | C *** *** *** *** *** *** *** *** *** *** *** |
|---|
| 459 | C LOOP FOR THE SIMPLIFIED NEWTON ITERATION |
|---|
| 460 | C *** *** *** *** *** *** *** *** *** *** *** |
|---|
| 461 | NEWT=0 |
|---|
| 462 | THETA=ABS(THET) |
|---|
| 463 | IF (REJECT) THETA=2*ABS(THET) |
|---|
| 464 | 40 CONTINUE |
|---|
| 465 | IF (NEWT.GE.NIT) THEN |
|---|
| 466 | H=H/2.D0 |
|---|
| 467 | REJECT=.TRUE. |
|---|
| 468 | NEWTRE=.TRUE. |
|---|
| 469 | IF (CALJAC) GOTO 20 |
|---|
| 470 | GOTO 10 |
|---|
| 471 | END IF |
|---|
| 472 | |
|---|
| 473 | C --- COMPUTE THE RIGHT-HAND SIDE |
|---|
| 474 | DO 41 I=1,N |
|---|
| 475 | H1(I)=G1(I)-Z(I) |
|---|
| 476 | 41 CONT(I)=Y(I)+Z(I) |
|---|
| 477 | CALL FCN(N,XCH,CONT,F1) |
|---|
| 478 | NFCN=NFCN+1 |
|---|
| 479 | |
|---|
| 480 | C --- KppSolve THE LINEAR SYSTEMS |
|---|
| 481 | IF (IMPLCT) THEN |
|---|
| 482 | print *, 'IMPLCT 2' |
|---|
| 483 | ELSE |
|---|
| 484 | DO 345 I=1,N |
|---|
| 485 | 345 F1(I)=H1(I)*FAC1+F1(I) |
|---|
| 486 | C CALL SOL(N,LE,E,F1,IP) |
|---|
| 487 | CALL KppSolve(E, F1) |
|---|
| 488 | END IF |
|---|
| 489 | NEWT=NEWT+1 |
|---|
| 490 | DYNO=0.D0 |
|---|
| 491 | C --- NORM 2 --- |
|---|
| 492 | DO 57 I=1,N |
|---|
| 493 | 57 DYNO=DYNO+(F1(I)*SCAL(I))**2 |
|---|
| 494 | DYNO=DSQRT(DYNO/N) |
|---|
| 495 | C --- NORM INF --- |
|---|
| 496 | C DO 57 I=1,N |
|---|
| 497 | C 57 DYNO=DMAX1( DYNO, DABS(F1(I)*SCAL(I)) ) |
|---|
| 498 | |
|---|
| 499 | |
|---|
| 500 | C --- BAD CONVERGENCE OR NUMBER OF ITERATIONS TO LARGE |
|---|
| 501 | IF (NEWT.GE.2.AND.NEWT.LT.NIT) THEN |
|---|
| 502 | THETA=DYNO/DYNOLD |
|---|
| 503 | IF (THETA.LT.0.99D0) THEN |
|---|
| 504 | FACCON=THETA/(1.0D0-THETA) |
|---|
| 505 | DYTH=FACCON*DYNO*THETA**(NIT-1-NEWT) |
|---|
| 506 | QNEWT=DMAX1(1.0D-4,DMIN1(16.0D0,DYTH/FNEWT)) |
|---|
| 507 | IF (QNEWT.GE.1.0D0) THEN |
|---|
| 508 | H=.8D0*H*QNEWT**(-1.0D0/(NIT-NEWT)) |
|---|
| 509 | REJECT=.TRUE. |
|---|
| 510 | NEWTRE=.TRUE. |
|---|
| 511 | IF (CALJAC) GOTO 20 |
|---|
| 512 | GOTO 10 |
|---|
| 513 | END IF |
|---|
| 514 | ELSE |
|---|
| 515 | NEWTRE=.TRUE. |
|---|
| 516 | GOTO 78 |
|---|
| 517 | END IF |
|---|
| 518 | END IF |
|---|
| 519 | DYNOLD=DYNO |
|---|
| 520 | DO 58 I=1,N |
|---|
| 521 | 58 Z(I)=Z(I)+F1(I) |
|---|
| 522 | NSOL=NSOL+1 |
|---|
| 523 | IF (FACCON*DYNO.GT.FNEWT) GOTO 40 |
|---|
| 524 | |
|---|
| 525 | C --- END OF SIMPILFIED NEWTON |
|---|
| 526 | IF (ISTAGE.EQ.1) THEN |
|---|
| 527 | DO I=1,N |
|---|
| 528 | Z1(I) = Z(I) |
|---|
| 529 | END DO |
|---|
| 530 | FACCO1=FACCON |
|---|
| 531 | END IF |
|---|
| 532 | IF (ISTAGE.EQ.2) THEN |
|---|
| 533 | DO I=1,N |
|---|
| 534 | Z2(I) = Z(I) |
|---|
| 535 | END DO |
|---|
| 536 | FACCO2=FACCON |
|---|
| 537 | END IF |
|---|
| 538 | IF (ISTAGE.EQ.3) THEN |
|---|
| 539 | DO I=1,N |
|---|
| 540 | Z3(I) = Z(I) |
|---|
| 541 | END DO |
|---|
| 542 | FACCO3=FACCON |
|---|
| 543 | END IF |
|---|
| 544 | IF (ISTAGE.EQ.4) THEN |
|---|
| 545 | DO I=1,N |
|---|
| 546 | Z4(I) = Z(I) |
|---|
| 547 | END DO |
|---|
| 548 | FACCO4=FACCON |
|---|
| 549 | END IF |
|---|
| 550 | IF (ISTAGE.EQ.5) THEN |
|---|
| 551 | DO I=1,N |
|---|
| 552 | Z5(I) = Z(I) |
|---|
| 553 | END DO |
|---|
| 554 | FACCO5=FACCON |
|---|
| 555 | END IF |
|---|
| 556 | 59 CONTINUE |
|---|
| 557 | |
|---|
| 558 | |
|---|
| 559 | C *** *** *** *** *** *** *** |
|---|
| 560 | C ERROR ESTIMATION |
|---|
| 561 | C *** *** *** *** *** *** *** |
|---|
| 562 | NSTEP=NSTEP+1 |
|---|
| 563 | IF (IMPLCT) THEN |
|---|
| 564 | print *,'IMPLCT 3' |
|---|
| 565 | ELSE |
|---|
| 566 | DO 461 I=1,N |
|---|
| 567 | 461 CONT(I)=FAC1*(Z5(I)-YHAT(I)) |
|---|
| 568 | END IF |
|---|
| 569 | |
|---|
| 570 | CALL KppSolve(E, CONT) |
|---|
| 571 | |
|---|
| 572 | ERR=0.D0 |
|---|
| 573 | C ---- NORM 2 --- |
|---|
| 574 | DO 64 I=1,N |
|---|
| 575 | 64 ERR=ERR+(CONT(I)*SCAL(I))**2 |
|---|
| 576 | ERR=DMAX1(DSQRT(ERR/N),1.D-10) |
|---|
| 577 | |
|---|
| 578 | C ---- NORM INF --- |
|---|
| 579 | C DO 64 I=1,N |
|---|
| 580 | c 64 ERR=DMAX1( ERR, DABS( CONT(I)*SCAL(I) ) ) |
|---|
| 581 | |
|---|
| 582 | C --- COMPUTATION OF HNEW |
|---|
| 583 | C --- WE REQUIRE .25<=HNEW/H<=10. |
|---|
| 584 | FAC=DMIN1(SAFE,CFAC/(NEWT+2*NIT)) |
|---|
| 585 | QUOT=DMAX1(.25D0,DMIN1(10.D0,(ERR)**.25D0/FAC)) |
|---|
| 586 | HNEW= H/QUOT |
|---|
| 587 | |
|---|
| 588 | C *** *** *** *** *** *** *** |
|---|
| 589 | C IS THE ERROR SMALL ENOUGH ? |
|---|
| 590 | C *** *** *** *** *** *** *** |
|---|
| 591 | IF (ERR.LT.1.D0) THEN |
|---|
| 592 | C --- STEP IS ACCEPTED |
|---|
| 593 | FIRST=.FALSE. |
|---|
| 594 | NACCPT=NACCPT+1 |
|---|
| 595 | HOLD=H |
|---|
| 596 | XOLD=X |
|---|
| 597 | C --- COEFFICIENTS FOR CONTINUOUS SOLUTION |
|---|
| 598 | DO 74 I=1,N |
|---|
| 599 | Z1I=Z1(I) |
|---|
| 600 | Z2I=Z2(I) |
|---|
| 601 | Z3I=Z3(I) |
|---|
| 602 | Z4I=Z4(I) |
|---|
| 603 | Z5I=Z5(I) |
|---|
| 604 | CONT(I)=Y(I) |
|---|
| 605 | Y(I)=Y(I)+Z5I |
|---|
| 606 | CONT(I+NN) =D11*Z1I+D12*Z2I+D13*Z3I+D14*Z4I+D15*Z5I |
|---|
| 607 | CONT(I+NN2)=D21*Z1I+D22*Z2I+D23*Z3I+D24*Z4I+D25*Z5I |
|---|
| 608 | CONT(I+NN3)=D31*Z1I+D32*Z2I+D33*Z3I+D34*Z4I+D35*Z5I |
|---|
| 609 | CONT(I+NN4)=D41*Z1I+D42*Z2I+D43*Z3I+D44*Z4I+D45*Z5I |
|---|
| 610 | YHAT(I)=Z5I |
|---|
| 611 | IF (ITOL.EQ.0) THEN |
|---|
| 612 | SCAL(I)=1.D0/( AbsTol(1)+RelTol(1)*DABS(Y(I)) ) |
|---|
| 613 | ELSE |
|---|
| 614 | SCAL(I)=1.D0/( AbsTol(I)+RelTol(I)*DABS(Y(I)) ) |
|---|
| 615 | END IF |
|---|
| 616 | 74 CONTINUE |
|---|
| 617 | X=XPH |
|---|
| 618 | CALJAC=.FALSE. |
|---|
| 619 | IF ((X-XEND)*POSNEG+UROUND.GT.0.D0) THEN |
|---|
| 620 | H=HOPT |
|---|
| 621 | IDID=1 |
|---|
| 622 | RETURN |
|---|
| 623 | END IF |
|---|
| 624 | IF (IJAC.EQ.0) CALL FCN(N,X,Y,Y0) |
|---|
| 625 | NFCN=NFCN+1 |
|---|
| 626 | HNEW=POSNEG*DMIN1(DABS(HNEW),HMAX1) |
|---|
| 627 | HOPT=HNEW |
|---|
| 628 | IF (REJECT) HNEW=POSNEG*DMIN1(DABS(HNEW),DABS(H)) |
|---|
| 629 | REJECT=.FALSE. |
|---|
| 630 | NEWTRE=.FALSE. |
|---|
| 631 | IF ((X+HNEW/QUOT1-XEND)*POSNEG.GT.0.D0) THEN |
|---|
| 632 | H=XEND-X |
|---|
| 633 | ELSE |
|---|
| 634 | QT=HNEW/H |
|---|
| 635 | IF (THETA.LE.THET.AND.QT.GE.QUOT1.AND.QT.LE.QUOT2) GOTO 30 |
|---|
| 636 | H = HNEW |
|---|
| 637 | END IF |
|---|
| 638 | IF (THETA.LE.THET) GOTO 20 |
|---|
| 639 | GOTO 10 |
|---|
| 640 | |
|---|
| 641 | ELSE |
|---|
| 642 | C --- STEP IS REJECTED |
|---|
| 643 | REJECT=.TRUE. |
|---|
| 644 | IF (FIRST) THEN |
|---|
| 645 | H=H/10.D0 |
|---|
| 646 | ELSE |
|---|
| 647 | H=HNEW |
|---|
| 648 | END IF |
|---|
| 649 | IF (NACCPT.GE.1) NREJCT=NREJCT+1 |
|---|
| 650 | IF (CALJAC) GOTO 20 |
|---|
| 651 | GOTO 10 |
|---|
| 652 | END IF |
|---|
| 653 | |
|---|
| 654 | C --- UNEXPECTED STEP-REJECTION |
|---|
| 655 | 78 CONTINUE |
|---|
| 656 | IF (IER.NE.0) THEN |
|---|
| 657 | WRITE (6,*) ' MATRIX IS SINGULAR, IER=',IER,' X=',X,' H=',H |
|---|
| 658 | NSING=NSING+1 |
|---|
| 659 | IF (NSING.GE.6) GOTO 79 |
|---|
| 660 | END IF |
|---|
| 661 | H=H*0.5D0 |
|---|
| 662 | REJECT=.TRUE. |
|---|
| 663 | IF (CALJAC) GOTO 20 |
|---|
| 664 | GOTO 10 |
|---|
| 665 | |
|---|
| 666 | C --- FAIL EXIT |
|---|
| 667 | 79 WRITE (6,979) X,H,IER |
|---|
| 668 | 979 FORMAT(' EXIT OF SDIRK4 AT X=',D14.7,' H=',D14.7,' IER=',I4) |
|---|
| 669 | IDID=-1 |
|---|
| 670 | RETURN |
|---|
| 671 | END |
|---|
| 672 | C |
|---|
| 673 | |
|---|
| 674 | |
|---|
| 675 | SUBROUTINE FUNC_CHEM(N, T, Y, P) |
|---|
| 676 | INCLUDE 'KPP_ROOT_Parameters.h' |
|---|
| 677 | INCLUDE 'KPP_ROOT_Global.h' |
|---|
| 678 | INTEGER N |
|---|
| 679 | KPP_REAL T, Told |
|---|
| 680 | KPP_REAL Y(NVAR), P(NVAR) |
|---|
| 681 | Told = TIME |
|---|
| 682 | TIME = T |
|---|
| 683 | CALL Update_SUN() |
|---|
| 684 | CALL Update_RCONST() |
|---|
| 685 | CALL Fun( Y, FIX, RCONST, P ) |
|---|
| 686 | TIME = Told |
|---|
| 687 | RETURN |
|---|
| 688 | END |
|---|
| 689 | |
|---|
| 690 | |
|---|
| 691 | SUBROUTINE JAC_CHEM(N, T, Y, J) |
|---|
| 692 | INCLUDE 'KPP_ROOT_Parameters.h' |
|---|
| 693 | INCLUDE 'KPP_ROOT_Global.h' |
|---|
| 694 | INTEGER N |
|---|
| 695 | KPP_REAL Told, T |
|---|
| 696 | KPP_REAL Y(NVAR), J(LU_NONZERO) |
|---|
| 697 | Told = TIME |
|---|
| 698 | TIME = T |
|---|
| 699 | CALL Update_SUN() |
|---|
| 700 | CALL Update_RCONST() |
|---|
| 701 | CALL Jac_SP( Y, FIX, RCONST, J ) |
|---|
| 702 | TIME = Told |
|---|
| 703 | RETURN |
|---|
| 704 | END |
|---|
| 705 | |
|---|
