chem/KPP/kpp/kpp-2.1/int/oldies/rodas3_ddm.f

   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