[2696] | 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 | CALL RODAS3(NVAR,TIN,TOUT,STEPMIN,STEPMAX,STEPMIN, |
---|
| 17 | + VAR,ATOL,RTOL, |
---|
| 18 | + Info,FUNC_CHEM,JAC_CHEM) |
---|
| 19 | |
---|
| 20 | RETURN |
---|
| 21 | END |
---|
| 22 | |
---|
| 23 | |
---|
| 24 | |
---|
| 25 | SUBROUTINE RODAS3(N,T,Tnext,Hmin,Hmax,Hstart, |
---|
| 26 | + y,AbsTol,RelTol, |
---|
| 27 | + Info,FUNC_CHEM,JAC_CHEM) |
---|
| 28 | IMPLICIT NONE |
---|
| 29 | INCLUDE 'KPP_ROOT_params.h' |
---|
| 30 | INCLUDE 'KPP_ROOT_sparse.h' |
---|
| 31 | |
---|
| 32 | C Stiffly accurate Rosenbrock 3(2), with |
---|
| 33 | C stiffly accurate embedded formula for error control. |
---|
| 34 | C |
---|
| 35 | C All the arguments aggree with the KPP syntax. |
---|
| 36 | C |
---|
| 37 | C INPUT ARGUMENTS: |
---|
| 38 | C y = Vector of (NVAR) concentrations, contains the |
---|
| 39 | C initial values on input |
---|
| 40 | C [T, Tnext] = the integration interval |
---|
| 41 | C Hmin, Hmax = lower and upper bounds for the selected step-size. |
---|
| 42 | C Note that for Step = Hmin the current computed |
---|
| 43 | C solution is unconditionally accepted by the error |
---|
| 44 | C control mechanism. |
---|
| 45 | C AbsTol, RelTol = (NVAR) dimensional vectors of |
---|
| 46 | C componentwise absolute and relative tolerances. |
---|
| 47 | C FUNC_CHEM = name of routine of derivatives. KPP syntax. |
---|
| 48 | C See the header below. |
---|
| 49 | C JAC_CHEM = name of routine that computes the Jacobian, in |
---|
| 50 | C sparse format. KPP syntax. See the header below. |
---|
| 51 | C Info(1) = 1 for autonomous system |
---|
| 52 | C = 0 for nonautonomous system |
---|
| 53 | C |
---|
| 54 | C OUTPUT ARGUMENTS: |
---|
| 55 | C y = the values of concentrations at Tend. |
---|
| 56 | C T = equals Tend on output. |
---|
| 57 | C Info(2) = # of FUNC_CHEM calls. |
---|
| 58 | C Info(3) = # of JAC_CHEM calls. |
---|
| 59 | C Info(4) = # of accepted steps. |
---|
| 60 | C Info(5) = # of rejected steps. |
---|
| 61 | C |
---|
| 62 | C Adrian Sandu, March 1996 |
---|
| 63 | C The Center for Global and Regional Environmental Research |
---|
| 64 | |
---|
| 65 | KPP_REAL K1(NVAR), K2(NVAR), K3(NVAR), K4(NVAR) |
---|
| 66 | KPP_REAL F1(NVAR), JAC(LU_NONZERO) |
---|
| 67 | KPP_REAL Hmin,Hmax,Hnew,Hstart,ghinv,uround |
---|
| 68 | KPP_REAL y(NVAR), ynew(NVAR) |
---|
| 69 | KPP_REAL AbsTol(NVAR), RelTol(NVAR) |
---|
| 70 | KPP_REAL T, Tnext, H, Hold, Tplus |
---|
| 71 | KPP_REAL ERR, factor, facmax |
---|
| 72 | KPP_REAL c43, tau, x1, x2, ytol, elo |
---|
| 73 | |
---|
| 74 | INTEGER n,nfcn,njac,Naccept,Nreject,i,j,ier |
---|
| 75 | INTEGER Info(5) |
---|
| 76 | LOGICAL IsReject,Autonomous |
---|
| 77 | EXTERNAL FUNC_CHEM, JAC_CHEM |
---|
| 78 | |
---|
| 79 | c Initialization of counters, etc. |
---|
| 80 | Autonomous = Info(1) .EQ. 1 |
---|
| 81 | uround = 1.d-15 |
---|
| 82 | c43 = - 8.d0/3.d0 |
---|
| 83 | H = DMAX1(1.d-8, Hstart) |
---|
| 84 | Hmin = DMAX1(Hmin,uround*(Tnext-T)) |
---|
| 85 | Hmax = DMIN1(Hmax,Tnext-T) |
---|
| 86 | Tplus = T |
---|
| 87 | IsReject = .false. |
---|
| 88 | Naccept = 0 |
---|
| 89 | Nreject = 0 |
---|
| 90 | Nfcn = 0 |
---|
| 91 | Njac = 0 |
---|
| 92 | |
---|
| 93 | |
---|
| 94 | C === Starting the time loop === |
---|
| 95 | 10 continue |
---|
| 96 | Tplus = T + H |
---|
| 97 | if ( Tplus .gt. Tnext ) then |
---|
| 98 | H = Tnext - T |
---|
| 99 | Tplus = Tnext |
---|
| 100 | end if |
---|
| 101 | |
---|
| 102 | CALL JAC_CHEM(NVAR, T, y, JAC) |
---|
| 103 | Njac = Njac+1 |
---|
| 104 | gHinv = -2.0d0/H |
---|
| 105 | do 20 j=1,NVAR |
---|
| 106 | JAC(LU_DIAG(j)) = JAC(LU_DIAG(j)) + gHinv |
---|
| 107 | 20 continue |
---|
| 108 | CALL KppDecomp (JAC, ier) |
---|
| 109 | |
---|
| 110 | if (ier.ne.0) then |
---|
| 111 | if ( H.gt.Hmin) then |
---|
| 112 | H = 5.0d-1*H |
---|
| 113 | go to 10 |
---|
| 114 | else |
---|
| 115 | print *,'IER <> 0, H=',H |
---|
| 116 | stop |
---|
| 117 | end if |
---|
| 118 | end if |
---|
| 119 | |
---|
| 120 | CALL FUNC_CHEM(NVAR, T, y, F1) |
---|
| 121 | |
---|
| 122 | C ====== NONAUTONOMOUS CASE =============== |
---|
| 123 | IF (.not. Autonomous) THEN |
---|
| 124 | tau = DSQRT( uround*DMAX1( 1.0d-5, DABS(T) ) ) |
---|
| 125 | CALL FUNC_CHEM(NVAR, T+tau, y, K2) |
---|
| 126 | nfcn=nfcn+1 |
---|
| 127 | do 30 j = 1,NVAR |
---|
| 128 | K3(j) = ( K2(j)-F1(j) )/tau |
---|
| 129 | 30 continue |
---|
| 130 | |
---|
| 131 | C ----- STAGE 1 (NONAUTONOMOUS) ----- |
---|
| 132 | x1 = 0.5*H |
---|
| 133 | do 40 j = 1,NVAR |
---|
| 134 | K1(j) = F1(j) + x1*K3(j) |
---|
| 135 | 40 continue |
---|
| 136 | CALL KppSolve (JAC, K1) |
---|
| 137 | |
---|
| 138 | C ----- STAGE 2 (NONAUTONOMOUS) ----- |
---|
| 139 | x1 = 4.d0/H |
---|
| 140 | x2 = 1.5d0*H |
---|
| 141 | do 50 j = 1,NVAR |
---|
| 142 | K2(j) = F1(j) - x1*K1(j) + x2*K3(j) |
---|
| 143 | 50 continue |
---|
| 144 | CALL KppSolve (JAC, K2) |
---|
| 145 | |
---|
| 146 | C ====== AUTONOMOUS CASE =============== |
---|
| 147 | ELSE |
---|
| 148 | C ----- STAGE 1 (AUTONOMOUS) ----- |
---|
| 149 | do 60 j = 1,NVAR |
---|
| 150 | K1(j) = F1(j) |
---|
| 151 | 60 continue |
---|
| 152 | CALL KppSolve (JAC, K1) |
---|
| 153 | |
---|
| 154 | C ----- STAGE 2 (AUTONOMOUS) ----- |
---|
| 155 | x1 = 4.d0/H |
---|
| 156 | do 70 j = 1,NVAR |
---|
| 157 | K2(j) = F1(j) - x1*K1(j) |
---|
| 158 | 70 continue |
---|
| 159 | CALL KppSolve (JAC, K2) |
---|
| 160 | END IF |
---|
| 161 | |
---|
| 162 | C ----- STAGE 3 ----- |
---|
| 163 | do 80 j = 1,NVAR |
---|
| 164 | ynew(j) = y(j) - 2.0d0*K1(j) |
---|
| 165 | 80 continue |
---|
| 166 | CALL FUNC_CHEM(NVAR, T+H, ynew, F1) |
---|
| 167 | nfcn=nfcn+1 |
---|
| 168 | do 90 j = 1,NVAR |
---|
| 169 | K3(j) = F1(j) + ( -K1(j) + K2(j) )/H |
---|
| 170 | 90 continue |
---|
| 171 | CALL KppSolve (JAC, K3) |
---|
| 172 | |
---|
| 173 | C ----- STAGE 4 ----- |
---|
| 174 | do 100 j = 1,NVAR |
---|
| 175 | ynew(j) = y(j) - 2.0d0*K1(j) - K3(j) |
---|
| 176 | 100 continue |
---|
| 177 | CALL FUNC_CHEM(NVAR, T+H, ynew, F1) |
---|
| 178 | nfcn=nfcn+1 |
---|
| 179 | do 110 j = 1,NVAR |
---|
| 180 | K4(j) = F1(j) + ( -K1(j) + K2(j) - C43*K3(j) )/H |
---|
| 181 | 110 continue |
---|
| 182 | CALL KppSolve (JAC, K4) |
---|
| 183 | |
---|
| 184 | C ---- The Solution --- |
---|
| 185 | |
---|
| 186 | do 120 j = 1,NVAR |
---|
| 187 | ynew(j) = y(j) - 2.0d0*K1(j) - K3(j) - K4(j) |
---|
| 188 | 120 continue |
---|
| 189 | |
---|
| 190 | |
---|
| 191 | C ====== Error estimation ======== |
---|
| 192 | |
---|
| 193 | ERR=0.d0 |
---|
| 194 | do 130 i=1,NVAR |
---|
| 195 | ytol = AbsTol(i) + RelTol(i)*DABS(ynew(i)) |
---|
| 196 | ERR = ERR + ( K4(i)/ytol )**2 |
---|
| 197 | 130 continue |
---|
| 198 | ERR = DMAX1( uround, DSQRT( ERR/NVAR ) ) |
---|
| 199 | |
---|
| 200 | C ======= Choose the stepsize =============================== |
---|
| 201 | elo = 3.0D0 ! estimator local order |
---|
| 202 | factor = DMAX1(2.0D-1,DMIN1(6.0D0,ERR**(1.0D0/elo)/.9D0)) |
---|
| 203 | Hnew = DMIN1(Hmax,DMAX1(Hmin, H/factor)) |
---|
| 204 | |
---|
| 205 | C ======= Rejected/Accepted Step ============================ |
---|
| 206 | |
---|
| 207 | IF ( (ERR.gt.1).and.(H.gt.Hmin) ) THEN |
---|
| 208 | IsReject = .true. |
---|
| 209 | H = DMIN1(H/10,Hnew) |
---|
| 210 | Nreject = Nreject+1 |
---|
| 211 | ELSE |
---|
| 212 | DO 140 i=1,NVAR |
---|
| 213 | y(i) = ynew(i) |
---|
| 214 | 140 CONTINUE |
---|
| 215 | T = Tplus |
---|
| 216 | IF (.NOT.IsReject) THEN |
---|
| 217 | H = Hnew ! Do not increase stepsize if previos step was rejected |
---|
| 218 | END IF |
---|
| 219 | IsReject = .false. |
---|
| 220 | Naccept = Naccept+1 |
---|
| 221 | END IF |
---|
| 222 | |
---|
| 223 | |
---|
| 224 | C ======= End of the time loop =============================== |
---|
| 225 | if ( T .lt. Tnext ) go to 10 |
---|
| 226 | |
---|
| 227 | |
---|
| 228 | |
---|
| 229 | C ======= Output Information ================================= |
---|
| 230 | Info(2) = Nfcn |
---|
| 231 | Info(3) = Njac |
---|
| 232 | Info(4) = Naccept |
---|
| 233 | Info(5) = Nreject |
---|
| 234 | Hstart = H |
---|
| 235 | |
---|
| 236 | RETURN |
---|
| 237 | END |
---|
| 238 | |
---|
| 239 | |
---|
| 240 | |
---|
| 241 | SUBROUTINE FUNC_CHEM(N, T, Y, P) |
---|
| 242 | INCLUDE 'KPP_ROOT_params.h' |
---|
| 243 | INCLUDE 'KPP_ROOT_global.h' |
---|
| 244 | INTEGER N |
---|
| 245 | KPP_REAL T, Told |
---|
| 246 | KPP_REAL Y(NVAR), P(NVAR) |
---|
| 247 | Told = TIME |
---|
| 248 | TIME = T |
---|
| 249 | CALL Update_SUN() |
---|
| 250 | CALL Update_RCONST() |
---|
| 251 | CALL Fun( Y, FIX, RCONST, P ) |
---|
| 252 | TIME = Told |
---|
| 253 | RETURN |
---|
| 254 | END |
---|
| 255 | |
---|
| 256 | |
---|
| 257 | SUBROUTINE JAC_CHEM(N, T, Y, J) |
---|
| 258 | INCLUDE 'KPP_ROOT_params.h' |
---|
| 259 | INCLUDE 'KPP_ROOT_global.h' |
---|
| 260 | INTEGER N |
---|
| 261 | KPP_REAL Told, T |
---|
| 262 | KPP_REAL Y(NVAR), J(LU_NONZERO) |
---|
| 263 | Told = TIME |
---|
| 264 | TIME = T |
---|
| 265 | CALL Update_SUN() |
---|
| 266 | CALL Update_RCONST() |
---|
| 267 | CALL Jac_SP( Y, FIX, RCONST, J ) |
---|
| 268 | TIME = Told |
---|
| 269 | RETURN |
---|
| 270 | END |
---|
| 271 | |
---|