source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int.modified_WCOPY/oldies/rodas3_ddm.f @ 2696

Last change on this file since 2696 was 2696, checked in by kanani, 6 years ago

Merge of branch palm4u into trunk

File size: 19.0 KB
Line 
1      SUBROUTINE INTEGRATE( NSENSIT, Y, TIN, TOUT )
2 
3      INCLUDE 'KPP_ROOT_params.h'
4      INCLUDE 'KPP_ROOT_global.h'
5
6C TIN - Start Time
7      KPP_REAL TIN
8C TOUT - End Time
9      KPP_REAL TOUT
10C Y - Concentrations and Sensitivities
11      KPP_REAL Y(NVAR*(NSENSIT+1))
12C ---  Note: Y contains: (1:NVAR) concentrations, followed by
13C ---                   (1:NVAR) sensitivities w.r.t. first parameter, followed by
14C ---                   etc.,  followed by         
15C ---                   (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'                                                                                                 
42C
43C       Stiffly accurate Rosenbrock 3(2), with
44C       stiffly accurate embedded formula for error control. 
45C
46C Direct decoupled computation of sensitivities.
47C The global variable DDMTYPE distinguishes between:
48C      DDMTYPE = 0 : sensitivities w.r.t. initial values
49C      DDMTYPE = 1 : sensitivities w.r.t. parameters
50C
51C  INPUT ARGUMENTS:
52C     y = Vector of:   (1:NVAR) concentrations, followed by
53C                      (1:NVAR) sensitivities w.r.t. first parameter, followed by
54C                       etc.,  followed by         
55C                      (1:NVAR) sensitivities w.r.t. NSENSIT's parameter
56C         (y contains initial values at input, final values at output)
57C     [T, Tnext] = the integration interval
58C     Hmin, Hmax = lower and upper bounds for the selected step-size.
59C          Note that for Step = Hmin the current computed
60C          solution is unconditionally accepted by the error
61C          control mechanism.
62C     AbsTol, RelTol = (NVAR) dimensional vectors of
63C          componentwise absolute and relative tolerances.
64C     FUNC_CHEM = name of routine of derivatives. KPP syntax.
65C          See the header below.
66C     JAC_CHEM = name of routine that computes the Jacobian, in
67C          sparse format. KPP syntax. See the header below.
68C     Info(1) = 1  for  Autonomous   system
69C             = 0  for nonAutonomous system
70C
71C  OUTPUT ARGUMENTS:
72C     y = the values of concentrations and sensitivities at Tend.
73C     T = equals TENDon output.
74C     Info(2) = # of FUNC_CHEM CALLs.
75C     Info(3) = # of JAC_CHEM CALLs.
76C     Info(4) = # of accepted steps.
77C     Info(5) = # of rejected steps.
78C 
79C    Adrian Sandu, December 2001
80C 
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
107C   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
143c     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 
158C === 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 
167C              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 
175C              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 
194C              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)
204C 
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                                                                                                                           
215C ------------ 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)
226C               --- 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
230C               --- 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
243C               --- 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
250C               --- End of derivative w.r.t. parameters         
251          CALL KppSolve (AJAC, K1(i*NVAR+1))
252 100   CONTINUE
253 
254C ----------- STAGE 2 -------------------------
255C Note: uses the same function values as Stage 1
256C       DO 110 j = 1,NVAR*(NSENSIT+1)
257C         ynew(j) = y(j) + a21*K1(j)
258C 110   CONTINUE
259C       CALL FUNC_CHEM(NVAR, T+alpha2*H, ynew, Fv)
260C       IF (DDMTYPE .EQ. 1) THEN
261C         CALL DFUNDPAR(NVAR, NSENSIT, T+alpha2*H, ynew, DFDP)
262C       END IF   
263C       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)
275C               --- 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
279C               --- 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
295C               --- 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
302C               --- End of derivative w.r.t. parameters         
303          CALL KppSolve (AJAC, K2(i*NVAR+1))
304 160   CONTINUE
305 
306 
307C ------------ 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)
328C               --- 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
332C               --- 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
348C               --- 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         
355C               --- End of derivative w.r.t. parameters
356          CALL KppSolve (AJAC, K3(i*NVAR+1))
357 220   CONTINUE
358 
359C ------------ 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)
381C               --- 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
385C               --- 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
401C --- 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
412C ---- The Solution ---
413       DO 280 j = 1,NVAR*(NSENSIT+1)
414C         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 
419C ====== 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)))
424C          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 
430C ======= 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 
436C ======= 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 
454C ======= End of the time loop ===============================
455      IF ( T .lt. Tnext ) GO TO 10
456 
457 
458 
459C ======= 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)
488C ---  Computes the partial derivatives of FUNC_CHEM w.r.t. parameters 
489      INCLUDE 'KPP_ROOT_params.h'
490      INCLUDE 'KPP_ROOT_global.h'
491C ---  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()
502C
503      IF (DDMTYPE .EQ. 0) THEN
504C ---  Note: the values below are for sensitivities w.r.t. initial values;
505C ---       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     
512C ---  Example: the call below is for sensitivities w.r.t. rate coefficients;
513C ---       JCOEFF(1:NSENSIT) are the indices of the NSENSIT rate coefficients
514C ---       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)
538C ---  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'
541C ---  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()
553C
554      IF (DDMTYPE .EQ. 0) THEN
555C ---  Note: the values below are for sensitivities w.r.t. initial values;
556C ---       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     
563C ---  Example: the call below is for sensitivities w.r.t. rate coefficients;
564C ---       JCOEFF(1:NSENSIT) are the indices of the NSENSIT rate coefficients
565C ---       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                                                                                                                           
Note: See TracBrowser for help on using the repository browser.