1 SUBROUTINE INTEGRATE
( NSENSIT
, Y
, TIN
, TOUT
)
3 INCLUDE
'KPP_ROOT_params.h'
4 INCLUDE
'KPP_ROOT_global.h'
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
19 EXTERNAL FUNC_CHEM
, JAC_CHEM
23 CALL ROS3_DDM
(NVAR
,NSENSIT
,TIN
,TOUT
,STEPMIN
,STEPMAX
,
24 + STEPMIN
,Y
,ATOL
,RTOL
,
25 + Info
,FUNC_CHEM
,JAC_CHEM
)
31 SUBROUTINE ROS3_DDM
(N
,NSENSIT
,T
,Tnext
,Hmin
,Hmax
,Hstart
,
33 + Info
,FUNC_CHEM
,JAC_CHEM
)
35 INCLUDE
'KPP_ROOT_params.h'
36 INCLUDE
'KPP_ROOT_global.h'
37 INCLUDE
'KPP_ROOT_sparse.h'
39 C L-stable Rosenbrock 3(2), with
40 C strongly A-stable embedded formula for error control.
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
47 C All the arguments aggree with the KPP syntax.
50 C y = Vector of: (1:NVAR) concentrations, followed by
51 C (1:NVAR) sensitivities w.r.t. first parameter, 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
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
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.
77 C Adrian Sandu, April 1996
78 C The Center for Global and Regional Environmental Research
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
)
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
97 LOGICAL IsReject
,Autonomous
98 EXTERNAL FUNC_CHEM
, JAC_CHEM
, HESS_CHEM
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
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
118 g1
= .43586652150845899941601945119356d
+00
119 g2
= .24291996454816804366592249683314d
+00
120 g3
= .21851380027664058511513169485832d
+01
123 c Initialization of counters, etc.
124 Autonomous
= Info
(1) .EQ
. 1
126 dround
= DSQRT
(uround
)
127 IF (Hmax
.le
.0.D0
) THEN
130 H
= DMAX1
(1.d
-8, Hstart
)
139 C === Starting the time loop ===
142 C ====== Initial Function, Jacobian, and Hessian values ===============
143 CALL FUNC_CHEM
(NVAR
, T
, y
, Fv
)
145 CALL JAC_CHEM
(NVAR
, T
, y
, JAC
)
147 CALL HESS_CHEM
( NVAR
, T
, y
, HESS
)
148 IF (DDMTYPE
.EQ
. 1) THEN
149 CALL DFUNDPAR
(NVAR
, NSENSIT
, T
, y
, DFDP
)
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
)
158 DFDT
(j
) = ( K2
(j
)-Fv
(j
) )/tau
160 CALL JAC_CHEM
(NVAR
, T
+tau
, y
, AJAC
)
161 DO 30 j
= 1,LU_NONZERO
162 DJDT
(j
) = ( AJAC
(j
)-JAC
(j
) )/tau
165 CALL Jac_SP_Vec
(DJDT
,y
(i*NVAR
+1),DFDT
(i*NVAR
+1))
171 IF ( Tplus
.gt
. Tnext
) then
176 gHinv
= 1.0d0
/(gamma*H
)
181 AJAC
(LU_DIAG
(j
)) = AJAC
(LU_DIAG
(j
)) + gHinv
183 CALL KppDecomp
(AJAC
, ier
)
190 PRINT
*,'IER <> 0, H=',H
197 C ------------------------------- STAGE 1 --------------------------------------
201 IF (.NOT
.Autonomous
) THEN
204 K1
(j
) = K1
(j
) + x1*DFDT
(j
)
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
)
212 C --- End of derivative w.r.t. parameters
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
)
218 K1
(i*NVAR
+j
) = K1
(i*NVAR
+j
) + Hv
(j
)
220 IF (.NOT
. Autonomous
) THEN
222 K1
(i*NVAR
+j
) = K1
(i*NVAR
+j
) + x1*DFDT
(i*NVAR
+j
)
225 C --- If derivative w.r.t. parameters
226 IF (DDMTYPE
.EQ
. 1) THEN
228 K1
(i*NVAR
+j
) = K1
(i*NVAR
+j
) + DFDP
((i
-1)*NVAR
+j
)
229 & + DJDP
((i
-1)*NVAR
+j
)
232 C --- End of derivative w.r.t. parameters
233 CALL KppSolve
(AJAC
, K1
(i*NVAR
+1))
236 C ------------------------------- STAGE 2 --------------------------------------
237 DO 120 j
= 1,NVAR*
(NSENSIT
+1)
238 ynew
(j
) = y
(j
) + a21*K1
(j
)
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
)
247 K2
(j
) = Fv
(j
) + x1*K1
(j
)
249 IF (.NOT
.Autonomous
) THEN
252 K2
(j
) = K2
(j
) + x2*DFDT
(j
)
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
)
260 C --- End of derivative w.r.t. parameters
262 CALL JAC_CHEM
(NVAR
, T
+alpha2*H
, ynew
, JAC
)
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
)
268 K2
(i*NVAR
+j
) = K2
(i*NVAR
+j
) + x1*K1
(i*NVAR
+j
)
271 IF (.NOT
. Autonomous
) THEN
273 K2
(i*NVAR
+j
) = K2
(i*NVAR
+j
) + x2*DFDT
(i*NVAR
+j
)
276 C --- If derivative w.r.t. parameters
277 IF (DDMTYPE
.EQ
. 1) THEN
279 K2
(i*NVAR
+j
) = K2
(i*NVAR
+j
) + DFDP
((i
-1)*NVAR
+j
)
280 & + DJDP
((i
-1)*NVAR
+j
)
283 C --- End of derivative w.r.t. parameters
284 CALL KppSolve
(AJAC
, K2
(i*NVAR
+1))
287 C ------------------------------- STAGE 3 --------------------------------------
291 K3
(j
) = Fv
(j
) + x1*K1
(j
) + x2*K2
(j
)
293 IF (.NOT
.Autonomous
) THEN
296 K3
(j
) = K3
(j
) + x3*DFDT
(j
)
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
)
304 C --- End of derivative w.r.t. parameters
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
)
310 K3
(i*NVAR
+j
) = K3
(i*NVAR
+j
) +x1*K1
(i*NVAR
+j
)
311 & + x2*K2
(i*NVAR
+j
) + Hv
(j
)
313 IF (.NOT
. Autonomous
) THEN
315 K3
(i*NVAR
+j
) = K3
(i*NVAR
+j
) + x3*DFDT
(i*NVAR
+j
)
318 C --- If derivative w.r.t. parameters
319 IF (DDMTYPE
.EQ
. 1) THEN
321 K3
(i*NVAR
+j
) = K3
(i*NVAR
+j
) + DFDP
((i
-1)*NVAR
+j
)
322 & + DJDP
((i
-1)*NVAR
+j
)
325 C --- End of derivative w.r.t. parameters
326 CALL KppSolve
(AJAC
, K3
(i*NVAR
+1))
329 C ------------------------------ The Solution ---
331 DO 230 j
= 1,NVAR*
(NSENSIT
+1)
332 ynew
(j
) = y
(j
) + b1*K1
(j
) + b2*K2
(j
) + b3*K3
(j
)
336 C ====== Error estimation ========
340 ytol
= AbsTol
(i
) + RelTol
(i
)*DABS
(ynew
(i
))
341 ERR
=ERR
+((d1*K1
(i
)+d2*K2
(i
)+d3*K3
(i
))/ytol
)**2
343 ERR
= DMAX1
( uround
, DSQRT
( ERR
/NVAR
) )
345 C ======= Choose the stepsize ===============================
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
))
351 C ======= Rejected/Accepted Step ============================
353 IF ( (ERR
.gt
.1).and
.(H
.gt
.Hmin
) ) THEN
359 DO 250 j
=1,NVAR*
(NSENSIT
+1)
363 IF (.NOT
.IsReject
) THEN
364 H
= Hnew
! Do not increase stepsize
IF previos step was rejected
371 C ======= End of the time loop ===============================
372 IF ( T
.lt
. Tnext
) GO TO 10
375 C ======= Output Information =================================
387 SUBROUTINE FUNC_CHEM
(N
, T
, Y
, P
)
388 INCLUDE
'KPP_ROOT_params.h'
389 INCLUDE
'KPP_ROOT_global.h'
392 KPP_REAL Y
(NVAR
), P
(NVAR
)
397 CALL Fun
( Y
, FIX
, RCONST
, P
)
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
409 INTEGER NCOEFF
, JCOEFF
(NREACT
)
410 COMMON /DDMRCOEFF
/ NCOEFF
, JCOEFF
413 KPP_REAL Y
(NVAR
), P
(NVAR*NSENSIT
)
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
424 P
(i
+NVAR*
(j
-1)) = 0.0D0
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
)
437 SUBROUTINE JAC_CHEM
(N
, T
, Y
, J
)
438 INCLUDE
'KPP_ROOT_params.h'
439 INCLUDE
'KPP_ROOT_global.h'
442 KPP_REAL Y
(NVAR
), J
(LU_NONZERO
)
447 CALL Jac_SP
( Y
, FIX
, RCONST
, J
)
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
459 INTEGER NCOEFF
, JCOEFF
(NREACT
)
460 COMMON /DDMRCOEFF
/ NCOEFF
, JCOEFF
463 KPP_REAL Y
(NVAR
), U
(NVAR
)
464 KPP_REAL P
(NVAR*NSENSIT
)
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
475 P
(i
+NVAR*
(j
-1)) = 0.0D0
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
)
489 SUBROUTINE HESS_CHEM
(N
, T
, Y
, HESS
)
490 INCLUDE
'KPP_ROOT_params.h'
491 INCLUDE
'KPP_ROOT_global.h'
494 KPP_REAL Y
(NVAR
), HESS
(NHESS
)
499 CALL Hessian
( Y
, FIX
, RCONST
, HESS
)