chem/KPP/kpp/kpp-2.1/int/kpp_sdirk.f90

   1 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
   2 !  SDIRK - Singly-Diagonally-Implicit Runge-Kutta method                  !
   3 !                       (L-stable, 5 stages, order 4)                     !
   4 !  By default the code employs the KPP sparse linear algebra routines     !
   5 !  Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK)        !
   6 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
   7 ! A. Sandu - version of July 10, 2005
   8
   9 MODULE KPP_ROOT_Integrator
  10
  11   USE KPP_ROOT_Precision
  12   USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME
  13   USE KPP_ROOT_Parameters, ONLY: NVAR, NSPEC, NFIX, LU_NONZERO
  14   USE KPP_ROOT_JacobianSP, ONLY: LU_DIAG
  15   USE KPP_ROOT_LinearAlgebra, ONLY: KppDecomp, KppSolve, &
  16                Set2zero, WLAMCH, WAXPY, WCOPY
  17
  18   IMPLICIT NONE
  19   PUBLIC
  20   SAVE
  21
  22   !~~~>  Statistics on the work performed by the SDIRK method
  23   INTEGER :: Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng
  24   INTEGER, PARAMETER :: ifun=1, ijac=2, istp=3, iacc=4,  &
  25     irej=5, idec=6, isol=7, isng=8, itexit=1, ihexit=2
  26   !  SDIRK method coefficients
  27   KPP_REAL :: rkAlpha(5,4), rkBeta(5,4), rkD(4,5),  &
  28                    rkGamma, rkA(5,5), rkB(5), rkC(5)
  29
  30   ! description of the error numbers IERR
  31   CHARACTER(LEN=50), PARAMETER, DIMENSION(-8:1) :: IERR_NAMES = (/ &
  32     'Matrix is repeatedly singular                     ', & ! -8
  33     'Step size too small: T + 10*H = T or H < Roundoff ', & ! -7
  34     'No of steps exceeds maximum bound                 ', & ! -6
  35     'Improper tolerance values                         ', & ! -5
  36     'FacMin/FacMax/FacRej must be positive             ', & ! -4
  37     'Hmin/Hmax/Hstart must be positive                 ', & ! -3
  38     'Improper value for maximal no of Newton iterations', & ! -2
  39     'Improper value for maximal no of steps            ', & ! -1
  40     '                                                  ', & !  0 (not used)
  41     'Success                                           ' /) !  1
  42
  43 CONTAINS
  44
  45 SUBROUTINE INTEGRATE( TIN, TOUT, &
  46   ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U, IERR_U )
  47
  48    USE KPP_ROOT_Parameters
  49    USE KPP_ROOT_Global
  50    IMPLICIT NONE
  51
  52    KPP_REAL, INTENT(IN) :: TIN  ! Start Time
  53    KPP_REAL, INTENT(IN) :: TOUT ! End Time
  54    ! Optional input parameters and statistics
  55    INTEGER,  INTENT(IN),  OPTIONAL :: ICNTRL_U(20)
  56    KPP_REAL, INTENT(IN),  OPTIONAL :: RCNTRL_U(20)
  57    INTEGER,  INTENT(OUT), OPTIONAL :: ISTATUS_U(20)
  58    KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20)
  59    INTEGER,  INTENT(OUT), OPTIONAL :: IERR_U
  60
  61    !INTEGER, SAVE :: Ntotal = 0
  62    KPP_REAL :: RCNTRL(20), RSTATUS(20)
  63    INTEGER  :: ICNTRL(20), ISTATUS(20), IERR
  64
  65    ICNTRL(:)  = 0
  66    RCNTRL(:)  = 0.0_dp
  67    ISTATUS(:) = 0
  68    RSTATUS(:) = 0.0_dp
  69
  70    ! If optional parameters are given, and if they are >0,
  71    ! then they overwrite default settings.
  72    IF (PRESENT(ICNTRL_U)) THEN
  73      WHERE(ICNTRL_U(:) > 0) ICNTRL(:) = ICNTRL_U(:)
  74    END IF
  75    IF (PRESENT(RCNTRL_U)) THEN
  76      WHERE(RCNTRL_U(:) > 0) RCNTRL(:) = RCNTRL_U(:)
  77    END IF
  78
  79    CALL SDIRK( NVAR,TIN,TOUT,VAR,RTOL,ATOL,          &
  80                RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR )
  81
  82 ! mz_rs_20050716: IERR and ISTATUS(istp) are returned to the user who then
  83 ! decides what to do about it, i.e. either stop the run or ignore it.
  84 !!$   IF (IERR < 0) THEN
  85 !!$        PRINT *,'SDIRK: Unsuccessful exit at T=',TIN,' (IERR=',IERR,')'
  86 !!$   ENDIF
  87 !!$   Ntotal = Ntotal + Nstp
  88 !!$   PRINT*,'NSTEPS=',Nstp, '(',Ntotal,')'
  89
  90    STEPMIN = RSTATUS(ihexit) ! Save last step
  91
  92    ! if optional parameters are given for output they to return information
  93    IF (PRESENT(ISTATUS_U)) ISTATUS_U(:) = ISTATUS(1:20)
  94    IF (PRESENT(RSTATUS_U)) RSTATUS_U(:) = RSTATUS(1:20)
  95    IF (PRESENT(IERR_U))    IERR_U       = IERR
  96
  97    END SUBROUTINE INTEGRATE
  98
  99 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 100       SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol,     &
 101                        RCNTRL, ICNTRL, RSTATUS, ISTATUS, IDID)
 102 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 103 !
 104 !    Solves the system y'=F(t,y) using a Singly-Diagonally-Implicit
 105 !    Runge-Kutta (SDIRK) method.
 106 !
 107 !    For details on SDIRK methods and their implementation consult:
 108 !      E. Hairer and G. Wanner
 109 !      "Solving ODEs II. Stiff and differential-algebraic problems".
 110 !      Springer series in computational mathematics, Springer-Verlag, 1996.
 111 !    This code is based on the SDIRK4 routine in the above book.
 112 !
 113 !    (C)  Adrian Sandu, July 2005
 114 !    Virginia Polytechnic Institute and State University
 115 !    Contact: sandu@cs.vt.edu
 116 !    This implementation is part of KPP - the Kinetic PreProcessor
 117 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 118 !
 119 !~~~>   INPUT ARGUMENTS:
 120 !
 121 !-     Y(NVAR)    = vector of initial conditions (at T=Tinitial)
 122 !-    [Tinitial,Tfinal]  = time range of integration
 123 !     (if Tinitial>Tfinal the integration is performed backwards in time)
 124 !-    RelTol, AbsTol = user precribed accuracy
 125 !- SUBROUTINE ode_Fun( T, Y, Ydot ) = ODE function,
 126 !                       returns Ydot = Y' = F(T,Y)
 127 !- SUBROUTINE ode_Fun( T, Y, Ydot ) = Jacobian of the ODE function,
 128 !                       returns Jcb = dF/dY
 129 !-    ICNTRL(1:20)    = integer inputs parameters
 130 !-    RCNTRL(1:20)    = real inputs parameters
 131 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 132 !
 133 !~~~>     OUTPUT ARGUMENTS:
 134 !
 135 !-    Y(NVAR)         -> vector of final states (at T->Tfinal)
 136 !-    ISTATUS(1:20)   -> integer output parameters
 137 !-    RSTATUS(1:20)   -> real output parameters
 138 !-    IDID            -> job status upon return
 139 !                        success (positive value) or
 140 !                        failure (negative value)
 141 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 142 !
 143 !~~~>     INPUT PARAMETERS:
 144 !
 145 !    Note: For input parameters equal to zero the default values of the
 146 !       corresponding variables are used.
 147 !
 148 !    Note: For input parameters equal to zero the default values of the
 149 !          corresponding variables are used.
 150 !~~~>
 151 !    ICNTRL(1) = not used
 152 !
 153 !    ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors
 154 !              = 1: AbsTol, RelTol are scalars
 155 !
 156 !    ICNTRL(3) = not used
 157 !
 158 !    ICNTRL(4)  -> maximum number of integration steps
 159 !        For ICNTRL(4)=0 the default value of 100000 is used
 160 !
 161 !    ICNTRL(5)  -> maximum number of Newton iterations
 162 !        For ICNTRL(5)=0 the default value of 8 is used
 163 !
 164 !    ICNTRL(6)  -> starting values of Newton iterations:
 165 !        ICNTRL(6)=0 : starting values are interpolated (the default)
 166 !        ICNTRL(6)=1 : starting values are zero
 167 !
 168 !~~~>  Real parameters
 169 !
 170 !    RCNTRL(1)  -> Hmin, lower bound for the integration step size
 171 !                  It is strongly recommended to keep Hmin = ZERO
 172 !    RCNTRL(2)  -> Hmax, upper bound for the integration step size
 173 !    RCNTRL(3)  -> Hstart, starting value for the integration step size
 174 !
 175 !    RCNTRL(4)  -> FacMin, lower bound on step decrease factor (default=0.2)
 176 !    RCNTRL(5)  -> FacMax, upper bound on step increase factor (default=6)
 177 !    RCNTRL(6)  -> FacRej, step decrease factor after multiple rejections
 178 !                 (default=0.1)
 179 !    RCNTRL(7)  -> FacSafe, by which the new step is slightly smaller
 180 !                  than the predicted value  (default=0.9)
 181 !    RCNTRL(8)  -> ThetaMin. If Newton convergence rate smaller
 182 !                  than ThetaMin the Jacobian is not recomputed;
 183 !                  (default=0.001)
 184 !    RCNTRL(9)  -> NewtonTol, stopping criterion for Newton's method
 185 !                  (default=0.03)
 186 !    RCNTRL(10) -> Qmin
 187 !    RCNTRL(11) -> Qmax. If Qmin < Hnew/Hold < Qmax, then the
 188 !                  step size is kept constant and the LU factorization
 189 !                  reused (default Qmin=1, Qmax=1.2)
 190 !
 191 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 192 !
 193 !~~~>     OUTPUT PARAMETERS:
 194 !
 195 !    Note: each call to Rosenbrock adds the current no. of fcn calls
 196 !      to previous value of ISTATUS(1), and similar for the other params.
 197 !      Set ISTATUS(1:10) = 0 before call to avoid this accumulation.
 198 !
 199 !    ISTATUS(1) = No. of function calls
 200 !    ISTATUS(2) = No. of jacobian calls
 201 !    ISTATUS(3) = No. of steps
 202 !    ISTATUS(4) = No. of accepted steps
 203 !    ISTATUS(5) = No. of rejected steps (except at the beginning)
 204 !    ISTATUS(6) = No. of LU decompositions
 205 !    ISTATUS(7) = No. of forward/backward substitutions
 206 !    ISTATUS(8) = No. of singular matrix decompositions
 207 !
 208 !    RSTATUS(1)  -> Texit, the time corresponding to the
 209 !                   computed Y upon return
 210 !    RSTATUS(2)  -> Hexit, last predicted step before exit
 211 !    For multiple restarts, use Hexit as Hstart in the following run
 212 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 213       IMPLICIT NONE
 214
 215 ! Arguments
 216       INTEGER, INTENT(IN)          :: N, ICNTRL(20)
 217       KPP_REAL, INTENT(IN)    :: Tinitial, Tfinal,    &
 218                     RelTol(NVAR), AbsTol(NVAR), RCNTRL(20)
 219       KPP_REAL, INTENT(INOUT) :: Y(NVAR)
 220       INTEGER, INTENT(OUT)         :: IDID
 221       INTEGER, INTENT(INOUT)       :: ISTATUS(20)
 222       KPP_REAL, INTENT(OUT)   :: RSTATUS(20)
 223
 224 ! Local variables
 225       KPP_REAL :: Hmin, Hmax, Hstart, Roundoff,    &
 226                        FacMin, Facmax, FacSafe, FacRej, &
 227                        ThetaMin, NewtonTol, Qmin, Qmax, &
 228                        Texit, Hexit
 229       INTEGER :: ITOL, NewtonMaxit, Max_no_steps, i
 230       KPP_REAL, PARAMETER :: ZERO = 0.0d0
 231
 232 !~~~>  Initialize statistics
 233    Nfun = ISTATUS(ifun)
 234    Njac = ISTATUS(ijac)
 235    Nstp = ISTATUS(istp)
 236    Nacc = ISTATUS(iacc)
 237    Nrej = ISTATUS(irej)
 238    Ndec = ISTATUS(idec)
 239    Nsol = ISTATUS(isol)
 240    Nsng = ISTATUS(isng)
 241 !   Nfun=0; Njac=0; Nstp=0; Nacc=0
 242 !   Nrej=0; Ndec=0; Nsol=0; Nsng=0
 243
 244    IDID = 0
 245
 246 !~~~>  For Scalar tolerances (ICNTRL(2).NE.0)  the code uses AbsTol(1) and RelTol(1)
 247 !   For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:NVAR) and RelTol(1:NVAR)
 248    IF (ICNTRL(2) == 0) THEN
 249       ITOL = 1
 250    ELSE
 251       ITOL = 0
 252    END IF
 253
 254 !~~~>   The maximum number of time steps admitted
 255    IF (ICNTRL(3) == 0) THEN
 256       Max_no_steps = 100000
 257    ELSEIF (Max_no_steps > 0) THEN
 258       Max_no_steps=ICNTRL(3)
 259    ELSE
 260       PRINT * ,'User-selected ICNTRL(3)=',ICNTRL(3)
 261       CALL SDIRK_ErrorMsg(-1,Tinitial,ZERO,IDID)
 262    END IF
 263
 264
 265 !~~~> The maximum number of Newton iterations admitted
 266    IF(ICNTRL(4) == 0)THEN
 267       NewtonMaxit=8
 268    ELSE
 269       NewtonMaxit=ICNTRL(4)
 270       IF(NewtonMaxit <= 0)THEN
 271           PRINT * ,'User-selected ICNTRL(4)=',ICNTRL(4)
 272           CALL SDIRK_ErrorMsg(-2,Tinitial,ZERO,IDID)
 273       END IF
 274    END IF
 275
 276 !~~~>  Unit roundoff (1+Roundoff>1)
 277       Roundoff = WLAMCH('E')
 278
 279 !~~~>  Lower bound on the step size: (positive value)
 280    IF (RCNTRL(1) == ZERO) THEN
 281       Hmin = ZERO
 282    ELSEIF (RCNTRL(1) > ZERO) THEN
 283       Hmin = RCNTRL(1)
 284    ELSE
 285       PRINT * , 'User-selected RCNTRL(1)=', RCNTRL(1)
 286       CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,IDID)
 287    END IF
 288
 289 !~~~>  Upper bound on the step size: (positive value)
 290    IF (RCNTRL(2) == ZERO) THEN
 291       Hmax = ABS(Tfinal-Tinitial)
 292    ELSEIF (RCNTRL(2) > ZERO) THEN
 293       Hmax = MIN(ABS(RCNTRL(2)),ABS(Tfinal-Tinitial))
 294    ELSE
 295       PRINT * , 'User-selected RCNTRL(2)=', RCNTRL(2)
 296       CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,IDID)
 297    END IF
 298
 299 !~~~>  Starting step size: (positive value)
 300    IF (RCNTRL(3) == ZERO) THEN
 301       Hstart = MAX(Hmin,Roundoff)
 302    ELSEIF (RCNTRL(3) > ZERO) THEN
 303       Hstart = MIN(ABS(RCNTRL(3)),ABS(Tfinal-Tinitial))
 304    ELSE
 305       PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3)
 306       CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,IDID)
 307    END IF
 308
 309 !~~~>  Step size can be changed s.t.  FacMin < Hnew/Hexit < FacMax
 310    IF (RCNTRL(4) == ZERO) THEN
 311       FacMin = 0.2_dp
 312    ELSEIF (RCNTRL(4) > ZERO) THEN
 313       FacMin = RCNTRL(4)
 314    ELSE
 315       PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4)
 316       CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,IDID)
 317    END IF
 318    IF (RCNTRL(5) == ZERO) THEN
 319       FacMax = 10.0_dp
 320    ELSEIF (RCNTRL(5) > ZERO) THEN
 321       FacMax = RCNTRL(5)
 322    ELSE
 323       PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5)
 324       CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,IDID)
 325    END IF
 326 !~~~>   FacRej: Factor to decrease step after 2 succesive rejections
 327    IF (RCNTRL(6) == ZERO) THEN
 328       FacRej = 0.1_dp
 329    ELSEIF (RCNTRL(6) > ZERO) THEN
 330       FacRej = RCNTRL(6)
 331    ELSE
 332       PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6)
 333       CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,IDID)
 334    END IF
 335 !~~~>   FacSafe: Safety Factor in the computation of new step size
 336    IF (RCNTRL(7) == ZERO) THEN
 337       FacSafe = 0.9_dp
 338    ELSEIF (RCNTRL(7) > ZERO) THEN
 339       FacSafe = RCNTRL(7)
 340    ELSE
 341       PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7)
 342       CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,IDID)
 343    END IF
 344
 345 !~~~> DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED;
 346       IF(RCNTRL(8) == 0.D0)THEN
 347          ThetaMin = 1.0d-3
 348       ELSE
 349          ThetaMin = RCNTRL(8)
 350       END IF
 351
 352 !~~~> STOPPING CRITERION FOR NEWTON'S METHOD
 353       IF(RCNTRL(9) == 0.0d0)THEN
 354          NewtonTol = 3.0d-2
 355       ELSE
 356          NewtonTol =RCNTRL(9)
 357       END IF
 358
 359 !~~~> Qmin AND Qmax: IF Qmin < Hnew/Hold < Qmax, STEP SIZE = CONST.
 360       IF(RCNTRL(10) == 0.D0)THEN
 361          Qmin=1.D0
 362       ELSE
 363          Qmin=RCNTRL(10)
 364       END IF
 365       IF(RCNTRL(11) == 0.D0)THEN
 366          Qmax=1.2D0
 367       ELSE
 368          Qmax=RCNTRL(11)
 369       END IF
 370
 371 !~~~>  Check if tolerances are reasonable
 372       IF (ITOL == 0) THEN
 373          IF (AbsTol(1) <= 0.D0.OR.RelTol(1) <= 10.D0*Roundoff) THEN
 374             PRINT * , ' Scalar AbsTol = ',AbsTol(1)
 375             PRINT * , ' Scalar RelTol = ',RelTol(1)
 376             CALL SDIRK_ErrorMsg(-5,Tinitial,ZERO,IDID)
 377          END IF
 378       ELSE
 379          DO i=1,N
 380             IF (AbsTol(i) <= 0.D0.OR.RelTol(i) <= 10.D0*Roundoff) THEN
 381               PRINT * , ' AbsTol(',i,') = ',AbsTol(i)
 382               PRINT * , ' RelTol(',i,') = ',RelTol(i)
 383               CALL SDIRK_ErrorMsg(-5,Tinitial,ZERO,IDID)
 384             END IF
 385          END DO
 386       END IF
 387
 388     IF (IDID < 0) RETURN
 389
 390
 391 !~~~> CALL TO CORE INTEGRATOR
 392     CALL SDIRK_Integrator( N,Tinitial,Tfinal,Y,Hmin,Hmax,Hstart, &
 393           RelTol,AbsTol,ITOL, Max_no_steps, NewtonMaxit,         &
 394           Roundoff, FacMin, FacMax, FacRej, FacSafe, ThetaMin,   &
 395           NewtonTol, Qmin, Qmax, Hexit, Texit, IDID )
 396
 397 !~~~>  Collect run statistics
 398    ISTATUS(ifun) = Nfun
 399    ISTATUS(ijac) = Njac
 400    ISTATUS(istp) = Nstp
 401    ISTATUS(iacc) = Nacc
 402    ISTATUS(irej) = Nrej
 403    ISTATUS(idec) = Ndec
 404    ISTATUS(isol) = Nsol
 405    ISTATUS(isng) = Nsng
 406 !~~~> Last T and H
 407    RSTATUS(itexit) = Texit
 408    RSTATUS(ihexit) = Hexit
 409
 410 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 411    CONTAINS  !  PROCEDURES INTERNAL TO SDIRK
 412 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 413
 414
 415 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 416    SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Hmin,Hmax,Hstart, &
 417                 RelTol,AbsTol,ITOL, Max_no_steps, NewtonMaxit,        &
 418                 Roundoff, FacMin, FacMax, FacRej, FacSafe, ThetaMin,  &
 419                 NewtonTol, Qmin, Qmax, Hexit, Texit, IDID )
 420 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 421 !     CORE INTEGRATOR FOR SDIRK4
 422 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 423
 424       USE KPP_ROOT_Parameters
 425       IMPLICIT NONE
 426
 427 !~~~> Arguments:
 428       INTEGER :: N
 429       KPP_REAL, INTENT(INOUT) :: Y(NVAR)
 430       KPP_REAL, INTENT(IN) :: Tinitial, Tfinal, Hmin, Hmax, Hstart, &
 431                         RelTol(NVAR), AbsTol(NVAR), Roundoff,            &
 432                         FacMin, FacMax, FacRej, FacSafe, ThetaMin,       &
 433                         NewtonTol, Qmin, Qmax
 434       KPP_REAL, INTENT(OUT) :: Hexit, Texit
 435       INTEGER, INTENT(IN)  :: ITOL, Max_no_steps, NewtonMaxit
 436       INTEGER, INTENT(OUT) :: IDID
 437
 438 !~~~> Local variables:
 439       KPP_REAL :: Z(NVAR,5), FV(NVAR,5), CONT(NVAR,4),      &
 440                        NewtonFactor(5), SCAL(NVAR), RHS(NVAR),   &
 441                        G(NVAR), Yhat(NVAR), TMP(NVAR),           &
 442                        T, H, Hold, Theta, Hratio, Hmax1, W,      &
 443                        HGammaInv, DYTH, QNEWT, ERR, Fac, Hnew,   &
 444                        Tdirection, NewtonErr, NewtonErrOld
 445       INTEGER :: i, j, IER, istage, NewtonIter, IP(NVAR)
 446       LOGICAL :: Reject, FIRST, NewtonReject, FreshJac, SkipJacUpdate, SkipLU
 447
 448 #ifdef FULL_ALGEBRA
 449       KPP_REAL FJAC(NVAR,NVAR),  E(NVAR,NVAR)
 450 #else
 451       KPP_REAL FJAC(LU_NONZERO), E(LU_NONZERO)
 452 #endif
 453       KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
 454
 455 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 456 !  INITIALISATIONS
 457 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 458       CALL SDIRK_Coefficients
 459       T = Tinitial
 460       Tdirection = SIGN(1.D0,Tfinal-Tinitial)
 461       Hmax1=MIN(ABS(Hmax),ABS(Tfinal-Tinitial))
 462       H = MAX(ABS(Hmin),ABS(Hstart))
 463       IF (ABS(H) <= 10.D0*Roundoff) H=1.0D-6
 464       H=MIN(ABS(H),Hmax1)
 465       H=SIGN(H,Tdirection)
 466       Hold=H
 467       NewtonReject=.FALSE.
 468       SkipLU =.FALSE.
 469       FreshJac = .FALSE.
 470       SkipJacUpdate = .FALSE.
 471       Reject=.FALSE.
 472       FIRST=.TRUE.
 473       NewtonFactor(1:5)=ONE
 474
 475       CALL SDIRK_ErrorScale(ITOL, AbsTol, RelTol, Y, SCAL)
 476
 477 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 478 !~~~>  Time loop begins
 479 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 480 Tloop: DO WHILE ( (Tfinal-T)*Tdirection - Roundoff > ZERO )
 481
 482
 483 !~~~>  Compute E = 1/(h*gamma)-Jac and its LU decomposition
 484       IF ( SkipLU ) THEN ! This time around skip the Jac update and LU
 485          SkipLU = .FALSE.; FreshJac = .FALSE.; SkipJacUpdate = .FALSE.
 486       ELSE
 487          CALL SDIRK_PrepareMatrix ( H, T, Y, FJAC, &
 488                    FreshJac, SkipJacUpdate, E, IP, Reject, IER )
 489          IF (IER /= 0) THEN
 490              CALL SDIRK_ErrorMsg(-8,T,H,IDID); RETURN
 491          END IF
 492       END IF
 493
 494       IF (Nstp>Max_no_steps) THEN
 495              CALL SDIRK_ErrorMsg(-6,T,H,IDID); RETURN
 496       END IF
 497       IF ( (T+0.1d0*H == T) .OR. (ABS(H) <= Roundoff) ) THEN
 498              CALL SDIRK_ErrorMsg(-7,T,H,IDID); RETURN
 499       END IF
 500
 501       HGammaInv = ONE/(H*rkGamma)
 502
 503 !~~~>    NEWTON ITERATION
 504 stages:DO istage=1,5
 505
 506       NewtonFactor(istage) = MAX(NewtonFactor(istage),Roundoff)**0.8d0
 507
 508 !~~~>  STARTING VALUES FOR NEWTON ITERATION
 509       CALL Set2zero(N,G)
 510       CALL Set2zero(N,Z(1,istage))
 511       IF (istage==1) THEN
 512           IF (FIRST.OR.NewtonReject) THEN
 513               CALL Set2zero(N,Z(1,istage))
 514           ELSE
 515               W=ONE+rkGamma*H/Hold
 516               DO i=1,N
 517                 Z(i,istage)=W*(CONT(i,1)+W*(CONT(i,2)+W*(CONT(i,3)+W*CONT(i,4))))-Yhat(i)
 518               END DO
 519           END IF
 520       ELSE
 521           DO j = 1, istage-1
 522             ! Gj(:) = sum_j Beta(i,j)*Zj(:) = H * sum_j A(i,j)*Fun(Zj(:))
 523             CALL WAXPY(N,rkBeta(istage,j),Z(1,j),1,G,1)
 524             ! CALL WAXPY(N,H*rkA(istage,j),FV(1,j),1,G,1)
 525             ! Zi(:) = sum_j Alpha(i,j)*Zj(:)
 526             CALL WAXPY(N,rkAlpha(istage,j),Z(1,j),1,Z(1,istage),1)
 527           END DO
 528           IF (istage==5) CALL WCOPY(N,Z(1,istage),1,Yhat,1)  ! Yhat(:) <- Z5(:)
 529       END IF
 530
 531 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 532 !  LOOP FOR THE SIMPLIFIED NEWTON ITERATION
 533 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 534             NewtonIter=0
 535             Theta=ABS(ThetaMin)
 536             IF (Reject) Theta=2*ABS(ThetaMin)
 537             NewtonErr = 1.0e+6 ! To force-enter Newton loop
 538
 539   Newton:   DO WHILE (NewtonFactor(istage)*NewtonErr > NewtonTol)
 540
 541             IF (NewtonIter >= NewtonMaxit) THEN
 542                 H=H*0.5d0
 543                 Reject=.TRUE.
 544                 NewtonReject=.TRUE.
 545                 CYCLE Tloop
 546             END IF
 547             NewtonIter=NewtonIter+1
 548
 549 !~~~>     COMPUTE THE RIGHT-HAND SIDE
 550             TMP(1:N) = Y(1:N) + Z(1:N,istage)
 551             CALL FUN_CHEM(T+rkC(istage)*H,TMP,RHS)
 552             TMP(1:N) = G(1:N) - Z(1:N,istage)
 553             CALL WAXPY(N,HGammaInv,TMP,1,RHS,1) ! RHS(:) <- RHS(:) + HGammaInv*(G(:)-Z(:))
 554
 555 !~~~>     SOLVE THE LINEAR SYSTEMS
 556 #ifdef FULL_ALGEBRA
 557             CALL DGETRS( 'N', N, 1, E, N, IP, RHS, N, IER )
 558 #else
 559             CALL KppSolve(E, RHS)
 560 #endif
 561             Nsol=Nsol+1
 562
 563 !~~~>     CHECK CONVERGENCE OR IF NUMBER OF ITERATIONS TOO LARGE
 564             CALL SDIRK_ErrorNorm(N, RHS, SCAL, NewtonErr)
 565             IF ( (NewtonIter >= 2) .AND. (NewtonIter < NewtonMaxit) ) THEN
 566                 Theta = NewtonErr/NewtonErrOld
 567                 IF (Theta < 0.99d0) THEN
 568                     NewtonFactor(istage)=Theta/(ONE-Theta)
 569                     DYTH = NewtonFactor(istage)*NewtonErr* &
 570                            Theta**(NewtonMaxit-1-NewtonIter)
 571                     QNEWT = MAX(1.0d-4,MIN(16.0d0,DYTH/NewtonTol))
 572                     IF (QNEWT >= ONE) THEN
 573                          H=.8D0*H*QNEWT**(-ONE/(NewtonMaxit-NewtonIter))
 574                          Reject=.TRUE.
 575                          NewtonReject=.TRUE.
 576                          CYCLE Tloop ! go back to the beginning of DO step
 577                     END IF
 578                 ELSE
 579                     NewtonReject=.TRUE.
 580                     H=H*0.5d0
 581                     Reject=.TRUE.
 582                     CYCLE Tloop ! go back to the beginning of DO step
 583                 END IF
 584             END IF
 585             NewtonErrOld = NewtonErr
 586             CALL WAXPY(N,ONE,RHS,1,Z(1,istage),1) ! Z(:) <-- Z(:)+RHS(:)
 587
 588             END DO Newton
 589
 590 !~~> END OF SIMPLIFIED NEWTON ITERATION
 591             ! Save function values
 592             TMP(1:N) = Y(1:N) + Z(1:N,istage)
 593             CALL FUN_CHEM(T+rkC(istage)*H,TMP,FV(1,istage))
 594
 595    END DO stages
 596
 597
 598 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 599 !  ERROR ESTIMATION
 600 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 601       Nstp=Nstp+1
 602       TMP(1:N)=HGammaInv*(Z(1:N,5)-Yhat(1:N))
 603
 604 #ifdef FULL_ALGEBRA
 605       CALL DGETRS( 'N', N, 1, E, N, IP, TMP, N, IER )
 606 #else
 607       CALL KppSolve(E, TMP)
 608 #endif
 609
 610       CALL SDIRK_ErrorNorm(N, TMP, SCAL, ERR)
 611
 612 !~~~> COMPUTATION OF Hnew: WE REQUIRE FacMin <= Hnew/H <= FacMax
 613       !Safe = FacSafe*DBLE(1+2*NewtonMaxit)/DBLE(NewtonIter+2*NewtonMaxit)
 614       Fac  = MAX(FacMin,MIN(FacMax,(ERR)**(-0.25d0)*FacSafe))
 615       Hnew = H*Fac
 616
 617 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 618 !  ACCEPT/Reject STEP
 619 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 620 accept:  IF ( ERR < ONE ) THEN !~~~> STEP IS ACCEPTED
 621
 622          FIRST=.FALSE.
 623          Nacc=Nacc+1
 624          Hold=H
 625
 626 !~~~> COEFFICIENTS FOR CONTINUOUS SOLUTION
 627          CALL WAXPY(N,ONE,Z(1,5),1,Y,1) ! Y(:) <-- Y(:)+Z5(:)
 628          CALL WCOPY(N,Z(1,5),1,Yhat,1)  ! Yhat <-- Z5
 629
 630          DO i=1,4  ! CONTi <-- Sum_j rkD(i,j)*Zj
 631            CALL Set2zero(N,CONT(1,i))
 632            DO j = 1,5
 633              CALL WAXPY(N,rkD(i,j),Z(1,j),1,CONT(1,i),1)
 634            END DO
 635          END DO
 636
 637          CALL SDIRK_ErrorScale(ITOL, AbsTol, RelTol, Y, SCAL)
 638
 639          T=T+H
 640          FreshJac=.FALSE.
 641
 642          Hnew = Tdirection*MIN(ABS(Hnew),Hmax1)
 643          Hexit = Hnew
 644          IF (Reject) Hnew=Tdirection*MIN(ABS(Hnew),ABS(H))
 645          Reject = .FALSE.
 646          NewtonReject = .FALSE.
 647          IF ((T+Hnew/Qmin-Tfinal)*Tdirection > 0.D0) THEN
 648             H = Tfinal-T
 649          ELSE
 650             Hratio=Hnew/H
 651             ! If step not changed too much, keep it as is;
 652             !    do not update Jacobian and reuse LU
 653             IF ( (Theta <= ThetaMin) .AND. (Hratio >= Qmin) &
 654                                      .AND. (Hratio <= Qmax) ) THEN
 655                SkipJacUpdate = .TRUE.
 656                SkipLU        = .TRUE.
 657             ELSE
 658                H = Hnew
 659             END IF
 660          END IF
 661          ! If convergence is fast enough, do not update Jacobian
 662          IF (Theta <= ThetaMin) SkipJacUpdate = .TRUE.
 663
 664       ELSE accept !~~~> STEP IS REJECTED
 665
 666          Reject=.TRUE.
 667          IF (FIRST) THEN
 668              H=H*FacRej
 669          ELSE
 670              H=Hnew
 671          END IF
 672          IF (Nacc >= 1) Nrej=Nrej+1
 673
 674       END IF accept
 675
 676       END DO Tloop
 677
 678       ! Successful return
 679       Texit = T
 680       IDID  = 1
 681
 682       END SUBROUTINE SDIRK_Integrator
 683
 684
 685 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 686       SUBROUTINE SDIRK_ErrorScale(ITOL, AbsTol, RelTol, Y, SCAL)
 687 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 688       IMPLICIT NONE
 689       INTEGER :: i, ITOL
 690       KPP_REAL :: AbsTol(NVAR), RelTol(NVAR), &
 691                        Y(NVAR), SCAL(NVAR)
 692       IF (ITOL == 0) THEN
 693         DO i=1,NVAR
 694           SCAL(i) = 1.0d0 / ( AbsTol(1)+RelTol(1)*ABS(Y(i)) )
 695         END DO
 696       ELSE
 697         DO i=1,NVAR
 698           SCAL(i) = 1.0d0 / ( AbsTol(i)+RelTol(i)*ABS(Y(i)) )
 699         END DO
 700       END IF
 701       END SUBROUTINE SDIRK_ErrorScale
 702
 703
 704 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 705       SUBROUTINE SDIRK_ErrorNorm(N, Y, SCAL, ERR)
 706 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 707 !
 708       INTEGER :: i, N
 709       KPP_REAL :: Y(N), SCAL(N), ERR
 710       ERR=0.0d0
 711       DO i=1,N
 712            ERR = ERR+(Y(i)*SCAL(i))**2
 713       END DO
 714       ERR = MAX( SQRT(ERR/DBLE(N)), 1.0d-10 )
 715 !
 716       END SUBROUTINE SDIRK_ErrorNorm
 717
 718
 719 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 720  SUBROUTINE SDIRK_ErrorMsg(Code,T,H,IERR)
 721 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 722 !    Handles all error messages
 723 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 724
 725    KPP_REAL, INTENT(IN) :: T, H
 726    INTEGER, INTENT(IN)  :: Code
 727    INTEGER, INTENT(OUT) :: IERR
 728
 729    IERR = Code
 730    PRINT * , &
 731      'Forced exit from SDIRK due to the following error:'
 732    IF ((Code>=-8).AND.(Code<=-1)) THEN
 733      PRINT *, IERR_NAMES(Code)
 734    ELSE
 735      PRINT *, 'Unknown Error code: ', Code
 736    ENDIF
 737
 738    PRINT *, "T=", T, "and H=", H
 739
 740  END SUBROUTINE SDIRK_ErrorMsg
 741
 742
 743 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 744       SUBROUTINE SDIRK_PrepareMatrix ( H, T, Y, FJAC, &
 745                    FreshJac, SkipJacUpdate, E, IP, Reject, IER )
 746 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 747 !
 748       IMPLICIT NONE
 749
 750       KPP_REAL, INTENT(INOUT) :: H
 751       KPP_REAL, INTENT(IN)    :: T, Y(NVAR)
 752       LOGICAL, INTENT(INOUT)       :: FreshJac, SkipJacUpdate
 753       INTEGER, INTENT(OUT)         :: IER, IP(NVAR)
 754       LOGICAL, INTENT(INOUT)       :: Reject
 755 #ifdef FULL_ALGEBRA
 756       KPP_REAL, INTENT(INOUT) :: FJAC(NVAR,NVAR)
 757       KPP_REAL, INTENT(OUT)   :: E(NVAR,NVAR)
 758 #else
 759       KPP_REAL, INTENT(INOUT) :: FJAC(LU_NONZERO)
 760       KPP_REAL, INTENT(OUT)   :: E(LU_NONZERO)
 761 #endif
 762       KPP_REAL                :: HGammaInv
 763       INTEGER                      :: i, j, ConsecutiveSng
 764       KPP_REAL, PARAMETER     :: ONE = 1.0d0
 765
 766  20   CONTINUE
 767
 768 !~~~>  COMPUTE THE JACOBIAN
 769       IF (SkipJacUpdate) THEN
 770           SkipJacUpdate = .FALSE.
 771       ELSE IF ( .NOT.FreshJac  ) THEN
 772           CALL JAC_CHEM( T, Y, FJAC )
 773           FreshJac = .TRUE.
 774       END IF
 775
 776 !~~~>  Compute the matrix E = 1/(H*GAMMA)*Jac, and its decomposition
 777       ConsecutiveSng = 0
 778       IER = 1
 779
 780 Hloop: DO WHILE (IER /= 0)
 781
 782       HGammaInv = ONE/(H*rkGamma)
 783
 784 #ifdef FULL_ALGEBRA
 785       DO j=1,NVAR
 786          DO i=1,NVAR
 787             E(i,j)=-FJAC(i,j)
 788          END DO
 789          E(j,j)=E(j,j)+HGammaInv
 790       END DO
 791       CALL DGETRF( NVAR, NVAR, E, NVAR, IP, IER )
 792 #else
 793       DO i = 1,LU_NONZERO
 794             E(i) = -FJAC(i)
 795       END DO
 796       DO i = 1,NVAR
 797           j = LU_DIAG(i); E(j) = E(j) + HGammaInv
 798       END DO
 799       CALL KppDecomp ( E, IER)
 800       IP(1) = 1
 801 #endif
 802       Ndec=Ndec+1
 803
 804       IF (IER /= 0) THEN
 805           WRITE (6,*) ' MATRIX IS SINGULAR, IER=',IER,' T=',T,' H=',H
 806           Nsng = Nsng+1; ConsecutiveSng = ConsecutiveSng + 1
 807           IF (ConsecutiveSng >= 6) RETURN ! Failure
 808           H=H*0.5d0
 809           Reject=.TRUE.
 810           !~~~> Update Jacobian if not fresh
 811           IF ( .NOT.FreshJac ) GOTO 20
 812       END IF
 813
 814       END DO Hloop
 815
 816       END SUBROUTINE SDIRK_PrepareMatrix
 817
 818
 819 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 820       SUBROUTINE SDIRK_Coefficients
 821 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 822           rkGamma=4.0d0/15.0d0
 823
 824           rkA(1,1)= 4.d0/15.d0
 825           rkA(2,1)= 1.d0/2.d0
 826           rkA(2,2)= 4.d0/15.d0
 827           rkA(3,1)= 51069.d0/144200.d0
 828           rkA(3,2)=-7809.d0/144200.d0
 829           rkA(3,3)= 4.d0/15.d0
 830           rkA(4,1)= 12047244770625658.d0/141474406359725325.d0
 831           rkA(4,2)=-3057890203562191.d0/47158135453241775.d0
 832           rkA(4,3)= 2239631894905804.d0/28294881271945065.d0
 833           rkA(4,4)= 4.d0/15.d0
 834           rkA(5,1)= 181513.d0/86430.d0
 835           rkA(5,2)=-89074.d0/116015.d0
 836           rkA(5,3)= 83636.d0/34851.d0
 837           rkA(5,4)=-69863904375173.d0/23297141763930.d0
 838           rkA(5,5)= 4.d0/15.d0
 839
 840           rkB(1)= 181513.d0/86430.d0
 841           rkB(2)=-89074.d0/116015.d0
 842           rkB(3)= 83636.d0/34851.d0
 843           rkB(4)=-69863904375173.d0/23297141763930.d0
 844           rkB(5)= 4/15.d0
 845
 846           rkC(1)=4.d0/15.d0
 847           rkC(2)=23.d0/30.d0
 848           rkC(3)=17.d0/30.d0
 849           rkC(4)=707.d0/1931.d0
 850           rkC(5)=1.d0
 851
 852           rkBeta(2,1)=15.0d0/8.0d0
 853           rkBeta(3,1)=1577061.0d0/922880.0d0
 854           rkBeta(3,2)=-23427.0d0/115360.0d0
 855           rkBeta(4,1)=647163682356923881.0d0/2414496535205978880.0d0
 856           rkBeta(4,2)=-593512117011179.0d0/3245291041943520.0d0
 857           rkBeta(4,3)=559907973726451.0d0/1886325418129671.0d0
 858           rkBeta(5,1)=724545451.0d0/796538880.0d0
 859           rkBeta(5,2)=-830832077.0d0/267298560.0d0
 860           rkBeta(5,3)=30957577.0d0/2509272.0d0
 861           rkBeta(5,4)=-69863904375173.0d0/6212571137048.0d0
 862
 863           rkAlpha(2,1)= 23.d0/8.d0
 864           rkAlpha(3,1)= 0.9838473040915402d0
 865           rkAlpha(3,2)= 0.3969226768377252d0
 866           rkAlpha(4,1)= 0.6563374010466914d0
 867           rkAlpha(4,2)= 0.0d0
 868           rkAlpha(4,3)= 0.3372498196189311d0
 869           rkAlpha(5,1)=7752107607.0d0/11393456128.0d0
 870           rkAlpha(5,2)=-17881415427.0d0/11470078208.0d0
 871           rkAlpha(5,3)=2433277665.0d0/179459416.0d0
 872           rkAlpha(5,4)=-96203066666797.0d0/6212571137048.0d0
 873
 874           rkD(1,1)= 24.74416644927758d0
 875           rkD(1,2)= -4.325375951824688d0
 876           rkD(1,3)= 41.39683763286316d0
 877           rkD(1,4)= -61.04144619901784d0
 878           rkD(1,5)= -3.391332232917013d0
 879           rkD(2,1)= -51.98245719616925d0
 880           rkD(2,2)= 10.52501981094525d0
 881           rkD(2,3)= -154.2067922191855d0
 882           rkD(2,4)= 214.3082125319825d0
 883           rkD(2,5)= 14.71166018088679d0
 884           rkD(3,1)= 33.14347947522142d0
 885           rkD(3,2)= -19.72986789558523d0
 886           rkD(3,3)= 230.4878502285804d0
 887           rkD(3,4)= -287.6629744338197d0
 888           rkD(3,5)= -18.99932366302254d0
 889           rkD(4,1)= -5.905188728329743d0
 890           rkD(4,2)= 13.53022403646467d0
 891           rkD(4,3)= -117.6778956422581d0
 892           rkD(4,4)= 134.3962081008550d0
 893           rkD(4,5)= 8.678995715052762d0
 894
 895       END SUBROUTINE SDIRK_Coefficients
 896
 897 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 898    END SUBROUTINE SDIRK ! AND ALL ITS INTERNAL PROCEDURES
 899 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 900
 901 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 902       SUBROUTINE FUN_CHEM( T, Y, P )
 903 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 904
 905       USE KPP_ROOT_Global, ONLY: NVAR
 906       USE KPP_ROOT_Function
 907
 908       INTEGER N
 909       KPP_REAL   T !, Told
 910       KPP_REAL   Y(NVAR), P(NVAR)
 911
 912       !Told = TIME
 913       !TIME = T
 914       !CALL Update_SUN()
 915       !CALL Update_RCONST()
 916
 917       CALL Fun( Y,  FIX, RCONST, P )
 918
 919       !TIME = Told
 920       Nfun=Nfun+1
 921
 922       END SUBROUTINE FUN_CHEM
 923
 924
 925 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 926       SUBROUTINE JAC_CHEM( T, Y, JV )
 927 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 928
 929       USE KPP_ROOT_Global, ONLY: NVAR
 930       USE KPP_ROOT_Jacobian
 931
 932       INTEGER N
 933       KPP_REAL   T !, Told
 934       KPP_REAL   Y(NVAR)
 935 #ifdef FULL_ALGEBRA
 936       KPP_REAL :: JS(LU_NONZERO), JV(NVAR,NVAR)
 937       INTEGER :: i, j
 938 #else
 939       KPP_REAL :: JV(LU_NONZERO)
 940 #endif
 941
 942       !Told = TIME
 943       !TIME = T
 944       !CALL Update_SUN()
 945       !CALL Update_RCONST()
 946
 947 #ifdef FULL_ALGEBRA
 948       CALL Jac_SP(Y, FIX, RCONST, JS)
 949       DO j=1,NVAR
 950          DO j=1,NVAR
 951           JV(i,j) = 0.0D0
 952          END DO
 953       END DO
 954       DO i=1,LU_NONZERO
 955          JV(LU_IROW(i),LU_ICOL(i)) = JS(i)
 956       END DO
 957 #else
 958       CALL Jac_SP(Y, FIX, RCONST, JV)
 959 #endif
 960       !TIME = Told
 961       Njac = Njac+1
 962
 963       END SUBROUTINE JAC_CHEM
 964
 965 END MODULE KPP_ROOT_Integrator