| 1 | SUBROUTINE INTEGRATE( TIN, TOUT ) |
|---|
| 2 | |
|---|
| 3 | INCLUDE 'KPP_ROOT_params.h' |
|---|
| 4 | INCLUDE 'KPP_ROOT_global.h' |
|---|
| 5 | |
|---|
| 6 | C TIN - Start Time |
|---|
| 7 | KPP_REAL TIN |
|---|
| 8 | C TOUT - End Time |
|---|
| 9 | KPP_REAL TOUT |
|---|
| 10 | |
|---|
| 11 | INTEGER INFO(5) |
|---|
| 12 | |
|---|
| 13 | EXTERNAL FUNC_CHEM, JAC_CHEM |
|---|
| 14 | |
|---|
| 15 | INFO(1) = Autonomous |
|---|
| 16 | |
|---|
| 17 | CALL ROS4(NVAR,TIN,TOUT,STEPMIN,STEPMAX, |
|---|
| 18 | + STEPMIN,VAR,ATOL,RTOL, |
|---|
| 19 | + Info,FUNC_CHEM,JAC_CHEM) |
|---|
| 20 | |
|---|
| 21 | RETURN |
|---|
| 22 | END |
|---|
| 23 | |
|---|
| 24 | |
|---|
| 25 | |
|---|
| 26 | |
|---|
| 27 | SUBROUTINE ROS4(N,T,Tnext,Hmin,Hmax,Hstart, |
|---|
| 28 | + y,AbsTol,RelTol, |
|---|
| 29 | + Info,FUNC_CHEM,JAC_CHEM) |
|---|
| 30 | IMPLICIT NONE |
|---|
| 31 | INCLUDE 'KPP_ROOT_params.h' |
|---|
| 32 | INCLUDE 'KPP_ROOT_sparse.h' |
|---|
| 33 | C |
|---|
| 34 | C Four Stages, Fourth Order L-stable Rosenbrock Method, |
|---|
| 35 | C with embedded L-stable, third order method for error control |
|---|
| 36 | C Simplified version of E. Hairer's atmros4; the coefficients are slightly |
|---|
| 37 | C different |
|---|
| 38 | C |
|---|
| 39 | C |
|---|
| 40 | C INPUT ARGUMENTS: |
|---|
| 41 | C y = Vector of (NVAR) concentrations, contains the |
|---|
| 42 | C initial values on input |
|---|
| 43 | C [T, Tnext] = the integration interval |
|---|
| 44 | C Hmin, Hmax = lower and upper bounds for the selected step-size. |
|---|
| 45 | C Note that for Step = Hmin the current computed |
|---|
| 46 | C solution is unconditionally accepted by the error |
|---|
| 47 | C control mechanism. |
|---|
| 48 | C AbsTol, RelTol = (NVAR) dimensional vectors of |
|---|
| 49 | C componentwise absolute and relative tolerances. |
|---|
| 50 | C FUNC_CHEM = name of routine of derivatives. KPP syntax. |
|---|
| 51 | C See the header below. |
|---|
| 52 | C JAC_CHEM = name of routine that computes the Jacobian, in |
|---|
| 53 | C sparse format. KPP syntax. See the header below. |
|---|
| 54 | C Info(1) = 1 for Autonomous system |
|---|
| 55 | C = 0 for nonAutonomous system |
|---|
| 56 | C |
|---|
| 57 | C OUTPUT ARGUMENTS: |
|---|
| 58 | C y = the values of concentrations at Tend. |
|---|
| 59 | C T = equals TENDon output. |
|---|
| 60 | C Info(2) = # of FUNC_CHEM CALLs. |
|---|
| 61 | C Info(3) = # of JAC_CHEM CALLs. |
|---|
| 62 | C Info(4) = # of accepted steps. |
|---|
| 63 | C Info(5) = # of rejected steps. |
|---|
| 64 | C Hstart = The last accepted stepsize |
|---|
| 65 | C |
|---|
| 66 | C Adrian Sandu, December 2001 |
|---|
| 67 | C |
|---|
| 68 | KPP_REAL K1(NVAR), K2(NVAR), K3(NVAR), K4(NVAR) |
|---|
| 69 | KPP_REAL F1(NVAR) |
|---|
| 70 | KPP_REAL DFDT(NVAR) |
|---|
| 71 | KPP_REAL JAC(LU_NONZERO) |
|---|
| 72 | KPP_REAL Hmin,Hmax,Hstart |
|---|
| 73 | KPP_REAL y(NVAR), ynew(NVAR) |
|---|
| 74 | KPP_REAL AbsTol(NVAR), RelTol(NVAR) |
|---|
| 75 | KPP_REAL T, Tnext, H, Hnew, Tplus |
|---|
| 76 | KPP_REAL elo,ghinv,uround |
|---|
| 77 | KPP_REAL ERR, factor, facmax |
|---|
| 78 | KPP_REAL w, e, dround, tau |
|---|
| 79 | KPP_REAL hgam1, hgam2, hgam3, hgam4 |
|---|
| 80 | KPP_REAL hc21, hc31, hc32, hc41, hc42, hc43 |
|---|
| 81 | |
|---|
| 82 | INTEGER n,nfcn,njac,Naccept,Nreject,i,j,ier |
|---|
| 83 | INTEGER Info(5) |
|---|
| 84 | LOGICAL IsReject, Autonomous |
|---|
| 85 | EXTERNAL FUNC_CHEM, JAC_CHEM |
|---|
| 86 | |
|---|
| 87 | |
|---|
| 88 | C The method coefficients |
|---|
| 89 | DOUBLE PRECISION gamma, gamma2, gamma3, gamma4 |
|---|
| 90 | PARAMETER ( gamma = 0.5728200000000000D+00 ) |
|---|
| 91 | PARAMETER ( gamma2 = -0.1769193891319233D+01 ) |
|---|
| 92 | PARAMETER ( gamma3 = 0.7592633437920482D+00 ) |
|---|
| 93 | PARAMETER ( gamma4 = -0.1049021087100450D+00 ) |
|---|
| 94 | DOUBLE PRECISION a21, a31, a32, a41, a42, a43 |
|---|
| 95 | PARAMETER ( a21 = 0.2000000000000000D+01 ) |
|---|
| 96 | PARAMETER ( a31 = 0.1867943637803922D+01 ) |
|---|
| 97 | PARAMETER ( a32 = 0.2344449711399156D+00 ) |
|---|
| 98 | DOUBLE PRECISION alpha2, alpha3 |
|---|
| 99 | PARAMETER ( alpha2 = 0.1145640000000000D+01 ) |
|---|
| 100 | PARAMETER ( alpha3 = 0.6552168638155900D+00 ) |
|---|
| 101 | DOUBLE PRECISION c21, c31, c32, c41, c42, c43 |
|---|
| 102 | PARAMETER ( c21 = -0.7137615036412310D+01 ) |
|---|
| 103 | PARAMETER ( c31 = 0.2580708087951457D+01 ) |
|---|
| 104 | PARAMETER ( c32 = 0.6515950076447975D+00 ) |
|---|
| 105 | PARAMETER ( c41 = -0.2137148994382534D+01 ) |
|---|
| 106 | PARAMETER ( c42 = -0.3214669691237626D+00 ) |
|---|
| 107 | PARAMETER ( c43 = -0.6949742501781779D+00 ) |
|---|
| 108 | DOUBLE PRECISION b1, b2, b3, b4 |
|---|
| 109 | PARAMETER ( b1 = 0.2255570073418735D+01 ) |
|---|
| 110 | PARAMETER ( b2 = 0.2870493262186792D+00 ) |
|---|
| 111 | PARAMETER ( b3 = 0.4353179431840180D+00 ) |
|---|
| 112 | PARAMETER ( b4 = 0.1093502252409163D+01 ) |
|---|
| 113 | DOUBLE PRECISION d1, d2, d3, d4 |
|---|
| 114 | PARAMETER ( d1 = -0.2815431932141155D+00 ) |
|---|
| 115 | PARAMETER ( d2 = -0.7276199124938920D-01 ) |
|---|
| 116 | PARAMETER ( d3 = -0.1082196201495311D+00 ) |
|---|
| 117 | PARAMETER ( d4 = -0.1093502252409163D+01 ) |
|---|
| 118 | |
|---|
| 119 | c Initialization of counters, etc. |
|---|
| 120 | Autonomous = Info(1) .EQ. 1 |
|---|
| 121 | uround = 1.d-15 |
|---|
| 122 | dround = DSQRT(uround) |
|---|
| 123 | IF (Hmax.le.0.D0) THEN |
|---|
| 124 | Hmax = DABS(Tnext-T) |
|---|
| 125 | END IF |
|---|
| 126 | H = DMAX1(1.d-8, Hstart) |
|---|
| 127 | Tplus = T |
|---|
| 128 | IsReject = .false. |
|---|
| 129 | Naccept = 0 |
|---|
| 130 | Nreject = 0 |
|---|
| 131 | Nfcn = 0 |
|---|
| 132 | Njac = 0 |
|---|
| 133 | |
|---|
| 134 | C === Starting the time loop === |
|---|
| 135 | 10 CONTINUE |
|---|
| 136 | |
|---|
| 137 | Tplus = T + H |
|---|
| 138 | IF ( Tplus .gt. Tnext ) THEN |
|---|
| 139 | H = Tnext - T |
|---|
| 140 | Tplus = Tnext |
|---|
| 141 | END IF |
|---|
| 142 | |
|---|
| 143 | C Initial Function and Jacobian values |
|---|
| 144 | CALL FUNC_CHEM( T, y, F1 ) |
|---|
| 145 | CALL JAC_CHEM( T, y, JAC ) |
|---|
| 146 | |
|---|
| 147 | C The time derivative for non-Autonomous case |
|---|
| 148 | IF (.not. Autonomous) THEN |
|---|
| 149 | tau = DSIGN(dround*DMAX1( 1.0d-6, DABS(T) ), T) |
|---|
| 150 | CALL FUNC_CHEM( T+tau, y, K2 ) |
|---|
| 151 | nfcn=nfcn+1 |
|---|
| 152 | DO 20 j = 1,NVAR |
|---|
| 153 | DFDT(j) = ( K2(j)-F1(j) )/tau |
|---|
| 154 | 20 CONTINUE |
|---|
| 155 | END IF |
|---|
| 156 | |
|---|
| 157 | C Form the Prediction matrix and compute its LU factorization |
|---|
| 158 | Njac = Njac+1 |
|---|
| 159 | ghinv = 1.0d0/(gamma*H) |
|---|
| 160 | DO 30 j=1,LU_NONZERO |
|---|
| 161 | JAC(j) = -JAC(j) |
|---|
| 162 | 30 CONTINUE |
|---|
| 163 | DO 40 j=1,NVAR |
|---|
| 164 | JAC(LU_DIAG(j)) = JAC(LU_DIAG(j)) + ghinv |
|---|
| 165 | 40 CONTINUE |
|---|
| 166 | CALL KppDecomp (JAC, ier) |
|---|
| 167 | C |
|---|
| 168 | IF (ier.ne.0) THEN |
|---|
| 169 | IF ( H.gt.Hmin) THEN |
|---|
| 170 | H = 5.0d-1*H |
|---|
| 171 | GO TO 10 |
|---|
| 172 | ELSE |
|---|
| 173 | PRINT *,'ROS4: Singular factorization at T=',T,'; H=',H |
|---|
| 174 | STOP |
|---|
| 175 | END IF |
|---|
| 176 | END IF |
|---|
| 177 | |
|---|
| 178 | |
|---|
| 179 | C ------------ STAGE 1------------------------- |
|---|
| 180 | DO 50 j = 1,NVAR |
|---|
| 181 | K1(j) = F1(j) |
|---|
| 182 | 50 CONTINUE |
|---|
| 183 | IF (.NOT. Autonomous) THEN |
|---|
| 184 | hgam1 = H*gamma |
|---|
| 185 | DO 60 j=1,NVAR |
|---|
| 186 | K1(j) = K1(j) + hgam1*DFDT(j) |
|---|
| 187 | 60 CONTINUE |
|---|
| 188 | END IF |
|---|
| 189 | CALL KppSolve (JAC, K1) |
|---|
| 190 | |
|---|
| 191 | C ----------- STAGE 2 ------------------------- |
|---|
| 192 | DO 70 j = 1,NVAR |
|---|
| 193 | ynew(j) = y(j) + a21*K1(j) |
|---|
| 194 | 70 CONTINUE |
|---|
| 195 | CALL FUNC_CHEM( T+alpha2*H, ynew, F1) |
|---|
| 196 | nfcn=nfcn+1 |
|---|
| 197 | hc21 = c21/H |
|---|
| 198 | DO 80 j = 1,NVAR |
|---|
| 199 | K2(j) = F1(j) + hc21*K1(j) |
|---|
| 200 | 80 CONTINUE |
|---|
| 201 | IF (.NOT. Autonomous) THEN |
|---|
| 202 | hgam2 = H*gamma2 |
|---|
| 203 | DO 90 j=1,NVAR |
|---|
| 204 | K2(j) = K2(j) + hgam2*DFDT(j) |
|---|
| 205 | 90 CONTINUE |
|---|
| 206 | END IF |
|---|
| 207 | CALL KppSolve (JAC, K2) |
|---|
| 208 | |
|---|
| 209 | |
|---|
| 210 | C ------------ STAGE 3 ------------------------- |
|---|
| 211 | DO 100 j = 1,NVAR |
|---|
| 212 | ynew(j) = y(j) + a31*K1(j) + a32*K2(j) |
|---|
| 213 | 100 CONTINUE |
|---|
| 214 | CALL FUNC_CHEM( T+alpha3*H, ynew, F1) |
|---|
| 215 | nfcn=nfcn+1 |
|---|
| 216 | hc31 = c31/H |
|---|
| 217 | hc32 = c32/H |
|---|
| 218 | DO 110 j = 1,NVAR |
|---|
| 219 | K3(j) = F1(j) + hc31*K1(j) + hc32*K2(j) |
|---|
| 220 | 110 CONTINUE |
|---|
| 221 | IF (.NOT. Autonomous) THEN |
|---|
| 222 | hgam3 = H*gamma3 |
|---|
| 223 | DO 120 j=1,NVAR |
|---|
| 224 | K3(j) = K3(j) + hgam3*DFDT(j) |
|---|
| 225 | 120 CONTINUE |
|---|
| 226 | END IF |
|---|
| 227 | CALL KppSolve (JAC, K3) |
|---|
| 228 | |
|---|
| 229 | C ------------ STAGE 4 ------------------------- |
|---|
| 230 | C Note: uses the same function value as stage 3 |
|---|
| 231 | hc41 = c41/H |
|---|
| 232 | hc42 = c42/H |
|---|
| 233 | hc43 = c43/H |
|---|
| 234 | DO 140 j = 1,NVAR |
|---|
| 235 | K4(j) = F1(j) + hc41*K1(j) + hc42*K2(j) + hc43*K3(j) |
|---|
| 236 | 140 CONTINUE |
|---|
| 237 | IF (.NOT. Autonomous) THEN |
|---|
| 238 | hgam4 = H*gamma4 |
|---|
| 239 | DO 150 j=1,NVAR |
|---|
| 240 | K4(j) = K4(j) + hgam4*DFDT(j) |
|---|
| 241 | 150 CONTINUE |
|---|
| 242 | END IF |
|---|
| 243 | CALL KppSolve (JAC, K4) |
|---|
| 244 | |
|---|
| 245 | |
|---|
| 246 | |
|---|
| 247 | C ---- The Solution --- |
|---|
| 248 | DO 160 j = 1,NVAR |
|---|
| 249 | ynew(j) = y(j) + b1*K1(j) + b2*K2(j) + b3*K3(j) + b4*K4(j) |
|---|
| 250 | 160 CONTINUE |
|---|
| 251 | |
|---|
| 252 | |
|---|
| 253 | C ====== Error estimation ======== |
|---|
| 254 | |
|---|
| 255 | ERR=0.d0 |
|---|
| 256 | DO 170 j = 1,NVAR |
|---|
| 257 | w = AbsTol(j) + RelTol(j)*DMAX1(DABS(y(j)),DABS(ynew(j))) |
|---|
| 258 | e = d1*K1(j) + d2*K2(j) + d3*K3(j) + d4*K4(j) |
|---|
| 259 | ERR = ERR + ( e/w )**2 |
|---|
| 260 | 170 CONTINUE |
|---|
| 261 | ERR = DMAX1( uround, DSQRT( ERR/NVAR ) ) |
|---|
| 262 | |
|---|
| 263 | C ======= Choose the stepsize =============================== |
|---|
| 264 | |
|---|
| 265 | elo = 4.0D0 ! estimator local order |
|---|
| 266 | factor = DMAX1(2.0D-1,DMIN1(6.0D0,ERR**(1.0D0/elo)/.9D0)) |
|---|
| 267 | Hnew = DMIN1(Hmax,DMAX1(Hmin, H/factor)) |
|---|
| 268 | |
|---|
| 269 | C ======= Rejected/Accepted Step ============================ |
|---|
| 270 | |
|---|
| 271 | IF ( (ERR.gt.1).and.(H.gt.Hmin) ) THEN |
|---|
| 272 | IsReject = .true. |
|---|
| 273 | H = DMIN1(H/10,Hnew) |
|---|
| 274 | Nreject = Nreject+1 |
|---|
| 275 | ELSE |
|---|
| 276 | DO 180 i=1,NVAR |
|---|
| 277 | y(i) = ynew(i) |
|---|
| 278 | 180 CONTINUE |
|---|
| 279 | T = Tplus |
|---|
| 280 | IF (.NOT.IsReject) THEN |
|---|
| 281 | H = Hnew ! Do not increase stepsize if previos step was rejected |
|---|
| 282 | END IF |
|---|
| 283 | IsReject = .false. |
|---|
| 284 | Naccept = Naccept+1 |
|---|
| 285 | END IF |
|---|
| 286 | |
|---|
| 287 | C ======= End of the time loop =============================== |
|---|
| 288 | IF ( T .lt. Tnext ) GO TO 10 |
|---|
| 289 | |
|---|
| 290 | C ======= Output Information ================================= |
|---|
| 291 | Info(2) = Nfcn |
|---|
| 292 | Info(3) = Njac |
|---|
| 293 | Info(4) = Naccept |
|---|
| 294 | Info(5) = Nreject |
|---|
| 295 | Hstart = H |
|---|
| 296 | |
|---|
| 297 | RETURN |
|---|
| 298 | END |
|---|
| 299 | |
|---|
| 300 | |
|---|
| 301 | SUBROUTINE FUNC_CHEM( T, Y, P ) |
|---|
| 302 | INCLUDE 'KPP_ROOT_params.h' |
|---|
| 303 | INCLUDE 'KPP_ROOT_global.h' |
|---|
| 304 | KPP_REAL T, Told |
|---|
| 305 | KPP_REAL Y(NVAR), P(NVAR) |
|---|
| 306 | Told = TIME |
|---|
| 307 | TIME = T |
|---|
| 308 | CALL Update_SUN() |
|---|
| 309 | CALL Update_RCONST() |
|---|
| 310 | CALL Fun( Y, FIX, RCONST, P ) |
|---|
| 311 | TIME = Told |
|---|
| 312 | RETURN |
|---|
| 313 | END |
|---|
| 314 | |
|---|
| 315 | |
|---|
| 316 | SUBROUTINE JAC_CHEM( T, Y, J ) |
|---|
| 317 | INCLUDE 'KPP_ROOT_params.h' |
|---|
| 318 | INCLUDE 'KPP_ROOT_global.h' |
|---|
| 319 | KPP_REAL Told, T |
|---|
| 320 | KPP_REAL Y(NVAR), J(LU_NONZERO) |
|---|
| 321 | Told = TIME |
|---|
| 322 | TIME = T |
|---|
| 323 | CALL Update_SUN() |
|---|
| 324 | CALL Update_RCONST() |
|---|
| 325 | CALL Jac_SP( Y, FIX, RCONST, J ) |
|---|
| 326 | TIME = Told |
|---|
| 327 | RETURN |
|---|
| 328 | END |
|---|
| 329 | |
|---|
| 330 | |
|---|
| 331 | |
|---|
| 332 | |
|---|
| 333 | |
|---|
| 334 | |
|---|
