| 1 | SUBROUTINE INTEGRATE( NSENSIT, Y, 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 | C Y - Concentrations and Sensitivities |
|---|
| 11 | KPP_REAL Y(NVAR*(NSENSIT+1)) |
|---|
| 12 | C --- Note: Y contains: (1:NVAR) concentrations, followed by |
|---|
| 13 | C --- (1:NVAR) sensitivities w.r.t. first parameter, followed by |
|---|
| 14 | C --- etc., followed by |
|---|
| 15 | C --- (1:NVAR) sensitivities w.r.t. NSENSIT's parameter |
|---|
| 16 | |
|---|
| 17 | INTEGER INFO(5) |
|---|
| 18 | |
|---|
| 19 | EXTERNAL FUNC_CHEM, JAC_CHEM |
|---|
| 20 | |
|---|
| 21 | INFO(1) = Autonomous |
|---|
| 22 | |
|---|
| 23 | CALL ROS3_DDM(NVAR,NSENSIT,TIN,TOUT,STEPMIN,STEPMAX, |
|---|
| 24 | + STEPMIN,Y,ATOL,RTOL, |
|---|
| 25 | + Info,FUNC_CHEM,JAC_CHEM) |
|---|
| 26 | |
|---|
| 27 | RETURN |
|---|
| 28 | END |
|---|
| 29 | |
|---|
| 30 | |
|---|
| 31 | SUBROUTINE ROS3_DDM(N,NSENSIT,T,Tnext,Hmin,Hmax,Hstart, |
|---|
| 32 | + y,AbsTol,RelTol, |
|---|
| 33 | + Info,FUNC_CHEM,JAC_CHEM) |
|---|
| 34 | IMPLICIT NONE |
|---|
| 35 | INCLUDE 'KPP_ROOT_params.h' |
|---|
| 36 | INCLUDE 'KPP_ROOT_global.h' |
|---|
| 37 | INCLUDE 'KPP_ROOT_sparse.h' |
|---|
| 38 | |
|---|
| 39 | C L-stable Rosenbrock 3(2), with |
|---|
| 40 | C strongly A-stable embedded formula for error control. |
|---|
| 41 | C |
|---|
| 42 | C Direct decoupled computation of sensitivities. |
|---|
| 43 | C The global variable DDMTYPE distinguishes between: |
|---|
| 44 | C DDMTYPE = 0 : sensitivities w.r.t. initial values |
|---|
| 45 | C DDMTYPE = 1 : sensitivities w.r.t. parameters |
|---|
| 46 | C |
|---|
| 47 | C All the arguments aggree with the KPP syntax. |
|---|
| 48 | C |
|---|
| 49 | C INPUT ARGUMENTS: |
|---|
| 50 | C y = Vector of: (1:NVAR) concentrations, followed by |
|---|
| 51 | C (1:NVAR) sensitivities w.r.t. first parameter, followed by |
|---|
| 52 | C etc., followed by |
|---|
| 53 | C (1:NVAR) sensitivities w.r.t. NSENSIT's parameter |
|---|
| 54 | C (y contains initial values at input, final values at output) |
|---|
| 55 | C [T, Tnext] = the integration interval |
|---|
| 56 | C Hmin, Hmax = lower and upper bounds for the selected step-size. |
|---|
| 57 | C Note that for Step = Hmin the current computed |
|---|
| 58 | C solution is unconditionally accepted by the error |
|---|
| 59 | C control mechanism. |
|---|
| 60 | C AbsTol, RelTol = (NVAR) dimensional vectors of |
|---|
| 61 | C componentwise absolute and relative tolerances. |
|---|
| 62 | C FUNC_CHEM = name of routine of derivatives. KPP syntax. |
|---|
| 63 | C See the header below. |
|---|
| 64 | C JAC_CHEM = name of routine that computes the Jacobian, in |
|---|
| 65 | C sparse format. KPP syntax. See the header below. |
|---|
| 66 | C Info(1) = 1 for Autonomous system |
|---|
| 67 | C = 0 for nonAutonomous system |
|---|
| 68 | C |
|---|
| 69 | C OUTPUT ARGUMENTS: |
|---|
| 70 | C y = the values of concentrations at TEND. |
|---|
| 71 | C T = equals TEND on output. |
|---|
| 72 | C Info(2) = # of FUNC_CHEM calls. |
|---|
| 73 | C Info(3) = # of JAC_CHEM calls. |
|---|
| 74 | C Info(4) = # of accepted steps. |
|---|
| 75 | C Info(5) = # of rejected steps. |
|---|
| 76 | C |
|---|
| 77 | C Adrian Sandu, April 1996 |
|---|
| 78 | C The Center for Global and Regional Environmental Research |
|---|
| 79 | INTEGER NSENSIT |
|---|
| 80 | KPP_REAL y(NVAR*(NSENSIT+1)), ynew(NVAR*(NSENSIT+1)) |
|---|
| 81 | KPP_REAL K1(NVAR*(NSENSIT+1)) |
|---|
| 82 | KPP_REAL K2(NVAR*(NSENSIT+1)) |
|---|
| 83 | KPP_REAL K3(NVAR*(NSENSIT+1)) |
|---|
| 84 | KPP_REAL DFDT(NVAR*(NSENSIT+1)) |
|---|
| 85 | KPP_REAL DFDP(NVAR*NSENSIT), DFDPDT(NVAR*NSENSIT) |
|---|
| 86 | KPP_REAL DJDP(NVAR*NSENSIT) |
|---|
| 87 | KPP_REAL JAC(LU_NONZERO), AJAC(LU_NONZERO) |
|---|
| 88 | KPP_REAL DJDT(LU_NONZERO) |
|---|
| 89 | KPP_REAL Fv(NVAR), Hv(NVAR) |
|---|
| 90 | KPP_REAL HESS(NHESS) |
|---|
| 91 | KPP_REAL Hmin,Hmax,Hstart,ghinv,uround |
|---|
| 92 | KPP_REAL AbsTol(NVAR), RelTol(NVAR) |
|---|
| 93 | KPP_REAL T, Tnext, Tplus, H, Hnew, elo |
|---|
| 94 | KPP_REAL ERR, factor, facmax |
|---|
| 95 | INTEGER n,nfcn,njac,Naccept,Nreject,i,j,ier |
|---|
| 96 | INTEGER Info(5) |
|---|
| 97 | LOGICAL IsReject,Autonomous |
|---|
| 98 | EXTERNAL FUNC_CHEM, JAC_CHEM, HESS_CHEM |
|---|
| 99 | |
|---|
| 100 | KPP_REAL gamma, c21, c31,c32,b1,b2,b3,d1,d2,d3,a21,a31,a32 |
|---|
| 101 | KPP_REAL alpha2, alpha3, g1, g2, g3, x1, x2, x3, ytol |
|---|
| 102 | KPP_REAL dround, tau |
|---|
| 103 | |
|---|
| 104 | gamma= .43586652150845899941601945119356d+00 |
|---|
| 105 | c21= -.10156171083877702091975600115545d+01 |
|---|
| 106 | c31= .40759956452537699824805835358067d+01 |
|---|
| 107 | c32= .92076794298330791242156818474003d+01 |
|---|
| 108 | b1= .10000000000000000000000000000000d+01 |
|---|
| 109 | b2= .61697947043828245592553615689730d+01 |
|---|
| 110 | b3= -.42772256543218573326238373806514d+00 |
|---|
| 111 | d1= .50000000000000000000000000000000d+00 |
|---|
| 112 | d2= -.29079558716805469821718236208017d+01 |
|---|
| 113 | d3= .22354069897811569627360909276199d+00 |
|---|
| 114 | a21 = 1.d0 |
|---|
| 115 | a31 = 1.d0 |
|---|
| 116 | a32 = 0.d0 |
|---|
| 117 | alpha2 = gamma |
|---|
| 118 | g1= .43586652150845899941601945119356d+00 |
|---|
| 119 | g2= .24291996454816804366592249683314d+00 |
|---|
| 120 | g3= .21851380027664058511513169485832d+01 |
|---|
| 121 | |
|---|
| 122 | |
|---|
| 123 | c Initialization of counters, etc. |
|---|
| 124 | Autonomous = Info(1) .EQ. 1 |
|---|
| 125 | uround = 1.d-15 |
|---|
| 126 | dround = DSQRT(uround) |
|---|
| 127 | IF (Hmax.le.0.D0) THEN |
|---|
| 128 | Hmax = DABS(Tnext-T) |
|---|
| 129 | END IF |
|---|
| 130 | H = DMAX1(1.d-8, Hstart) |
|---|
| 131 | Tplus = T |
|---|
| 132 | IsReject = .false. |
|---|
| 133 | Naccept = 0 |
|---|
| 134 | Nreject = 0 |
|---|
| 135 | Nfcn = 0 |
|---|
| 136 | Njac = 0 |
|---|
| 137 | |
|---|
| 138 | |
|---|
| 139 | C === Starting the time loop === |
|---|
| 140 | 10 CONTINUE |
|---|
| 141 | |
|---|
| 142 | C ====== Initial Function, Jacobian, and Hessian values =============== |
|---|
| 143 | CALL FUNC_CHEM(NVAR, T, y, Fv) |
|---|
| 144 | Nfcn = Nfcn + 1 |
|---|
| 145 | CALL JAC_CHEM(NVAR, T, y, JAC) |
|---|
| 146 | Njac = Njac + 1 |
|---|
| 147 | CALL HESS_CHEM( NVAR, T, y, HESS ) |
|---|
| 148 | IF (DDMTYPE .EQ. 1) THEN |
|---|
| 149 | CALL DFUNDPAR(NVAR, NSENSIT, T, y, DFDP) |
|---|
| 150 | END IF |
|---|
| 151 | |
|---|
| 152 | C ====== Time derivatives for NONAutonomousous CASE =============== |
|---|
| 153 | IF (.not. Autonomous) THEN |
|---|
| 154 | tau = DSIGN(dround*DMAX1( 1.0d-6, DABS(T) ), T) |
|---|
| 155 | CALL FUNC_CHEM(NVAR, T+tau, y, K2) |
|---|
| 156 | nfcn=nfcn+1 |
|---|
| 157 | DO 20 j = 1,NVAR |
|---|
| 158 | DFDT(j) = ( K2(j)-Fv(j) )/tau |
|---|
| 159 | 20 CONTINUE |
|---|
| 160 | CALL JAC_CHEM(NVAR, T+tau, y, AJAC) |
|---|
| 161 | DO 30 j = 1,LU_NONZERO |
|---|
| 162 | DJDT(j) = ( AJAC(j)-JAC(j) )/tau |
|---|
| 163 | 30 CONTINUE |
|---|
| 164 | DO 40 i=1,NSENSIT |
|---|
| 165 | CALL Jac_SP_Vec (DJDT,y(i*NVAR+1),DFDT(i*NVAR+1)) |
|---|
| 166 | 40 CONTINUE |
|---|
| 167 | END IF |
|---|
| 168 | |
|---|
| 169 | |
|---|
| 170 | Tplus = T + H |
|---|
| 171 | IF ( Tplus .gt. Tnext ) then |
|---|
| 172 | H = Tnext - T |
|---|
| 173 | Tplus = Tnext |
|---|
| 174 | END IF |
|---|
| 175 | |
|---|
| 176 | gHinv = 1.0d0/(gamma*H) |
|---|
| 177 | DO 50 j=1,LU_NONZERO |
|---|
| 178 | AJAC(j) = -JAC(j) |
|---|
| 179 | 50 CONTINUE |
|---|
| 180 | DO 60 j=1,NVAR |
|---|
| 181 | AJAC(LU_DIAG(j)) = AJAC(LU_DIAG(j)) + gHinv |
|---|
| 182 | 60 CONTINUE |
|---|
| 183 | CALL KppDecomp (AJAC, ier) |
|---|
| 184 | |
|---|
| 185 | IF (ier.NE.0) THEN |
|---|
| 186 | IF ( H.GT.Hmin) THEN |
|---|
| 187 | H = 5.0d-1*H |
|---|
| 188 | GO TO 10 |
|---|
| 189 | ELSE |
|---|
| 190 | PRINT *,'IER <> 0, H=',H |
|---|
| 191 | STOP |
|---|
| 192 | END IF |
|---|
| 193 | END IF |
|---|
| 194 | |
|---|
| 195 | Autonomous = .true. |
|---|
| 196 | |
|---|
| 197 | C ------------------------------- STAGE 1 -------------------------------------- |
|---|
| 198 | DO 70 j = 1,NVAR |
|---|
| 199 | K1(j) = Fv(j) |
|---|
| 200 | 70 CONTINUE |
|---|
| 201 | IF (.NOT.Autonomous) THEN |
|---|
| 202 | x1 = gamma*H |
|---|
| 203 | DO 80 j = 1,NVAR |
|---|
| 204 | K1(j) = K1(j) + x1*DFDT(j) |
|---|
| 205 | 80 CONTINUE |
|---|
| 206 | END IF |
|---|
| 207 | CALL KppSolve (AJAC, K1) |
|---|
| 208 | C --- If derivative w.r.t. parameters |
|---|
| 209 | IF (DDMTYPE .EQ. 1) THEN |
|---|
| 210 | CALL DJACDPAR(NVAR, NSENSIT, T, y, K1(1), DJDP) |
|---|
| 211 | END IF |
|---|
| 212 | C --- End of derivative w.r.t. parameters |
|---|
| 213 | |
|---|
| 214 | DO 110 i=1,NSENSIT |
|---|
| 215 | CALL Jac_SP_Vec (JAC,y(i*NVAR+1),K1(i*NVAR+1)) |
|---|
| 216 | CALL Hess_Vec ( HESS, y(i*NVAR+1), K1(1), Hv ) |
|---|
| 217 | DO 90 j=1,NVAR |
|---|
| 218 | K1(i*NVAR+j) = K1(i*NVAR+j) + Hv(j) |
|---|
| 219 | 90 CONTINUE |
|---|
| 220 | IF (.NOT. Autonomous) THEN |
|---|
| 221 | DO 100 j=1,NVAR |
|---|
| 222 | K1(i*NVAR+j) = K1(i*NVAR+j) + x1*DFDT(i*NVAR+j) |
|---|
| 223 | 100 CONTINUE |
|---|
| 224 | END IF |
|---|
| 225 | C --- If derivative w.r.t. parameters |
|---|
| 226 | IF (DDMTYPE .EQ. 1) THEN |
|---|
| 227 | DO 44 j = 1,NVAR |
|---|
| 228 | K1(i*NVAR+j) = K1(i*NVAR+j) + DFDP((i-1)*NVAR+j) |
|---|
| 229 | & + DJDP((i-1)*NVAR+j) |
|---|
| 230 | 44 CONTINUE |
|---|
| 231 | END IF |
|---|
| 232 | C --- End of derivative w.r.t. parameters |
|---|
| 233 | CALL KppSolve (AJAC, K1(i*NVAR+1)) |
|---|
| 234 | 110 CONTINUE |
|---|
| 235 | |
|---|
| 236 | C ------------------------------- STAGE 2 -------------------------------------- |
|---|
| 237 | DO 120 j = 1,NVAR*(NSENSIT+1) |
|---|
| 238 | ynew(j) = y(j) + a21*K1(j) |
|---|
| 239 | 120 CONTINUE |
|---|
| 240 | CALL FUNC_CHEM(NVAR, T + alpha2*H, ynew, Fv) |
|---|
| 241 | IF (DDMTYPE .EQ. 1) THEN |
|---|
| 242 | CALL DFUNDPAR(NVAR, NSENSIT, T+alpha3*H, ynew, DFDP) |
|---|
| 243 | END IF |
|---|
| 244 | nfcn=nfcn+1 |
|---|
| 245 | x1 = c21/H |
|---|
| 246 | DO 130 j = 1,NVAR |
|---|
| 247 | K2(j) = Fv(j) + x1*K1(j) |
|---|
| 248 | 130 CONTINUE |
|---|
| 249 | IF (.NOT.Autonomous) THEN |
|---|
| 250 | x2 = g2*H |
|---|
| 251 | DO 140 j = 1,NVAR |
|---|
| 252 | K2(j) = K2(j) + x2*DFDT(j) |
|---|
| 253 | 140 CONTINUE |
|---|
| 254 | END IF |
|---|
| 255 | CALL KppSolve (AJAC, K2) |
|---|
| 256 | C --- If derivative w.r.t. parameters |
|---|
| 257 | IF (DDMTYPE .EQ. 1) THEN |
|---|
| 258 | CALL DJACDPAR(NVAR, NSENSIT, T, y, K2(1), DJDP) |
|---|
| 259 | END IF |
|---|
| 260 | C --- End of derivative w.r.t. parameters |
|---|
| 261 | |
|---|
| 262 | CALL JAC_CHEM(NVAR, T+alpha2*H, ynew, JAC) |
|---|
| 263 | njac=njac+1 |
|---|
| 264 | DO 170 i=1,NSENSIT |
|---|
| 265 | CALL Jac_SP_Vec (JAC,ynew(i*NVAR+1),K2(i*NVAR+1)) |
|---|
| 266 | CALL Hess_Vec ( HESS, y(i*NVAR+1), K2(1), Hv ) |
|---|
| 267 | DO 150 j = 1,NVAR |
|---|
| 268 | K2(i*NVAR+j) = K2(i*NVAR+j) + x1*K1(i*NVAR+j) |
|---|
| 269 | & + Hv(j) |
|---|
| 270 | 150 CONTINUE |
|---|
| 271 | IF (.NOT. Autonomous) THEN |
|---|
| 272 | DO 160 j=1,NVAR |
|---|
| 273 | K2(i*NVAR+j) = K2(i*NVAR+j) + x2*DFDT(i*NVAR+j) |
|---|
| 274 | 160 CONTINUE |
|---|
| 275 | END IF |
|---|
| 276 | C --- If derivative w.r.t. parameters |
|---|
| 277 | IF (DDMTYPE .EQ. 1) THEN |
|---|
| 278 | DO 165 j = 1,NVAR |
|---|
| 279 | K2(i*NVAR+j) = K2(i*NVAR+j) + DFDP((i-1)*NVAR+j) |
|---|
| 280 | & + DJDP((i-1)*NVAR+j) |
|---|
| 281 | 165 CONTINUE |
|---|
| 282 | END IF |
|---|
| 283 | C --- End of derivative w.r.t. parameters |
|---|
| 284 | CALL KppSolve (AJAC, K2(i*NVAR+1)) |
|---|
| 285 | 170 CONTINUE |
|---|
| 286 | |
|---|
| 287 | C ------------------------------- STAGE 3 -------------------------------------- |
|---|
| 288 | x1 = c31/H |
|---|
| 289 | x2 = c32/H |
|---|
| 290 | DO 180 j = 1,NVAR |
|---|
| 291 | K3(j) = Fv(j) + x1*K1(j) + x2*K2(j) |
|---|
| 292 | 180 CONTINUE |
|---|
| 293 | IF (.NOT.Autonomous) THEN |
|---|
| 294 | x3 = g3*H |
|---|
| 295 | DO 190 j = 1,NVAR |
|---|
| 296 | K3(j) = K3(j) + x3*DFDT(j) |
|---|
| 297 | 190 CONTINUE |
|---|
| 298 | END IF |
|---|
| 299 | CALL KppSolve (AJAC, K3) |
|---|
| 300 | C --- If derivative w.r.t. parameters |
|---|
| 301 | IF (DDMTYPE .EQ. 1) THEN |
|---|
| 302 | CALL DJACDPAR(NVAR, NSENSIT, T, y, K3(1), DJDP) |
|---|
| 303 | END IF |
|---|
| 304 | C --- End of derivative w.r.t. parameters |
|---|
| 305 | |
|---|
| 306 | DO 220 i=1,NSENSIT |
|---|
| 307 | CALL Jac_SP_Vec (JAC,ynew(i*NVAR+1),K3(i*NVAR+1)) |
|---|
| 308 | CALL Hess_Vec ( HESS, y(i*NVAR+1), K3(1), Hv ) |
|---|
| 309 | DO 200 j = 1,NVAR |
|---|
| 310 | K3(i*NVAR+j) = K3(i*NVAR+j) +x1*K1(i*NVAR+j) |
|---|
| 311 | & + x2*K2(i*NVAR+j) + Hv(j) |
|---|
| 312 | 200 CONTINUE |
|---|
| 313 | IF (.NOT. Autonomous) THEN |
|---|
| 314 | DO 210 j=1,NVAR |
|---|
| 315 | K3(i*NVAR+j) = K3(i*NVAR+j) + x3*DFDT(i*NVAR+j) |
|---|
| 316 | 210 CONTINUE |
|---|
| 317 | END IF |
|---|
| 318 | C --- If derivative w.r.t. parameters |
|---|
| 319 | IF (DDMTYPE .EQ. 1) THEN |
|---|
| 320 | DO 215 j = 1,NVAR |
|---|
| 321 | K3(i*NVAR+j) = K3(i*NVAR+j) + DFDP((i-1)*NVAR+j) |
|---|
| 322 | & + DJDP((i-1)*NVAR+j) |
|---|
| 323 | 215 CONTINUE |
|---|
| 324 | END IF |
|---|
| 325 | C --- End of derivative w.r.t. parameters |
|---|
| 326 | CALL KppSolve (AJAC, K3(i*NVAR+1)) |
|---|
| 327 | 220 CONTINUE |
|---|
| 328 | |
|---|
| 329 | C ------------------------------ The Solution --- |
|---|
| 330 | |
|---|
| 331 | DO 230 j = 1,NVAR*(NSENSIT+1) |
|---|
| 332 | ynew(j) = y(j) + b1*K1(j) + b2*K2(j) + b3*K3(j) |
|---|
| 333 | 230 CONTINUE |
|---|
| 334 | |
|---|
| 335 | |
|---|
| 336 | C ====== Error estimation ======== |
|---|
| 337 | |
|---|
| 338 | ERR=0.d0 |
|---|
| 339 | DO 240 i=1,NVAR |
|---|
| 340 | ytol = AbsTol(i) + RelTol(i)*DABS(ynew(i)) |
|---|
| 341 | ERR=ERR+((d1*K1(i)+d2*K2(i)+d3*K3(i))/ytol)**2 |
|---|
| 342 | 240 CONTINUE |
|---|
| 343 | ERR = DMAX1( uround, DSQRT( ERR/NVAR ) ) |
|---|
| 344 | |
|---|
| 345 | C ======= Choose the stepsize =============================== |
|---|
| 346 | |
|---|
| 347 | elo = 3.0D0 ! estimator local order |
|---|
| 348 | factor = DMAX1(2.0D-1,DMIN1(6.0D0,ERR**(1.0D0/elo)/.9D0)) |
|---|
| 349 | Hnew = DMIN1(Hmax,DMAX1(Hmin, H/factor)) |
|---|
| 350 | |
|---|
| 351 | C ======= Rejected/Accepted Step ============================ |
|---|
| 352 | |
|---|
| 353 | IF ( (ERR.gt.1).and.(H.gt.Hmin) ) THEN |
|---|
| 354 | IsReject = .true. |
|---|
| 355 | H = DMIN1(H/10,Hnew) |
|---|
| 356 | Nreject = Nreject+1 |
|---|
| 357 | GO TO 10 |
|---|
| 358 | ELSE |
|---|
| 359 | DO 250 j=1,NVAR*(NSENSIT+1) |
|---|
| 360 | y(j) = ynew(j) |
|---|
| 361 | 250 CONTINUE |
|---|
| 362 | T = Tplus |
|---|
| 363 | IF (.NOT.IsReject) THEN |
|---|
| 364 | H = Hnew ! Do not increase stepsize IF previos step was rejected |
|---|
| 365 | END IF |
|---|
| 366 | IsReject = .false. |
|---|
| 367 | Naccept = Naccept+1 |
|---|
| 368 | END IF |
|---|
| 369 | |
|---|
| 370 | |
|---|
| 371 | C ======= End of the time loop =============================== |
|---|
| 372 | IF ( T .lt. Tnext ) GO TO 10 |
|---|
| 373 | |
|---|
| 374 | |
|---|
| 375 | C ======= Output Information ================================= |
|---|
| 376 | Info(2) = Nfcn |
|---|
| 377 | Info(3) = Njac |
|---|
| 378 | Info(4) = Naccept |
|---|
| 379 | Info(5) = Nreject |
|---|
| 380 | Hstart = H |
|---|
| 381 | |
|---|
| 382 | RETURN |
|---|
| 383 | END |
|---|
| 384 | |
|---|
| 385 | |
|---|
| 386 | |
|---|
| 387 | SUBROUTINE FUNC_CHEM(N, T, Y, P) |
|---|
| 388 | INCLUDE 'KPP_ROOT_params.h' |
|---|
| 389 | INCLUDE 'KPP_ROOT_global.h' |
|---|
| 390 | INTEGER N |
|---|
| 391 | KPP_REAL T, Told |
|---|
| 392 | KPP_REAL Y(NVAR), P(NVAR) |
|---|
| 393 | Told = TIME |
|---|
| 394 | TIME = T |
|---|
| 395 | CALL Update_SUN() |
|---|
| 396 | CALL Update_RCONST() |
|---|
| 397 | CALL Fun( Y, FIX, RCONST, P ) |
|---|
| 398 | TIME = Told |
|---|
| 399 | RETURN |
|---|
| 400 | END |
|---|
| 401 | |
|---|
| 402 | |
|---|
| 403 | SUBROUTINE DFUNDPAR(N, NSENSIT, T, Y, P) |
|---|
| 404 | C --- Computes the partial derivatives of FUNC_CHEM w.r.t. parameters |
|---|
| 405 | INCLUDE 'KPP_ROOT_params.h' |
|---|
| 406 | INCLUDE 'KPP_ROOT_global.h' |
|---|
| 407 | C --- NCOEFF, JCOEFF useful for derivatives w.r.t. rate coefficients |
|---|
| 408 | INTEGER N |
|---|
| 409 | INTEGER NCOEFF, JCOEFF(NREACT) |
|---|
| 410 | COMMON /DDMRCOEFF/ NCOEFF, JCOEFF |
|---|
| 411 | |
|---|
| 412 | KPP_REAL T, Told |
|---|
| 413 | KPP_REAL Y(NVAR), P(NVAR*NSENSIT) |
|---|
| 414 | Told = TIME |
|---|
| 415 | TIME = T |
|---|
| 416 | CALL Update_SUN() |
|---|
| 417 | CALL Update_RCONST() |
|---|
| 418 | C |
|---|
| 419 | IF (DDMTYPE .EQ. 0) THEN |
|---|
| 420 | C --- Note: the values below are for sensitivities w.r.t. initial values; |
|---|
| 421 | C --- they may have to be changed for other applications |
|---|
| 422 | DO j=1,NSENSIT |
|---|
| 423 | DO i=1,NVAR |
|---|
| 424 | P(i+NVAR*(j-1)) = 0.0D0 |
|---|
| 425 | END DO |
|---|
| 426 | END DO |
|---|
| 427 | ELSE |
|---|
| 428 | C --- Example: the call below is for sensitivities w.r.t. rate coefficients; |
|---|
| 429 | C --- JCOEFF(1:NSENSIT) are the indices of the NSENSIT rate coefficients |
|---|
| 430 | C --- w.r.t. which one differentiates |
|---|
| 431 | CALL dFun_dRcoeff( Y, FIX, NCOEFF, JCOEFF, P ) |
|---|
| 432 | END IF |
|---|
| 433 | TIME = Told |
|---|
| 434 | RETURN |
|---|
| 435 | END |
|---|
| 436 | |
|---|
| 437 | SUBROUTINE JAC_CHEM(N, T, Y, J) |
|---|
| 438 | INCLUDE 'KPP_ROOT_params.h' |
|---|
| 439 | INCLUDE 'KPP_ROOT_global.h' |
|---|
| 440 | INTEGER N |
|---|
| 441 | KPP_REAL Told, T |
|---|
| 442 | KPP_REAL Y(NVAR), J(LU_NONZERO) |
|---|
| 443 | Told = TIME |
|---|
| 444 | TIME = T |
|---|
| 445 | CALL Update_SUN() |
|---|
| 446 | CALL Update_RCONST() |
|---|
| 447 | CALL Jac_SP( Y, FIX, RCONST, J ) |
|---|
| 448 | TIME = Told |
|---|
| 449 | RETURN |
|---|
| 450 | END |
|---|
| 451 | |
|---|
| 452 | |
|---|
| 453 | SUBROUTINE DJACDPAR(N, NSENSIT, T, Y, U, P) |
|---|
| 454 | C --- Computes the partial derivatives of JAC w.r.t. parameters times user vector U |
|---|
| 455 | INCLUDE 'KPP_ROOT_params.h' |
|---|
| 456 | INCLUDE 'KPP_ROOT_global.h' |
|---|
| 457 | C --- NCOEFF, JCOEFF useful for derivatives w.r.t. rate coefficients |
|---|
| 458 | INTEGER N |
|---|
| 459 | INTEGER NCOEFF, JCOEFF(NREACT) |
|---|
| 460 | COMMON /DDMRCOEFF/ NCOEFF, JCOEFF |
|---|
| 461 | |
|---|
| 462 | KPP_REAL T, Told |
|---|
| 463 | KPP_REAL Y(NVAR), U(NVAR) |
|---|
| 464 | KPP_REAL P(NVAR*NSENSIT) |
|---|
| 465 | Told = TIME |
|---|
| 466 | TIME = T |
|---|
| 467 | CALL Update_SUN() |
|---|
| 468 | CALL Update_RCONST() |
|---|
| 469 | C |
|---|
| 470 | IF (DDMTYPE .EQ. 0) THEN |
|---|
| 471 | C --- Note: the values below are for sensitivities w.r.t. initial values; |
|---|
| 472 | C --- they may have to be changed for other applications |
|---|
| 473 | DO j=1,NSENSIT |
|---|
| 474 | DO i=1,NVAR |
|---|
| 475 | P(i+NVAR*(j-1)) = 0.0D0 |
|---|
| 476 | END DO |
|---|
| 477 | END DO |
|---|
| 478 | ELSE |
|---|
| 479 | C --- Example: the call below is for sensitivities w.r.t. rate coefficients; |
|---|
| 480 | C --- JCOEFF(1:NSENSIT) are the indices of the NSENSIT rate coefficients |
|---|
| 481 | C --- w.r.t. which one differentiates |
|---|
| 482 | CALL dJac_dRcoeff( Y, FIX, U, NCOEFF, JCOEFF, P ) |
|---|
| 483 | END IF |
|---|
| 484 | TIME = Told |
|---|
| 485 | RETURN |
|---|
| 486 | END |
|---|
| 487 | |
|---|
| 488 | |
|---|
| 489 | SUBROUTINE HESS_CHEM(N, T, Y, HESS) |
|---|
| 490 | INCLUDE 'KPP_ROOT_params.h' |
|---|
| 491 | INCLUDE 'KPP_ROOT_global.h' |
|---|
| 492 | INTEGER N |
|---|
| 493 | KPP_REAL Told, T |
|---|
| 494 | KPP_REAL Y(NVAR), HESS(NHESS) |
|---|
| 495 | Told = TIME |
|---|
| 496 | TIME = T |
|---|
| 497 | CALL Update_SUN() |
|---|
| 498 | CALL Update_RCONST() |
|---|
| 499 | CALL Hessian( Y, FIX, RCONST, HESS ) |
|---|
| 500 | TIME = Told |
|---|
| 501 | RETURN |
|---|
| 502 | END |
|---|
| 503 | |
|---|
