source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/kpp_sdirk4.f90 @ 4654

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

Merge of branch palm4u into trunk

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