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