?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Merge branch 'master' into devel
[wrffire.git] / wrfv2_fire / chem / KPP / kpp / kpp-2.1 / int / oldies / ros4_ddm.f
blob79a939c81f744d2fe235803f34797d560974525c
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 ROS4_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 ROS4_DDM(N,NSENSIT,T,Tnext,Hmin,Hmax,Hstart,
35 + y,AbsTol,RelTol,
36 + Info,FUNC_CHEM,JAC_CHEM)
37 IMPLICIT NONE
38 INCLUDE 'KPP_ROOT_params.h'
39 INCLUDE 'KPP_ROOT_global.h'
40 INCLUDE 'KPP_ROOT_sparse.h'
41 C
42 C Four Stages, Fourth Order L-stable Rosenbrock Method,
43 C with embedded L-stable, third order method for error control
44 C Simplified version of E. Hairer's atmros4; the coefficients are slightly different
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
83 INTEGER NSENSIT
84 KPP_REAL y(NVAR*(NSENSIT+1)), ynew(NVAR*(NSENSIT+1))
85 KPP_REAL K1(NVAR*(NSENSIT+1))
86 KPP_REAL K2(NVAR*(NSENSIT+1))
87 KPP_REAL K3(NVAR*(NSENSIT+1))
88 KPP_REAL K4(NVAR*(NSENSIT+1))
89 KPP_REAL Fv(NVAR), Hv(NVAR)
90 KPP_REAL DFDT(NVAR*(NSENSIT+1))
91 KPP_REAL DFDP(NVAR*NSENSIT), DFDPDT(NVAR*NSENSIT)
92 KPP_REAL DJDP(NVAR*NSENSIT)
93 KPP_REAL JAC(LU_NONZERO), AJAC(LU_NONZERO)
94 KPP_REAL DJDT(LU_NONZERO)
95 KPP_REAL HESS(NHESS)
96 KPP_REAL Hmin,Hmax,Hstart,ghinv,uround
97 KPP_REAL AbsTol(NVAR), RelTol(NVAR)
98 KPP_REAL T, Tnext, Tplus, H, Hnew, elo
99 KPP_REAL ERR, factor, facmax, dround, tau
100 KPP_REAL w, e, beta1, beta2, beta3, beta4
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
108 C The method coefficients
109 C DOUBLE PRECISION gamma, gamma2, gamma3, gamma4
110 C PARAMETER ( gamma = 0.57281606D0 )
111 C PARAMETER ( gamma2 = -1.769177067112013949170520D0 )
112 C PARAMETER ( gamma3 = 0.759293964293209853670967D0 )
113 C PARAMETER ( gamma4 = -0.104894621490955803206743D0 )
114 C DOUBLE PRECISION a21, a31, a32, a41, a42, a43
115 C PARAMETER ( a21 = 2.00000000000000000000000D0 )
116 C PARAMETER ( a31 = 1.86794814949823713234476D0 )
117 C PARAMETER ( a32 = 0.23444556851723885002322D0 )
118 C DOUBLE PRECISION alpha2, alpha3, alpha4
119 C PARAMETER ( alpha2 = 1.145632120D0 )
120 C PARAMETER ( alpha3 = 0.655214975973133829477748D0 )
121 C DOUBLE PRECISION c21, c31, c32, c41, c42, c43
122 C PARAMETER ( c21 = -7.137649943349979830369260D0 )
123 C PARAMETER ( c31 = 2.580923666509657714488050D0 )
124 C PARAMETER ( c32 = 0.651629887302032023387417D0 )
125 C PARAMETER ( c41 = -2.137115266506619116806370D0 )
126 C PARAMETER ( c42 = -0.321469531339951070769241D0 )
127 C PARAMETER ( c43 = -0.694966049282445225157329D0 )
128 C DOUBLE PRECISION m1, m2, m3, m4, mhat1, mhat2, mhat3, mhat4
129 C PARAMETER ( m1 = 2.255566228604565243728840D0 )
130 C PARAMETER ( m2 = 0.287055063194157607662630D0 )
131 C PARAMETER ( m3 = 0.435311963379983213402707D0 )
132 C PARAMETER ( m4 = 1.093507656403247803214820D0 )
133 C PARAMETER ( mhat1 = 2.068399160527583734258670D0 )
134 C PARAMETER ( mhat2 = 0.238681352067532797956493D0 )
135 C PARAMETER ( mhat3 = 0.363373345435391708261747D0 )
136 C PARAMETER ( mhat4 = 0.366557127936155144309163D0 )
137 C DOUBLE PRECISION e1, e2, e3, e4
138 c PARAMETER ( e1 = 1.8716706807698191283861888D-01 )
139 c PARAMETER ( e2 = 4.8373711126624835410225955D-02 )
140 c PARAMETER ( e3 = 7.1938617944591554120847832D-02 )
141 c PARAMETER ( e4 = 7.2695052846709262706070831D-01 )
142 C PARAMETER ( e1 = -0.2815431932141155D+00 )
143 C PARAMETER ( e2 = -0.7276199124938920D-01 )
144 C PARAMETER ( e3 = -0.1082196201495311D+00 )
145 C PARAMETER ( e4 = -0.1093502252409163D+01 )
146 C The method coefficients
147 DOUBLE PRECISION gamma, gamma2, gamma3, gamma4
148 PARAMETER ( gamma = 0.5728200000000000D+00 )
149 PARAMETER ( gamma2 = -0.1769193891319233D+01 )
150 PARAMETER ( gamma3 = 0.7592633437920482D+00 )
151 PARAMETER ( gamma4 = -0.1049021087100450D+00 )
152 DOUBLE PRECISION a21, a31, a32, a41, a42, a43
153 PARAMETER ( a21 = 0.2000000000000000D+01 )
154 PARAMETER ( a31 = 0.1867943637803922D+01 )
155 PARAMETER ( a32 = 0.2344449711399156D+00 )
156 DOUBLE PRECISION alpha2, alpha3
157 PARAMETER ( alpha2 = 0.1145640000000000D+01 )
158 PARAMETER ( alpha3 = 0.6552168638155900D+00 )
159 DOUBLE PRECISION c21, c31, c32, c41, c42, c43
160 PARAMETER ( c21 = -0.7137615036412310D+01 )
161 PARAMETER ( c31 = 0.2580708087951457D+01 )
162 PARAMETER ( c32 = 0.6515950076447975D+00 )
163 PARAMETER ( c41 = -0.2137148994382534D+01 )
164 PARAMETER ( c42 = -0.3214669691237626D+00 )
165 PARAMETER ( c43 = -0.6949742501781779D+00 )
166 DOUBLE PRECISION b1, b2, b3, b4
167 PARAMETER ( b1 = 0.2255570073418735D+01 )
168 PARAMETER ( b2 = 0.2870493262186792D+00 )
169 PARAMETER ( b3 = 0.4353179431840180D+00 )
170 PARAMETER ( b4 = 0.1093502252409163D+01 )
171 DOUBLE PRECISION d1, d2, d3, d4
172 PARAMETER ( d1 = -0.2815431932141155D+00 )
173 PARAMETER ( d2 = -0.7276199124938920D-01 )
174 PARAMETER ( d3 = -0.1082196201495311D+00 )
175 PARAMETER ( d4 = -0.1093502252409163D+01 )
176
177
178 c Initialization of counters, etc.
179 Autonomous = Info(1) .EQ. 1
180 uround = 1.d-15
181 dround = DSQRT(uround)
182 IF (Hmax.le.0.D0) THEN
183 Hmax = DABS(Tnext-T)
184 END IF
185 H = DMAX1(1.d-8, Hstart)
186 Tplus = T
187 IsReject = .false.
188 Naccept = 0
189 Nreject = 0
190 Nfcn = 0
191 Njac = 0
192
193 C === Starting the time loop ===
194 10 CONTINUE
195
196 Tplus = T + H
197 IF ( Tplus .gt. Tnext ) THEN
198 H = Tnext - T
199 Tplus = Tnext
200 END IF
201
202 C Initial Function, Jacobian, and Hessian Values
203 CALL FUNC_CHEM(NVAR, T, y, Fv)
204 CALL JAC_CHEM(NVAR, T, y, JAC)
205 CALL HESS_CHEM( NVAR, T, y, HESS )
206 IF (DDMTYPE .EQ. 1) THEN
207 CALL DFUNDPAR(NVAR, NSENSIT, T, y, DFDP)
208 END IF
209
210 C The time derivatives for non-Autonomous case
211 IF (.not. Autonomous) THEN
212 tau = DSIGN(dround*DMAX1( 1.0d0, DABS(T) ), T)
213 CALL FUNC_CHEM(NVAR, T+tau, y, K2)
214 CALL JAC_CHEM(NVAR, T+tau, y, AJAC)
215 nfcn=nfcn+1
216 DO 20 j = 1,NVAR
217 DFDT(j) = ( K2(j)-Fv(j) )/tau
218 20 CONTINUE
219 DO 30 j = 1,LU_NONZERO
220 DJDT(j) = ( AJAC(j)-JAC(j) )/tau
221 30 CONTINUE
222 DO 35 i=1,NSENSIT
223 CALL Jac_SP_Vec (DJDT,y(i*NVAR+1),DFDT(i*NVAR+1))
224 35 CONTINUE
225 END IF
226
227 11 CONTINUE ! From here we restart after a rejected step
228
229 C Form the Prediction matrix and compute its LU factorization
230 Njac = Njac+1
231 ghinv = 1.0d0/(gamma*H)
232 DO 40 j=1,LU_NONZERO
233 AJAC(j) = -JAC(j)
234 40 CONTINUE
235 DO 50 j=1,NVAR
236 AJAC(LU_DIAG(j)) = AJAC(LU_DIAG(j)) + ghinv
237 50 CONTINUE
238 CALL KppDecomp (AJAC, ier)
239 C
240 IF (ier.ne.0) THEN
241 IF ( H.gt.Hmin) THEN
242 H = 5.0d-1*H
243 GO TO 10
244 ELSE
245 PRINT *,'ROS4: Singular factorization at T=',T,'; H=',H
246 STOP
247 END IF
248 END IF
249
250 C ------------ STAGE 1-------------------------
251 DO 60 j = 1,NVAR
252 K1(j) = Fv(j)
253 60 CONTINUE
254 IF (.NOT. Autonomous) THEN
255 beta1 = H*gamma
256 DO 70 j=1,NVAR
257 K1(j) = K1(j) + beta1*DFDT(j)
258 70 CONTINUE
259 END IF
260 CALL KppSolve (AJAC, K1)
261 C --- If derivative w.r.t. parameters
262 IF (DDMTYPE .EQ. 1) THEN
263 CALL DJACDPAR(NVAR, NSENSIT, T, y, K1(1), DJDP)
264 END IF
265 C --- End of derivative w.r.t. parameters
266
267 DO 100 i=1,NSENSIT
268 CALL Jac_SP_Vec (JAC,y(i*NVAR+1),K1(i*NVAR+1))
269 CALL Hess_Vec ( HESS, y(i*NVAR+1), K1(1), Hv )
270 DO 80 j=1,NVAR
271 K1(i*NVAR+j) = K1(i*NVAR+j) + Hv(j)
272 80 CONTINUE
273 IF (.NOT. Autonomous) THEN
274 DO 90 j=1,NVAR
275 K1(i*NVAR+j) = K1(i*NVAR+j) + beta1*DFDT(i*NVAR+j)
276 90 CONTINUE
277 END IF
278 C --- If derivative w.r.t. parameters
279 IF (DDMTYPE .EQ. 1) THEN
280 DO 95 j = 1,NVAR
281 K1(i*NVAR+j) = K1(i*NVAR+j) + DFDP((i-1)*NVAR+j)
282 & + DJDP((i-1)*NVAR+j)
283 95 CONTINUE
284 END IF
285 C --- End of derivative w.r.t. parameters
286 CALL KppSolve (AJAC, K1(i*NVAR+1))
287 100 CONTINUE
288
289 C ----------- STAGE 2 -------------------------
290 DO 110 j = 1,NVAR*(NSENSIT+1)
291 ynew(j) = y(j) + a21*K1(j)
292 110 CONTINUE
293 CALL FUNC_CHEM(NVAR, T+alpha2*H, ynew, Fv)
294 IF (DDMTYPE .EQ. 1) THEN
295 CALL DFUNDPAR(NVAR, NSENSIT, T+alpha2*H, ynew, DFDP)
296 END IF
297 nfcn=nfcn+1
298 beta1 = c21/H
299 DO 120 j = 1,NVAR
300 K2(j) = Fv(j) + beta1*K1(j)
301 120 CONTINUE
302 IF (.NOT. Autonomous) THEN
303 beta2 = H*gamma2
304 DO 130 j=1,NVAR
305 K2(j) = K2(j) + beta2*DFDT(j)
306 130 CONTINUE
307 END IF
308 CALL KppSolve (AJAC, K2)
309 C --- If derivative w.r.t. parameters
310 IF (DDMTYPE .EQ. 1) THEN
311 CALL DJACDPAR(NVAR, NSENSIT, T, y, K2(1), DJDP)
312 END IF
313 C --- End of derivative w.r.t. parameters
314
315 CALL JAC_CHEM(NVAR, T+alpha2*H, ynew, JAC)
316 njac=njac+1
317 DO 160 i=1,NSENSIT
318 CALL Jac_SP_Vec (JAC,ynew(i*NVAR+1),K2(i*NVAR+1))
319 CALL Hess_Vec ( HESS, y(i*NVAR+1), K2(1), Hv )
320 DO 140 j = 1,NVAR
321 K2(i*NVAR+j) = K2(i*NVAR+j) + beta1*K1(i*NVAR+j)
322 & + Hv(j)
323 140 CONTINUE
324 IF (.NOT. Autonomous) THEN
325 DO 150 j=1,NVAR
326 K2(i*NVAR+j) = K2(i*NVAR+j) + beta2*DFDT(i*NVAR+j)
327 150 CONTINUE
328 END IF
329 C --- If derivative w.r.t. parameters
330 IF (DDMTYPE .EQ. 1) THEN
331 DO 155 j = 1,NVAR
332 K2(i*NVAR+j) = K2(i*NVAR+j) + DFDP((i-1)*NVAR+j)
333 & + DJDP((i-1)*NVAR+j)
334 155 CONTINUE
335 END IF
336 C --- End of derivative w.r.t. parameters
337 CALL KppSolve (AJAC, K2(i*NVAR+1))
338 160 CONTINUE
339
340
341 C ------------ STAGE 3 -------------------------
342 DO 170 j = 1,NVAR*(NSENSIT+1)
343 ynew(j) = y(j) + a31*K1(j) + a32*K2(j)
344 170 CONTINUE
345 CALL FUNC_CHEM(NVAR, T+alpha3*H, ynew, Fv)
346 IF (DDMTYPE .EQ. 1) THEN
347 CALL DFUNDPAR(NVAR, NSENSIT, T+alpha3*H, ynew, DFDP)
348 END IF
349 nfcn=nfcn+1
350 beta1 = c31/H
351 beta2 = c32/H
352 DO 180 j = 1,NVAR
353 K3(j) = Fv(j) + beta1*K1(j) + beta2*K2(j)
354 180 CONTINUE
355 IF (.NOT. Autonomous) THEN
356 beta3 = H*gamma3
357 DO 190 j=1,NVAR
358 K3(j) = K3(j) + beta3*DFDT(j)
359 190 CONTINUE
360 END IF
361 CALL KppSolve (AJAC, K3)
362 C --- If derivative w.r.t. parameters
363 IF (DDMTYPE .EQ. 1) THEN
364 CALL DJACDPAR(NVAR, NSENSIT, T, y, K3(1), DJDP)
365 END IF
366 C --- End of derivative w.r.t. parameters
367
368 CALL JAC_CHEM(NVAR, T+alpha3*H, ynew, JAC)
369 njac=njac+1
370 DO 220 i=1,NSENSIT
371 CALL Jac_SP_Vec (JAC,ynew(i*NVAR+1),K3(i*NVAR+1))
372 CALL Hess_Vec ( HESS, y(i*NVAR+1), K3(1), Hv )
373 DO 200 j = 1,NVAR
374 K3(i*NVAR+j) = K3(i*NVAR+j) + beta1*K1(i*NVAR+j)
375 & + beta2*K2(i*NVAR+j) + Hv(j)
376 200 CONTINUE
377 IF (.NOT. Autonomous) THEN
378 DO 210 j=1,NVAR
379 K3(i*NVAR+j) = K3(i*NVAR+j) + beta3*DFDT(i*NVAR+j)
380 210 CONTINUE
381 END IF
382 C --- If derivative w.r.t. parameters
383 IF (DDMTYPE .EQ. 1) THEN
384 DO 215 j = 1,NVAR
385 K3(i*NVAR+j) = K3(i*NVAR+j) + DFDP((i-1)*NVAR+j)
386 & + DJDP((i-1)*NVAR+j)
387 215 CONTINUE
388 END IF
389 C --- End of derivative w.r.t. parameters
390 CALL KppSolve (AJAC, K3(i*NVAR+1))
391 220 CONTINUE
392
393 C ------------ STAGE 4 -------------------------
394 C Note: uses the same function values as stage 3
395 beta1 = c41/H
396 beta2 = c42/H
397 beta3 = c43/H
398 DO 230 j = 1,NVAR
399 K4(j) = Fv(j) + beta1*K1(j) + beta2*K2(j) + beta3*K3(j)
400 230 CONTINUE
401 IF (.NOT. Autonomous) THEN
402 beta4 = H*gamma4
403 DO 240 j=1,NVAR
404 K4(j) = K4(j) + beta4*DFDT(j)
405 240 CONTINUE
406 END IF
407 CALL KppSolve (AJAC, K4)
408 C --- If derivative w.r.t. parameters
409 IF (DDMTYPE .EQ. 1) THEN
410 CALL DJACDPAR(NVAR, NSENSIT, T, y, K4(1), DJDP)
411 END IF
412 C --- End of derivative w.r.t. parameters
413
414 njac=njac+1
415 DO 270 i=1,NSENSIT
416 CALL Jac_SP_Vec (JAC,ynew(i*NVAR+1),K4(i*NVAR+1))
417 CALL Hess_Vec ( HESS, y(i*NVAR+1), K4(1), Hv )
418 DO 250 j = 1,NVAR
419 K4(i*NVAR+j) = K4(i*NVAR+j) + beta1*K1(i*NVAR+j)
420 & + beta2*K2(i*NVAR+j) + beta3*K3(i*NVAR+j)
421 & + Hv(j)
422 250 CONTINUE
423 IF (.NOT. Autonomous) THEN
424 DO 260 j=1,NVAR
425 K4(i*NVAR+j) = K4(i*NVAR+j) + beta4*DFDT(i*NVAR+j)
426 260 CONTINUE
427 END IF
428 C --- If derivative w.r.t. parameters
429 IF (DDMTYPE .EQ. 1) THEN
430 DO 265 j = 1,NVAR
431 K4(i*NVAR+j) = K4(i*NVAR+j) + DFDP((i-1)*NVAR+j)
432 & + DJDP((i-1)*NVAR+j)
433 265 CONTINUE
434 END IF
435 CALL KppSolve (AJAC, K4(i*NVAR+1))
436 270 CONTINUE
437
438
439 C ---- The Solution ---
440 DO 280 j = 1,NVAR*(NSENSIT+1)
441 ynew(j) = y(j) + b1*K1(j) + b2*K2(j) + b3*K3(j) + b4*K4(j)
442 280 CONTINUE
443
444
445 C ====== Error estimation -- can be extended to control sensitivities too ========
446
447 ERR = 0.d0
448 DO 290 i=1,NVAR
449 w = AbsTol(i) + RelTol(i)*DMAX1(DABS(ynew(i)),DABS(y(i)))
450 e = d1*K1(i) + d2*K2(i) + d3*K3(i) + d4*K4(i)
451 ERR = ERR + ( e/w )**2
452 290 CONTINUE
453 ERR = DMAX1( uround, DSQRT( ERR/NVAR ) )
454
455 C ======= Choose the stepsize ===============================
456
457 elo = 4.0D0 ! estimator local order
458 factor = DMAX1(2.0D-1,DMIN1(6.0D0,ERR**(1.0D0/elo)/.9D0))
459 Hnew = DMIN1(Hmax,DMAX1(Hmin, H/factor))
460
461 C ======= Rejected/Accepted Step ============================
462
463 IF ( (ERR.gt.1).and.(H.gt.Hmin) ) THEN
464 IsReject = .true.
465 H = DMIN1(H/10,Hnew)
466 Nreject = Nreject+1
467 ELSE
468 DO 300 i=1,NVAR*(NSENSIT+1)
469 y(i) = ynew(i)
470 300 CONTINUE
471 T = Tplus
472 IF (.NOT.IsReject) THEN
473 H = Hnew ! Do not increase stepsize if previos step was rejected
474 END IF
475 IsReject = .false.
476 Naccept = Naccept+1
477 END IF
478
479 C ======= End of the time loop ===============================
480 IF ( T .lt. Tnext ) GO TO 10
481
482
483
484 C ======= Output Information =================================
485 Info(2) = Nfcn
486 Info(3) = Njac
487 Info(4) = Naccept
488 Info(5) = Nreject
489 Hstart = H
490
491 RETURN
492 END
493
494
495
496 SUBROUTINE FUNC_CHEM(N, T, Y, P)
497 INCLUDE 'KPP_ROOT_params.h'
498 INCLUDE 'KPP_ROOT_global.h'
499 KPP_REAL T, Told
500 KPP_REAL Y(NVAR), P(NVAR)
501 Told = TIME
502 TIME = T
503 CALL Update_SUN()
504 CALL Update_RCONST()
505 CALL Fun( Y, FIX, RCONST, P )
506 TIME = Told
507 RETURN
508 END
509
510
511 SUBROUTINE DFUNDPAR(N, NSENSIT, T, Y, P)
512 C --- Computes the partial derivatives of FUNC_CHEM w.r.t. parameters
513 INCLUDE 'KPP_ROOT_params.h'
514 INCLUDE 'KPP_ROOT_global.h'
515 C --- NCOEFF, JCOEFF useful for derivatives w.r.t. rate coefficients
516 INTEGER NCOEFF, JCOEFF(NREACT)
517 COMMON /DDMRCOEFF/ NCOEFF, JCOEFF
518
519 KPP_REAL T, Told
520 KPP_REAL Y(NVAR), P(NVAR*NSENSIT)
521 Told = TIME
522 TIME = T
523 CALL Update_SUN()
524 CALL Update_RCONST()
525 C
526 IF (DDMTYPE .EQ. 0) THEN
527 C --- Note: the values below are for sensitivities w.r.t. initial values;
528 C --- they may have to be changed for other applications
529 DO j=1,NSENSIT
530 DO i=1,NVAR
531 P(i+NVAR*(j-1)) = 0.0D0
532 END DO
533 END DO
534 ELSE
535 C --- Example: the call below is for sensitivities w.r.t. rate coefficients;
536 C --- JCOEFF(1:NSENSIT) are the indices of the NSENSIT rate coefficients
537 C --- w.r.t. which one differentiates
538 CALL dFun_dRcoeff( Y, FIX, NCOEFF, JCOEFF, P )
539 END IF
540 TIME = Told
541 RETURN
542 END
543
544 SUBROUTINE JAC_CHEM(N, T, Y, J)
545 INCLUDE 'KPP_ROOT_params.h'
546 INCLUDE 'KPP_ROOT_global.h'
547 INTEGER N
548 KPP_REAL Told, T
549 KPP_REAL Y(NVAR), J(LU_NONZERO)
550 Told = TIME
551 TIME = T
552 CALL Update_SUN()
553 CALL Update_RCONST()
554 CALL Jac_SP( Y, FIX, RCONST, J )
555 TIME = Told
556 RETURN
557 END
558
559
560 SUBROUTINE DJACDPAR(N, NSENSIT, T, Y, U, P)
561 C --- Computes the partial derivatives of JAC w.r.t. parameters times user vector U
562 INCLUDE 'KPP_ROOT_params.h'
563 INCLUDE 'KPP_ROOT_global.h'
564 C --- NCOEFF, JCOEFF useful for derivatives w.r.t. rate coefficients
565 INTEGER N
566 INTEGER NCOEFF, JCOEFF(NREACT)
567 COMMON /DDMRCOEFF/ NCOEFF, JCOEFF
568
569 KPP_REAL T, Told
570 KPP_REAL Y(NVAR), U(NVAR)
571 KPP_REAL P(NVAR*NSENSIT)
572 Told = TIME
573 TIME = T
574 CALL Update_SUN()
575 CALL Update_RCONST()
576 C
577 IF (DDMTYPE .EQ. 0) THEN
578 C --- Note: the values below are for sensitivities w.r.t. initial values;
579 C --- they may have to be changed for other applications
580 DO j=1,NSENSIT
581 DO i=1,NVAR
582 P(i+NVAR*(j-1)) = 0.0D0
583 END DO
584 END DO
585 ELSE
586 C --- Example: the call below is for sensitivities w.r.t. rate coefficients;
587 C --- JCOEFF(1:NSENSIT) are the indices of the NSENSIT rate coefficients
588 C --- w.r.t. which one differentiates
589 CALL dJac_dRcoeff( Y, FIX, U, NCOEFF, JCOEFF, P )
590 END IF
591 TIME = Told
592 RETURN
593 END
594
595
596 SUBROUTINE HESS_CHEM(N, T, Y, HESS)
597 INCLUDE 'KPP_ROOT_params.h'
598 INCLUDE 'KPP_ROOT_global.h'
599 INTEGER N
600 KPP_REAL Told, T
601 KPP_REAL Y(NVAR), HESS(NHESS)
602 Told = TIME
603 TIME = T
604 CALL Update_SUN()
605 CALL Update_RCONST()
606 CALL Hessian( Y, FIX, RCONST, HESS )
607 TIME = Told
608 RETURN
609 END
610
611
612
613
614
615
A coupled weather-fire forecasting model built on top of Weather Research and Forecasting (WRF).