?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
added README_changes.txt
[wrffire.git] / wrfv2_fire / chem / KPP / kpp / kpp-2.1 / int / oldies / rodas3_ddm.f
blobf4906c331a87e918c0bb1072359df9f68ae2e6ef
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 RODAS3_DDM(NVAR,NSENSIT,TIN,TOUT,STEPMIN,STEPMAX,
24 + STEPMIN,Y,ATOL,RTOL,
25 + Info,FUNC_CHEM,JAC_CHEM)
26
27
28 RETURN
29 END
30
31
32
33
34 SUBROUTINE RODAS3_DDM(N,NSENSIT,T,Tnext,Hmin,Hmax,Hstart,
35 + y,AbsTol,RelTol,
36 + Info,FUNC_CHEM,JAC_CHEM)
37
38 IMPLICIT NONE
39 INCLUDE 'KPP_ROOT_params.h'
40 INCLUDE 'KPP_ROOT_global.h'
41 INCLUDE 'KPP_ROOT_sparse.h'
42 C
43 C Stiffly accurate Rosenbrock 3(2), with
44 C stiffly accurate embedded formula for error control.
45 C
46 C Direct decoupled computation of sensitivities.
47 C The global variable DDMTYPE distinguishes between:
48 C DDMTYPE = 0 : sensitivities w.r.t. initial values
49 C DDMTYPE = 1 : sensitivities w.r.t. parameters
50 C
51 C INPUT ARGUMENTS:
52 C y = Vector of: (1:NVAR) concentrations, followed by
53 C (1:NVAR) sensitivities w.r.t. first parameter, followed by
54 C etc., followed by
55 C (1:NVAR) sensitivities w.r.t. NSENSIT's parameter
56 C (y contains initial values at input, final values at output)
57 C [T, Tnext] = the integration interval
58 C Hmin, Hmax = lower and upper bounds for the selected step-size.
59 C Note that for Step = Hmin the current computed
60 C solution is unconditionally accepted by the error
61 C control mechanism.
62 C AbsTol, RelTol = (NVAR) dimensional vectors of
63 C componentwise absolute and relative tolerances.
64 C FUNC_CHEM = name of routine of derivatives. KPP syntax.
65 C See the header below.
66 C JAC_CHEM = name of routine that computes the Jacobian, in
67 C sparse format. KPP syntax. See the header below.
68 C Info(1) = 1 for Autonomous system
69 C = 0 for nonAutonomous system
70 C
71 C OUTPUT ARGUMENTS:
72 C y = the values of concentrations and sensitivities at Tend.
73 C T = equals TENDon output.
74 C Info(2) = # of FUNC_CHEM CALLs.
75 C Info(3) = # of JAC_CHEM CALLs.
76 C Info(4) = # of accepted steps.
77 C Info(5) = # of rejected steps.
78 C
79 C Adrian Sandu, December 2001
80 C
81
82 INTEGER NSENSIT
83 KPP_REAL y(NVAR*(NSENSIT+1)), ynew(NVAR*(NSENSIT+1))
84 KPP_REAL K1(NVAR*(NSENSIT+1))
85 KPP_REAL K2(NVAR*(NSENSIT+1))
86 KPP_REAL K3(NVAR*(NSENSIT+1))
87 KPP_REAL K4(NVAR*(NSENSIT+1))
88 KPP_REAL Fv(NVAR), Hv(NVAR)
89 KPP_REAL DFDT(NVAR*(NSENSIT+1))
90 KPP_REAL DJDP(NVAR*NSENSIT)
91 KPP_REAL DFDP(NVAR*NSENSIT), DFDPDT(NVAR*NSENSIT)
92 KPP_REAL JAC(LU_NONZERO), AJAC(LU_NONZERO)
93 KPP_REAL DJDT(LU_NONZERO)
94 KPP_REAL HESS(NHESS)
95 KPP_REAL Hmin,Hmax,Hstart,ghinv,uround
96 KPP_REAL AbsTol(NVAR), RelTol(NVAR)
97 KPP_REAL T, Tnext, Tplus, H, Hnew, elo
98 KPP_REAL ERR, factor, facmax
99 KPP_REAL w, e, beta1, beta2, beta3, beta4
100 KPP_REAL tau, x1, x2, ytol, dround
101
102 INTEGER n,nfcn,njac,Naccept,Nreject,i,j,ier
103 INTEGER Info(5)
104 LOGICAL IsReject, Autonomous
105 EXTERNAL FUNC_CHEM, JAC_CHEM, HESS_CHEM
106
107 C The method coefficients
108 DOUBLE PRECISION gamma, gamma2, gamma3, gamma4
109 PARAMETER ( gamma = 0.5D+00 )
110 PARAMETER ( gamma2 = 1.5D+00 )
111 PARAMETER ( gamma3 = 0.0D+00 )
112 PARAMETER ( gamma4 = 0.0D+00 )
113 DOUBLE PRECISION a21, a31, a32, a41, a42, a43
114 PARAMETER ( a21 = 0.0D+00 )
115 PARAMETER ( a31 = 2.0D+00 )
116 PARAMETER ( a32 = 0.0D+00 )
117 PARAMETER ( a41 = 2.0D+00 )
118 PARAMETER ( a42 = 0.0D+00 )
119 PARAMETER ( a43 = 1.0D+00 )
120 DOUBLE PRECISION alpha2, alpha3, alpha4
121 PARAMETER ( alpha2 = 0.0D0 )
122 PARAMETER ( alpha3 = 1.0D0 )
123 PARAMETER ( alpha4 = 1.0D0 )
124 DOUBLE PRECISION c21, c31, c32, c41, c42, c43
125 PARAMETER ( c21 = 4.0D0 )
126 PARAMETER ( c31 = 1.0D0 )
127 PARAMETER ( c32 = -1.0D0 )
128 PARAMETER ( c41 = 1.0D0 )
129 PARAMETER ( c42 = -1.0D0 )
130 PARAMETER ( c43 = -2.666666666666667D0 )
131 DOUBLE PRECISION b1, b2, b3, b4
132 PARAMETER ( b1 = 2.0D+00 )
133 PARAMETER ( b2 = 0.0D0 )
134 PARAMETER ( b3 = 1.0D0 )
135 PARAMETER ( b4 = 1.0D0 )
136 DOUBLE PRECISION d1, d2, d3, d4
137 PARAMETER ( d1 = 0.0D0 )
138 PARAMETER ( d2 = 0.0D0 )
139 PARAMETER ( d3 = 0.0D0 )
140 PARAMETER ( d4 = 1.0D0 )
141
142
143 c Initialization of counters, etc.
144 Autonomous = Info(1) .EQ. 1
145 uround = 1.d-15
146 dround = DSQRT(uround)
147 IF (Hmax.le.0.D0) THEN
148 Hmax = DABS(Tnext-T)
149 END IF
150 H = DMAX1(1.d-8, Hstart)
151 Tplus = T
152 IsReject = .false.
153 Naccept = 0
154 Nreject = 0
155 Nfcn = 0
156 Njac = 0
157
158 C === Starting the time loop ===
159 10 CONTINUE
160
161 Tplus = T + H
162 IF ( Tplus .gt. Tnext ) THEN
163 H = Tnext - T
164 Tplus = Tnext
165 END IF
166
167 C Initial Function, Jacobian, and Hessian Values
168 CALL FUNC_CHEM(NVAR, T, y, Fv)
169 CALL JAC_CHEM(NVAR, T, y, JAC)
170 CALL HESS_CHEM( NVAR, T, y, HESS )
171 IF (DDMTYPE .EQ. 1) THEN
172 CALL DFUNDPAR(NVAR, NSENSIT, T, y, DFDP)
173 END IF
174
175 C The time derivatives for non-Autonomous case
176 IF (.not. Autonomous) THEN
177 tau = DSIGN(dround*DMAX1( 1.0d0, DABS(T) ), T)
178 CALL FUNC_CHEM(NVAR, T+tau, y, K2)
179 CALL JAC_CHEM(NVAR, T+tau, y, AJAC)
180 nfcn=nfcn+1
181 DO 20 j = 1,NVAR
182 DFDT(j) = ( K2(j)-Fv(j) )/tau
183 20 CONTINUE
184 DO 30 j = 1,LU_NONZERO
185 DJDT(j) = ( AJAC(j)-JAC(j) )/tau
186 30 CONTINUE
187 DO 35 i=1,NSENSIT
188 CALL Jac_SP_Vec (DJDT,y(i*NVAR+1),DFDT(i*NVAR+1))
189 35 CONTINUE
190 END IF
191
192 11 CONTINUE ! From here we restart after a rejected step
193
194 C Form the Prediction matrix and compute its LU factorization
195 Njac = Njac+1
196 ghinv = 1.0d0/(gamma*H)
197 DO 40 j=1,LU_NONZERO
198 AJAC(j) = -JAC(j)
199 40 CONTINUE
200 DO 50 j=1,NVAR
201 AJAC(LU_DIAG(j)) = AJAC(LU_DIAG(j)) + ghinv
202 50 CONTINUE
203 CALL KppDecomp (AJAC, ier)
204 C
205 IF (ier.ne.0) THEN
206 IF ( H.gt.Hmin) THEN
207 H = 5.0d-1*H
208 GO TO 10
209 ELSE
210 PRINT *,'ROS4: Singular factorization at T=',T,'; H=',H
211 STOP
212 END IF
213 END IF
214
215 C ------------ STAGE 1-------------------------
216 DO 60 j = 1,NVAR
217 K1(j) = Fv(j)
218 60 CONTINUE
219 IF (.NOT. Autonomous) THEN
220 beta1 = H*gamma
221 DO 70 j=1,NVAR
222 K1(j) = K1(j) + beta1*DFDT(j)
223 70 CONTINUE
224 END IF
225 CALL KppSolve (AJAC, K1)
226 C --- If derivative w.r.t. parameters
227 IF (DDMTYPE .EQ. 1) THEN
228 CALL DJACDPAR(NVAR, NSENSIT, T, y, K1(1), DJDP)
229 END IF
230 C --- End of derivative w.r.t. parameters
231
232 DO 100 i=1,NSENSIT
233 CALL Jac_SP_Vec (JAC,y(i*NVAR+1),K1(i*NVAR+1))
234 CALL Hess_Vec ( HESS, y(i*NVAR+1), K1(1), Hv )
235 DO 80 j=1,NVAR
236 K1(i*NVAR+j) = K1(i*NVAR+j) + Hv(j)
237 80 CONTINUE
238 IF (.NOT. Autonomous) THEN
239 DO 90 j=1,NVAR
240 K1(i*NVAR+j) = K1(i*NVAR+j) + beta1*DFDT(i*NVAR+j)
241 90 CONTINUE
242 END IF
243 C --- If derivative w.r.t. parameters
244 IF (DDMTYPE .EQ. 1) THEN
245 DO 95 j = 1,NVAR
246 K1(i*NVAR+j) = K1(i*NVAR+j) + DFDP((i-1)*NVAR+j)
247 & + DJDP((i-1)*NVAR+j)
248 95 CONTINUE
249 END IF
250 C --- End of derivative w.r.t. parameters
251 CALL KppSolve (AJAC, K1(i*NVAR+1))
252 100 CONTINUE
253
254 C ----------- STAGE 2 -------------------------
255 C Note: uses the same function values as Stage 1
256 C DO 110 j = 1,NVAR*(NSENSIT+1)
257 C ynew(j) = y(j) + a21*K1(j)
258 C 110 CONTINUE
259 C CALL FUNC_CHEM(NVAR, T+alpha2*H, ynew, Fv)
260 C IF (DDMTYPE .EQ. 1) THEN
261 C CALL DFUNDPAR(NVAR, NSENSIT, T+alpha2*H, ynew, DFDP)
262 C END IF
263 C nfcn=nfcn+1
264 beta1 = c21/H
265 DO 120 j = 1,NVAR
266 K2(j) = Fv(j) + beta1*K1(j)
267 120 CONTINUE
268 IF (.NOT. Autonomous) THEN
269 beta2 = H*gamma2
270 DO 130 j=1,NVAR
271 K2(j) = K2(j) + beta2*DFDT(j)
272 130 CONTINUE
273 END IF
274 CALL KppSolve (AJAC, K2)
275 C --- If derivative w.r.t. parameters
276 IF (DDMTYPE .EQ. 1) THEN
277 CALL DJACDPAR(NVAR, NSENSIT, T, y, K2(1), DJDP)
278 END IF
279 C --- End of derivative w.r.t. parameters
280
281 CALL JAC_CHEM(NVAR, T+alpha2*H, ynew, JAC)
282 njac=njac+1
283 DO 160 i=1,NSENSIT
284 CALL Jac_SP_Vec (JAC,ynew(i*NVAR+1),K2(i*NVAR+1))
285 CALL Hess_Vec ( HESS, y(i*NVAR+1), K2(1), Hv )
286 DO 140 j = 1,NVAR
287 K2(i*NVAR+j) = K2(i*NVAR+j) + beta1*K1(i*NVAR+j)
288 & + Hv(j)
289 140 CONTINUE
290 IF (.NOT. Autonomous) THEN
291 DO 150 j=1,NVAR
292 K2(i*NVAR+j) = K2(i*NVAR+j) + beta2*DFDT(i*NVAR+j)
293 150 CONTINUE
294 END IF
295 C --- If derivative w.r.t. parameters
296 IF (DDMTYPE .EQ. 1) THEN
297 DO 155 j = 1,NVAR
298 K2(i*NVAR+j) = K2(i*NVAR+j) + DFDP((i-1)*NVAR+j)
299 & + DJDP((i-1)*NVAR+j)
300 155 CONTINUE
301 END IF
302 C --- End of derivative w.r.t. parameters
303 CALL KppSolve (AJAC, K2(i*NVAR+1))
304 160 CONTINUE
305
306
307 C ------------ STAGE 3 -------------------------
308 DO 170 j = 1,NVAR*(NSENSIT+1)
309 ynew(j) = y(j) + a31*K1(j) + a32*K2(j)
310 170 CONTINUE
311 CALL FUNC_CHEM(NVAR, T+alpha3*H, ynew, Fv)
312 IF (DDMTYPE .EQ. 1) THEN
313 CALL DFUNDPAR(NVAR, NSENSIT, T+alpha3*H, ynew, DFDP)
314 END IF
315 nfcn=nfcn+1
316 beta1 = c31/H
317 beta2 = c32/H
318 DO 180 j = 1,NVAR
319 K3(j) = Fv(j) + beta1*K1(j) + beta2*K2(j)
320 180 CONTINUE
321 IF (.NOT. Autonomous) THEN
322 beta3 = H*gamma3
323 DO 190 j=1,NVAR
324 K3(j) = K3(j) + beta3*DFDT(j)
325 190 CONTINUE
326 END IF
327 CALL KppSolve (AJAC, K3)
328 C --- If derivative w.r.t. parameters
329 IF (DDMTYPE .EQ. 1) THEN
330 CALL DJACDPAR(NVAR, NSENSIT, T, y, K3(1), DJDP)
331 END IF
332 C --- End of derivative w.r.t. parameters
333
334 CALL JAC_CHEM(NVAR, T+alpha3*H, ynew, JAC)
335 njac=njac+1
336 DO 220 i=1,NSENSIT
337 CALL Jac_SP_Vec (JAC,ynew(i*NVAR+1),K3(i*NVAR+1))
338 CALL Hess_Vec ( HESS, y(i*NVAR+1), K3(1), Hv )
339 DO 200 j = 1,NVAR
340 K3(i*NVAR+j) = K3(i*NVAR+j) + beta1*K1(i*NVAR+j)
341 & + beta2*K2(i*NVAR+j) + Hv(j)
342 200 CONTINUE
343 IF (.NOT. Autonomous) THEN
344 DO 210 j=1,NVAR
345 K3(i*NVAR+j) = K3(i*NVAR+j) + beta3*DFDT(i*NVAR+j)
346 210 CONTINUE
347 END IF
348 C --- If derivative w.r.t. parameters
349 IF (DDMTYPE .EQ. 1) THEN
350 DO 215 j = 1,NVAR
351 K3(i*NVAR+j) = K3(i*NVAR+j) + DFDP((i-1)*NVAR+j)
352 & + DJDP((i-1)*NVAR+j)
353 215 CONTINUE
354 END IF
355 C --- End of derivative w.r.t. parameters
356 CALL KppSolve (AJAC, K3(i*NVAR+1))
357 220 CONTINUE
358
359 C ------------ STAGE 4 -------------------------
360 DO 225 j = 1,NVAR*(NSENSIT+1)
361 ynew(j) = y(j) + a41*K1(j) + a42*K2(j) + a43*K3(j)
362 225 CONTINUE
363 CALL FUNC_CHEM(NVAR, T+alpha4*H, ynew, Fv)
364 IF (DDMTYPE .EQ. 1) THEN
365 CALL DFUNDPAR(NVAR, NSENSIT, T+alpha4*H, ynew, DFDP)
366 END IF
367 nfcn=nfcn+1
368 beta1 = c41/H
369 beta2 = c42/H
370 beta3 = c43/H
371 DO 230 j = 1,NVAR
372 K4(j) = Fv(j) + beta1*K1(j) + beta2*K2(j) + beta3*K3(j)
373 230 CONTINUE
374 IF (.NOT. Autonomous) THEN
375 beta4 = H*gamma4
376 DO 240 j=1,NVAR
377 K4(j) = K4(j) + beta4*DFDT(j)
378 240 CONTINUE
379 END IF
380 CALL KppSolve (AJAC, K4)
381 C --- If derivative w.r.t. parameters
382 IF (DDMTYPE .EQ. 1) THEN
383 CALL DJACDPAR(NVAR, NSENSIT, T, y, K4(1), DJDP)
384 END IF
385 C --- End of derivative w.r.t. parameters
386
387 njac=njac+1
388 DO 270 i=1,NSENSIT
389 CALL Jac_SP_Vec (JAC,ynew(i*NVAR+1),K4(i*NVAR+1))
390 CALL Hess_Vec ( HESS, y(i*NVAR+1), K4(1), Hv )
391 DO 250 j = 1,NVAR
392 K4(i*NVAR+j) = K4(i*NVAR+j) + beta1*K1(i*NVAR+j)
393 & + beta2*K2(i*NVAR+j) + beta3*K3(i*NVAR+j)
394 & + Hv(j)
395 250 CONTINUE
396 IF (.NOT. Autonomous) THEN
397 DO 260 j=1,NVAR
398 K4(i*NVAR+j) = K4(i*NVAR+j) + beta4*DFDT(i*NVAR+j)
399 260 CONTINUE
400 END IF
401 C --- If derivative w.r.t. parameters
402 IF (DDMTYPE .EQ. 1) THEN
403 DO 265 j = 1,NVAR
404 K4(i*NVAR+j) = K4(i*NVAR+j) + DFDP((i-1)*NVAR+j)
405 & + DJDP((i-1)*NVAR+j)
406 265 CONTINUE
407 END IF
408 CALL KppSolve (AJAC, K4(i*NVAR+1))
409 270 CONTINUE
410
411
412 C ---- The Solution ---
413 DO 280 j = 1,NVAR*(NSENSIT+1)
414 C ynew(j) = y(j) + b1*K1(j) + b2*K2(j) + b3*K3(j) + b4*K4(j)
415 ynew(j) = y(j) + 2*K1(j) + K3(j) + K4(j)
416 280 CONTINUE
417
418
419 C ====== Error estimation -- can be extended to control sensitivities too ========
420
421 ERR = 0.d0
422 DO 290 i=1,NVAR
423 w = AbsTol(i) + RelTol(i)*DMAX1(DABS(ynew(i)),DABS(y(i)))
424 C e = d1*K1(i) + d2*K2(i) + d3*K3(i) + d4*K4(i)
425 e = K4(i)
426 ERR = ERR + ( e/w )**2
427 290 CONTINUE
428 ERR = DMAX1( uround, DSQRT( ERR/NVAR ) )
429
430 C ======= Choose the stepsize ===============================
431
432 elo = 3.0D0 ! estimator local order
433 factor = DMAX1(2.0D-1,DMIN1(6.0D0,ERR**(1.0D0/elo)/.9D0))
434 Hnew = DMIN1(Hmax,DMAX1(Hmin, H/factor))
435
436 C ======= Rejected/Accepted Step ============================
437
438 IF ( (ERR.gt.1).and.(H.gt.Hmin) ) THEN
439 IsReject = .true.
440 H = DMIN1(H/10,Hnew)
441 Nreject = Nreject+1
442 ELSE
443 DO 300 i=1,NVAR*(NSENSIT+1)
444 y(i) = ynew(i)
445 300 CONTINUE
446 T = Tplus
447 IF (.NOT.IsReject) THEN
448 H = Hnew ! Do not increase stepsize if previos step was rejected
449 END IF
450 IsReject = .false.
451 Naccept = Naccept+1
452 END IF
453
454 C ======= End of the time loop ===============================
455 IF ( T .lt. Tnext ) GO TO 10
456
457
458
459 C ======= Output Information =================================
460 Info(2) = Nfcn
461 Info(3) = Njac
462 Info(4) = Naccept
463 Info(5) = Nreject
464 Hstart = H
465
466 RETURN
467 END
468
469
470
471 SUBROUTINE FUNC_CHEM(N, T, Y, P)
472 INCLUDE 'KPP_ROOT_params.h'
473 INCLUDE 'KPP_ROOT_global.h'
474 INTEGER N
475 KPP_REAL T, Told
476 KPP_REAL Y(NVAR), P(NVAR)
477 Told = TIME
478 TIME = T
479 CALL Update_SUN()
480 CALL Update_RCONST()
481 CALL Fun( Y, FIX, RCONST, P )
482 TIME = Told
483 RETURN
484 END
485
486
487 SUBROUTINE DFUNDPAR(N, NSENSIT, T, Y, P)
488 C --- Computes the partial derivatives of FUNC_CHEM w.r.t. parameters
489 INCLUDE 'KPP_ROOT_params.h'
490 INCLUDE 'KPP_ROOT_global.h'
491 C --- NCOEFF, JCOEFF useful for derivatives w.r.t. rate coefficients
492 INTEGER N
493 INTEGER NCOEFF, JCOEFF(NREACT)
494 COMMON /DDMRCOEFF/ NCOEFF, JCOEFF
495
496 KPP_REAL T, Told
497 KPP_REAL Y(NVAR), P(NVAR*NSENSIT)
498 Told = TIME
499 TIME = T
500 CALL Update_SUN()
501 CALL Update_RCONST()
502 C
503 IF (DDMTYPE .EQ. 0) THEN
504 C --- Note: the values below are for sensitivities w.r.t. initial values;
505 C --- they may have to be changed for other applications
506 DO j=1,NSENSIT
507 DO i=1,NVAR
508 P(i+NVAR*(j-1)) = 0.0D0
509 END DO
510 END DO
511 ELSE
512 C --- Example: the call below is for sensitivities w.r.t. rate coefficients;
513 C --- JCOEFF(1:NSENSIT) are the indices of the NSENSIT rate coefficients
514 C --- w.r.t. which one differentiates
515 CALL dFun_dRcoeff( Y, FIX, NCOEFF, JCOEFF, P )
516 END IF
517 TIME = Told
518 RETURN
519 END
520
521 SUBROUTINE JAC_CHEM(N, T, Y, J)
522 INCLUDE 'KPP_ROOT_params.h'
523 INCLUDE 'KPP_ROOT_global.h'
524 INTEGER N
525 KPP_REAL Told, T
526 KPP_REAL Y(NVAR), J(LU_NONZERO)
527 Told = TIME
528 TIME = T
529 CALL Update_SUN()
530 CALL Update_RCONST()
531 CALL Jac_SP( Y, FIX, RCONST, J )
532 TIME = Told
533 RETURN
534 END
535
536
537 SUBROUTINE DJACDPAR(N, NSENSIT, T, Y, U, P)
538 C --- Computes the partial derivatives of JAC w.r.t. parameters times user vector U
539 INCLUDE 'KPP_ROOT_params.h'
540 INCLUDE 'KPP_ROOT_global.h'
541 C --- NCOEFF, JCOEFF useful for derivatives w.r.t. rate coefficients
542 INTEGER NCOEFF, JCOEFF(NREACT)
543 COMMON /DDMRCOEFF/ NCOEFF, JCOEFF
544
545 INTEGER N
546 KPP_REAL T, Told
547 KPP_REAL Y(NVAR), U(NVAR)
548 KPP_REAL P(NVAR*NSENSIT)
549 Told = TIME
550 TIME = T
551 CALL Update_SUN()
552 CALL Update_RCONST()
553 C
554 IF (DDMTYPE .EQ. 0) THEN
555 C --- Note: the values below are for sensitivities w.r.t. initial values;
556 C --- they may have to be changed for other applications
557 DO j=1,NSENSIT
558 DO i=1,NVAR
559 P(i+NVAR*(j-1)) = 0.0D0
560 END DO
561 END DO
562 ELSE
563 C --- Example: the call below is for sensitivities w.r.t. rate coefficients;
564 C --- JCOEFF(1:NSENSIT) are the indices of the NSENSIT rate coefficients
565 C --- w.r.t. which one differentiates
566 CALL dJac_dRcoeff( Y, FIX, U, NCOEFF, JCOEFF, P )
567 END IF
568 TIME = Told
569 RETURN
570 END
571
572
573 SUBROUTINE HESS_CHEM(N, T, Y, HESS)
574 INCLUDE 'KPP_ROOT_params.h'
575 INCLUDE 'KPP_ROOT_global.h'
576 INTEGER N
577 KPP_REAL Told, T
578 KPP_REAL Y(NVAR), HESS(NHESS)
579 Told = TIME
580 TIME = T
581 CALL Update_SUN()
582 CALL Update_RCONST()
583 CALL Hessian( Y, FIX, RCONST, HESS )
584 TIME = Told
585 RETURN
586 END
587
588
589
590
591
592
A coupled weather-fire forecasting model built on top of Weather Research and Forecasting (WRF).