source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/sdirk_tlm.f90 @ 4072

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

Merge of branch palm4u into trunk

File size: 50.6 KB
Line 
1!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
2!  SDIRK-TLM - Tangent Linear Model of Singly-Diagonally-Implicit RK      !
3!            * Sdirk 2a, 2b: L-stable, 2 stages, order 2                  !
4!            * Sdirk 3a:     L-stable, 3 stages, order 2, adj-invariant   !
5!            * Sdirk 4a, 4b: L-stable, 5 stages, order 4                  !
6!  By default the code employs the KPP sparse linear algebra routines     !
7!  Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK)        !
8!                                                                         !
9!    (C)  Adrian Sandu, July 2005                                         !
10!    Virginia Polytechnic Institute and State University                  !
11!    Contact: sandu@cs.vt.edu                                             !
12!    Revised by Philipp Miehe and Adrian Sandu, May 2006                  !
13!    This implementation is part of KPP - the Kinetic PreProcessor        !
14!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
15
16MODULE KPP_ROOT_Integrator
17
18  USE KPP_ROOT_Precision
19  USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME
20  USE KPP_ROOT_Parameters, ONLY: NVAR, NSPEC, NFIX, LU_NONZERO
21  USE KPP_ROOT_JacobianSP, ONLY: LU_DIAG
22  USE KPP_ROOT_Jacobian, ONLY:   Jac_SP_Vec
23  USE KPP_ROOT_LinearAlgebra, ONLY: KppDecomp,    &
24               KppSolve, Set2zero, WLAMCH, WCOPY, WAXPY, WSCAL, WADD
25 
26  IMPLICIT NONE
27  PUBLIC
28  SAVE
29 
30!~~~>  Statistics on the work performed by the SDIRK method
31  INTEGER, PARAMETER :: Nfun=1, Njac=2, Nstp=3, Nacc=4,  &
32           Nrej=5, Ndec=6, Nsol=7, Nsng=8,               &
33           Ntexit=1, Nhexit=2, Nhnew=3
34                 
35CONTAINS
36
37!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
38SUBROUTINE INTEGRATE_TLM( NTLM, Y, Y_tlm, TIN, TOUT, ATOL_tlm,RTOL_tlm, &
39       ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U, IERR_U )
40!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
41
42    USE KPP_ROOT_Parameters, ONLY: NVAR,ind_O3
43    USE KPP_ROOT_Global,     ONLY: ATOL,RTOL,VAR
44
45    IMPLICIT NONE
46
47!~~~> Y - Concentrations
48    KPP_REAL :: Y(NVAR)
49!~~~> NTLM - No. of sensitivity coefficients
50    INTEGER NTLM
51!~~~> Y_tlm - Sensitivities of concentrations
52!     Note: Y_tlm (1:NVAR,j) contains sensitivities of
53!               Y(1:NVAR) w.r.t. the j-th parameter, j=1...NTLM
54    KPP_REAL :: Y_tlm(NVAR,NTLM)
55    KPP_REAL :: TIN  ! TIN - Start Time
56    KPP_REAL :: TOUT ! TOUT - End Time
57    KPP_REAL, INTENT(IN),  OPTIONAL :: RTOL_tlm(NVAR,NTLM),ATOL_tlm(NVAR,NTLM)
58    INTEGER,       INTENT(IN),  OPTIONAL :: ICNTRL_U(20)
59    KPP_REAL, INTENT(IN),  OPTIONAL :: RCNTRL_U(20)
60    INTEGER,       INTENT(OUT), OPTIONAL :: ISTATUS_U(20)
61    KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20)
62    INTEGER,       INTENT(OUT), OPTIONAL :: IERR_U
63
64    INTEGER :: IERR
65
66    KPP_REAL :: RCNTRL(20), RSTATUS(20), T1, T2
67    INTEGER       :: ICNTRL(20), ISTATUS(20)
68    INTEGER, SAVE :: Ntotal = 0
69
70    ICNTRL(1:20) = 0
71    RCNTRL(1:20) = 0.0_dp
72
73    !~~~> fine-tune the integrator:
74    ICNTRL(2)  = 0   ! 0=vector tolerances, 1=scalar tolerances
75    ICNTRL(5)  = 8   ! Max no. of Newton iterations
76    ICNTRL(6)  = 0   ! Starting values for Newton are interpolated (0) or zero (1)
77    ICNTRL(7)  = 0   ! How to solve TLM: 0=modified Newton, 1=direct
78    ICNTRL(9)  = 0   ! TLM Newton Iterations influence
79    ICNTRL(12) = 0   ! TLM Truncation Error influence
80
81    !~~~> if optional parameters are given, and if they are >0,
82    !     then use them to overwrite default settings
83    IF (PRESENT(ICNTRL_U)) THEN
84      WHERE(ICNTRL_U(:) > 0) ICNTRL(:) = ICNTRL_U(:)
85    END IF
86    IF (PRESENT(RCNTRL_U)) THEN
87      WHERE(RCNTRL_U(:) > 0) RCNTRL(:) = RCNTRL_U(:)
88    END IF
89
90
91    T1 = TIN; T2 = TOUT
92    CALL SdirkTLM(  NVAR, NTLM, T1, T2, Y, Y_tlm, RTOL, ATOL, &
93                    RTOL_tlm, ATOL_tlm, RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR  )
94
95    !~~~> Debug option: print number of steps
96    ! Ntotal = Ntotal + ISTATUS(Nstp)
97    ! PRINT*,'NSTEPS=',ISTATUS(Nstp),' (',Ntotal,')','  O3=', VAR(ind_O3)
98
99    ! if optional parameters are given for output
100    ! use them to store information in them
101    IF (PRESENT(ISTATUS_U)) ISTATUS_U(:) = ISTATUS(:)
102    IF (PRESENT(RSTATUS_U)) RSTATUS_U(:) = RSTATUS(:)
103    IF (PRESENT(IERR_U)) IERR_U = IERR
104
105    IF (IERR < 0) THEN
106      PRINT *,'SDIRK-TLM: Unsuccessful exit at T=', TIN,' (IERR=',IERR,')'
107    ENDIF
108
109  END SUBROUTINE INTEGRATE_TLM
110
111
112
113!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
114   SUBROUTINE SdirkTLM(N, NTLM, Tinitial, Tfinal, Y, Y_tlm, RelTol, AbsTol, &
115                          RelTol_tlm, AbsTol_tlm, &
116                          RCNTRL, ICNTRL, RSTATUS, ISTATUS, Ierr)
117!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
118!
119!    Solves the system y'=F(t,y) using a Singly-Diagonally-Implicit
120!    Runge-Kutta (SDIRK) method.
121!
122!    This implementation is based on the book and the code Sdirk4:
123!
124!      E. Hairer and G. Wanner
125!      "Solving ODEs II. Stiff and differential-algebraic problems".
126!      Springer series in computational mathematics, Springer-Verlag, 1996.
127!    This code is based on the SDIRK4 routine in the above book.
128!
129!    Methods:
130!            * Sdirk 2a, 2b: L-stable, 2 stages, order 2                 
131!            * Sdirk 3a:     L-stable, 3 stages, order 2, adjoint-invariant   
132!            * Sdirk 4a, 4b: L-stable, 5 stages, order 4                 
133!
134!    (C)  Adrian Sandu, July 2005
135!    Virginia Polytechnic Institute and State University
136!    Contact: sandu@cs.vt.edu
137!    Revised by Philipp Miehe and Adrian Sandu, May 2006                 
138!    This implementation is part of KPP - the Kinetic PreProcessor
139!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
140!
141!~~~>   INPUT ARGUMENTS:
142!
143!-     Y(NVAR)    = vector of initial conditions (at T=Tinitial)
144!-    [Tinitial,Tfinal]  = time range of integration
145!     (if Tinitial>Tfinal the integration is performed backwards in time)
146!-    RelTol, AbsTol = user precribed accuracy
147!- SUBROUTINE ode_Fun( T, Y, Ydot ) = ODE function,
148!                       returns Ydot = Y' = F(T,Y)
149!- SUBROUTINE ode_Fun( T, Y, Ydot ) = Jacobian of the ODE function,
150!                       returns Jcb = dF/dY
151!-    ICNTRL(1:20)    = integer inputs parameters
152!-    RCNTRL(1:20)    = real inputs parameters
153!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
154!
155!~~~>     OUTPUT ARGUMENTS:
156!
157!-    Y(NVAR)         -> vector of final states (at T->Tfinal)
158!-    ISTATUS(1:20)   -> integer output parameters
159!-    RSTATUS(1:20)   -> real output parameters
160!-    Ierr            -> job status upon return
161!                        success (positive value) or
162!                        failure (negative value)
163!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
164!
165!~~~>     INPUT PARAMETERS:
166!
167!    Note: For input parameters equal to zero the default values of the
168!       corresponding variables are used.
169!
170!    Note: For input parameters equal to zero the default values of the
171!          corresponding variables are used.
172!~~~> 
173!    ICNTRL(1) = not used
174!
175!    ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors
176!              = 1: AbsTol, RelTol are scalars
177!
178!    ICNTRL(3) = Method
179!
180!    ICNTRL(4)  -> maximum number of integration steps
181!        For ICNTRL(4)=0 the default value of 100000 is used
182!
183!    ICNTRL(5)  -> maximum number of Newton iterations
184!        For ICNTRL(5)=0 the default value of 8 is used
185!
186!    ICNTRL(6)  -> starting values of Newton iterations:
187!        ICNTRL(6)=0 : starting values are interpolated (the default)
188!        ICNTRL(6)=1 : starting values are zero
189!
190!    ICNTRL(7)  -> method to solve TLM equations:
191!        ICNTRL(7)=0 : modified Newton re-using LU (the default)
192!        ICNTRL(7)=1 : direct solution (additional one LU factorization per stage)
193!
194!    ICNTRL(9) -> switch for TLM Newton iteration error estimation strategy
195!               ICNTRL(9) = 0: base number of iterations as forward solution
196!               ICNTRL(9) = 1: use RTOL_tlm and ATOL_tlm to calculate
197!                               error estimation for TLM at Newton stages
198!
199!    ICNTRL(12) -> switch for TLM truncation error control
200!               ICNTRL(12) = 0: TLM error is not used
201!               ICNTRL(12) = 1: TLM error is computed and used
202!
203!
204!~~~>  Real parameters
205!
206!    RCNTRL(1)  -> Hmin, lower bound for the integration step size
207!                  It is strongly recommended to keep Hmin = ZERO
208!    RCNTRL(2)  -> Hmax, upper bound for the integration step size
209!    RCNTRL(3)  -> Hstart, starting value for the integration step size
210!
211!    RCNTRL(4)  -> FacMin, lower bound on step decrease factor (default=0.2)
212!    RCNTRL(5)  -> FacMax, upper bound on step increase factor (default=6)
213!    RCNTRL(6)  -> FacRej, step decrease factor after multiple rejections
214!                 (default=0.1)
215!    RCNTRL(7)  -> FacSafe, by which the new step is slightly smaller
216!                  than the predicted value  (default=0.9)
217!    RCNTRL(8)  -> ThetaMin. If Newton convergence rate smaller
218!                  than ThetaMin the Jacobian is not recomputed;
219!                  (default=0.001)
220!    RCNTRL(9)  -> NewtonTol, stopping criterion for Newton's method
221!                  (default=0.03)
222!    RCNTRL(10) -> Qmin
223!    RCNTRL(11) -> Qmax. If Qmin < Hnew/Hold < Qmax, then the
224!                  step size is kept constant and the LU factorization
225!                  reused (default Qmin=1, Qmax=1.2)
226!
227!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
228!
229!~~~>     OUTPUT PARAMETERS:
230!
231!    Note: each call to Rosenbrock adds the current no. of fcn calls
232!      to previous value of ISTATUS(1), and similar for the other params.
233!      Set ISTATUS(1:10) = 0 before call to avoid this accumulation.
234!
235!    ISTATUS(1) = No. of function calls
236!    ISTATUS(2) = No. of jacobian calls
237!    ISTATUS(3) = No. of steps
238!    ISTATUS(4) = No. of accepted steps
239!    ISTATUS(5) = No. of rejected steps (except at the beginning)
240!    ISTATUS(6) = No. of LU decompositions
241!    ISTATUS(7) = No. of forward/backward substitutions
242!    ISTATUS(8) = No. of singular matrix decompositions
243!
244!    RSTATUS(1)  -> Texit, the time corresponding to the
245!                   computed Y upon return
246!    RSTATUS(2)  -> Hexit,last accepted step before return
247!    RSTATUS(3)  -> Hnew, last predicted step before return
248!        For multiple restarts, use Hnew as Hstart in the following run
249!
250!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
251      IMPLICIT NONE
252
253! Arguments     
254      INTEGER, INTENT(IN)          :: N, NTLM, ICNTRL(20)
255      KPP_REAL, INTENT(IN)    :: Tinitial, Tfinal, &
256                    RelTol(N), AbsTol(N), RCNTRL(20), &
257                    RelTol_tlm(N,NTLM), AbsTol_tlm(N,NTLM)
258      KPP_REAL, INTENT(INOUT) :: Y(NVAR), Y_tlm(N,NTLM)
259      INTEGER, INTENT(OUT)         :: Ierr
260      INTEGER, INTENT(INOUT)       :: ISTATUS(20) 
261      KPP_REAL, INTENT(OUT)   :: RSTATUS(20)
262       
263!~~~>  SDIRK method coefficients, up to 5 stages
264      INTEGER, PARAMETER :: Smax = 5
265      INTEGER, PARAMETER :: S2A=1, S2B=2, S3A=3, S4A=4, S4B=5
266      KPP_REAL :: rkGamma, rkA(Smax,Smax), rkB(Smax), rkC(Smax), &
267                       rkD(Smax),  rkE(Smax), rkBhat(Smax), rkELO,    &
268                       rkAlpha(Smax,Smax), rkTheta(Smax,Smax)
269      INTEGER :: sdMethod, rkS ! The number of stages
270
271! Local variables     
272      KPP_REAL :: Hmin, Hmax, Hstart, Roundoff,    &
273                       FacMin, Facmax, FacSafe, FacRej, &
274                       ThetaMin, NewtonTol, Qmin, Qmax
275      INTEGER       :: ITOL, NewtonMaxit, Max_no_steps, i
276      LOGICAL       :: StartNewton, DirectTLM, TLMNewtonEst, TLMtruncErr
277      KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
278
279       
280      Ierr = 0
281      ISTATUS(1:20) = 0
282      RSTATUS(1:20) = ZERO
283
284!~~~>  For Scalar tolerances (ICNTRL(2).NE.0)  the code uses AbsTol(1) and RelTol(1)
285!   For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:NVAR) and RelTol(1:NVAR)
286      IF (ICNTRL(2) == 0) THEN
287         ITOL = 1
288      ELSE
289         ITOL = 0
290      END IF
291
292!~~~> ICNTRL(3) - method selection       
293      SELECT CASE (ICNTRL(3))
294      CASE (0,1)
295         CALL Sdirk2a
296      CASE (2)
297         CALL Sdirk2b
298      CASE (3)
299         CALL Sdirk3a
300      CASE (4)
301         CALL Sdirk4a
302      CASE (5)
303         CALL Sdirk4b
304      CASE DEFAULT
305         CALL Sdirk2a
306      END SELECT
307     
308!~~~>   The maximum number of time steps admitted
309      IF (ICNTRL(4) == 0) THEN
310         Max_no_steps = 200000
311      ELSEIF (ICNTRL(4) > 0) THEN
312         Max_no_steps=ICNTRL(4)
313      ELSE
314         PRINT * ,'User-selected ICNTRL(4)=',ICNTRL(4)
315         CALL SDIRK_ErrorMsg(-1,Tinitial,ZERO,Ierr)
316   END IF
317   !~~~> The maximum number of Newton iterations admitted
318      IF(ICNTRL(5) == 0)THEN
319         NewtonMaxit=8
320      ELSE
321         NewtonMaxit=ICNTRL(5)
322         IF(NewtonMaxit <= 0)THEN
323             PRINT * ,'User-selected ICNTRL(5)=',ICNTRL(5)
324             CALL SDIRK_ErrorMsg(-2,Tinitial,ZERO,Ierr)
325         END IF
326      END IF
327!~~~> StartNewton:  Use extrapolation for starting values of Newton iterations
328      IF (ICNTRL(6) == 0) THEN
329         StartNewton = .TRUE.
330      ELSE
331         StartNewton = .FALSE.
332      END IF     
333!~~~> Solve TLM equations directly or by Newton iterations
334      DirectTLM = (ICNTRL(7) == 1)
335!~~~> Newton iteration error control selection 
336      IF (ICNTRL(9) == 0) THEN
337         TLMNewtonEst = .FALSE.
338      ELSE
339         TLMNewtonEst = .TRUE.
340      END IF     
341!~~~> TLM truncation error control selection 
342      IF (ICNTRL(12) == 0) THEN
343         TLMtruncErr = .FALSE.
344      ELSE
345         TLMtruncErr = .TRUE.
346      END IF     
347           
348!~~~>  Unit roundoff (1+Roundoff>1)
349      Roundoff = WLAMCH('E')
350
351!~~~>  Lower bound on the step size: (positive value)
352      IF (RCNTRL(1) == ZERO) THEN
353         Hmin = ZERO
354      ELSEIF (RCNTRL(1) > ZERO) THEN
355         Hmin = RCNTRL(1)
356      ELSE
357         PRINT * , 'User-selected RCNTRL(1)=', RCNTRL(1)
358         CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr)
359      END IF
360   
361!~~~>  Upper bound on the step size: (positive value)
362      IF (RCNTRL(2) == ZERO) THEN
363         Hmax = ABS(Tfinal-Tinitial)
364      ELSEIF (RCNTRL(2) > ZERO) THEN
365         Hmax = MIN(ABS(RCNTRL(2)),ABS(Tfinal-Tinitial))
366      ELSE
367         PRINT * , 'User-selected RCNTRL(2)=', RCNTRL(2)
368         CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr)
369      END IF
370   
371!~~~>  Starting step size: (positive value)
372      IF (RCNTRL(3) == ZERO) THEN
373         Hstart = MAX(Hmin,Roundoff)
374      ELSEIF (RCNTRL(3) > ZERO) THEN
375         Hstart = MIN(ABS(RCNTRL(3)),ABS(Tfinal-Tinitial))
376      ELSE
377         PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3)
378         CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr)
379      END IF
380   
381!~~~>  Step size can be changed s.t.  FacMin < Hnew/Hexit < FacMax
382      IF (RCNTRL(4) == ZERO) THEN
383         FacMin = 0.2_dp
384      ELSEIF (RCNTRL(4) > ZERO) THEN
385         FacMin = RCNTRL(4)
386      ELSE
387         PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4)
388         CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr)
389      END IF
390      IF (RCNTRL(5) == ZERO) THEN
391         FacMax = 10.0_dp
392      ELSEIF (RCNTRL(5) > ZERO) THEN
393         FacMax = RCNTRL(5)
394      ELSE
395         PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5)
396         CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr)
397      END IF
398!~~~>   FacRej: Factor to decrease step after 2 succesive rejections
399      IF (RCNTRL(6) == ZERO) THEN
400         FacRej = 0.1_dp
401      ELSEIF (RCNTRL(6) > ZERO) THEN
402         FacRej = RCNTRL(6)
403      ELSE
404         PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6)
405         CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr)
406      END IF
407!~~~>   FacSafe: Safety Factor in the computation of new step size
408      IF (RCNTRL(7) == ZERO) THEN
409         FacSafe = 0.9_dp
410      ELSEIF (RCNTRL(7) > ZERO) THEN
411         FacSafe = RCNTRL(7)
412      ELSE
413         PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7)
414         CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr)
415      END IF
416
417!~~~> ThetaMin: decides whether the Jacobian should be recomputed
418      IF(RCNTRL(8) == 0.D0)THEN
419         ThetaMin = 1.0d-3
420      ELSE
421         ThetaMin = RCNTRL(8)
422      END IF
423
424!~~~> Stopping criterion for Newton's method
425      IF(RCNTRL(9) == ZERO)THEN
426         NewtonTol = 3.0d-2
427      ELSE
428         NewtonTol = RCNTRL(9)
429      END IF
430
431!~~~> Qmin, Qmax: IF Qmin < Hnew/Hold < Qmax, STEP SIZE = CONST.
432      IF(RCNTRL(10) == ZERO)THEN
433         Qmin=ONE
434      ELSE
435         Qmin=RCNTRL(10)
436      END IF
437      IF(RCNTRL(11) == ZERO)THEN
438         Qmax=1.2D0
439      ELSE
440         Qmax=RCNTRL(11)
441      END IF
442
443!~~~>  Check if tolerances are reasonable
444      IF (ITOL == 0) THEN
445         IF (AbsTol(1) <= ZERO .OR. RelTol(1) <= 10.D0*Roundoff) THEN
446            PRINT * , ' Scalar AbsTol = ',AbsTol(1)
447            PRINT * , ' Scalar RelTol = ',RelTol(1)
448            CALL SDIRK_ErrorMsg(-5,Tinitial,ZERO,Ierr)
449         END IF
450      ELSE
451         DO i=1,N
452            IF (AbsTol(i) <= 0.D0.OR.RelTol(i) <= 10.D0*Roundoff) THEN
453              PRINT * , ' AbsTol(',i,') = ',AbsTol(i)
454              PRINT * , ' RelTol(',i,') = ',RelTol(i)
455              CALL SDIRK_ErrorMsg(-5,Tinitial,ZERO,Ierr)
456            END IF
457         END DO
458      END IF
459   
460    IF (Ierr < 0) RETURN
461
462    CALL SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr )
463
464
465!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
466   CONTAINS  !  PROCEDURES INTERNAL TO SDIRK
467!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
468
469
470!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
471   SUBROUTINE SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr )
472!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
473
474      USE KPP_ROOT_Parameters
475      IMPLICIT NONE
476
477!~~~> Arguments:     
478      INTEGER, INTENT(IN) :: N, NTLM
479      KPP_REAL, INTENT(INOUT) :: Y(N), Y_tlm(N,NTLM)
480      KPP_REAL, INTENT(IN) :: Tinitial, Tfinal
481      INTEGER, INTENT(OUT) :: Ierr
482     
483!~~~> Local variables:     
484      KPP_REAL :: Z(NVAR,rkS), G(NVAR), TMP(NVAR),         &   
485                       NewtonRate, SCAL(NVAR), DZ(NVAR),        &
486                       T, H, Theta, Hratio, NewtonPredictedErr, &
487                       Qnewton, Err, Fac, Hnew, Tdirection,     &
488                       NewtonIncrement, NewtonIncrementOld, &
489                       SCAL_tlm(NVAR), Yerr(N), Yerr_tlm(N,NTLM), ThetaTLM
490      KPP_REAL :: Z_tlm(NVAR,rkS,NTLM)
491      INTEGER :: itlm, j, IER, istage, NewtonIter, saveNiter, NewtonIterTLM
492      INTEGER :: IP(NVAR), IP_tlm(NVAR)
493      LOGICAL :: Reject, FirstStep, SkipJac, SkipLU, NewtonDone
494     
495#ifdef FULL_ALGEBRA     
496      KPP_REAL, DIMENSION(NVAR,NVAR)  :: FJAC, E, Jac, E_tlm
497#else     
498      KPP_REAL, DIMENSION(LU_NONZERO) :: FJAC, E, Jac, E_tlm
499#endif     
500      KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
501
502
503!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
504!~~~>  Initializations
505!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
506
507      T = Tinitial
508      Tdirection = SIGN(ONE,Tfinal-Tinitial)
509      H = MAX(ABS(Hmin),ABS(Hstart))
510      IF (ABS(H) <= 10.D0*Roundoff) H=1.0D-6
511      H=MIN(ABS(H),Hmax)
512      H=SIGN(H,Tdirection)
513      SkipLU  = .FALSE.
514      SkipJac = .FALSE.
515      Reject =  .FALSE.
516      FirstStep=.TRUE.
517
518      CALL SDIRK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL)
519
520!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
521!~~~>  Time loop begins
522!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
523Tloop: DO WHILE ( (Tfinal-T)*Tdirection - Roundoff > ZERO )
524
525
526!~~~>  Compute E = 1/(h*gamma)-Jac and its LU decomposition
527      IF ( .NOT.SkipLU ) THEN ! This time around skip the Jac update and LU
528         CALL SDIRK_PrepareMatrix ( H, T, Y, FJAC, &
529                   SkipJac, SkipLU, E, IP, Reject, IER )
530         IF (IER /= 0) THEN
531             CALL SDIRK_ErrorMsg(-8,T,H,Ierr); RETURN
532         END IF
533      END IF     
534
535      IF (ISTATUS(Nstp) > Max_no_steps) THEN
536             CALL SDIRK_ErrorMsg(-6,T,H,Ierr); RETURN
537      END IF   
538      IF ( (T+0.1d0*H == T) .OR. (ABS(H) <= Roundoff) ) THEN
539             CALL SDIRK_ErrorMsg(-7,T,H,Ierr); RETURN
540      END IF   
541
542stages:DO istage = 1, rkS
543
544!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
545!~~~>  Simplified Newton iterations
546!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
547
548!~~~>  Starting values for Newton iterations
549       CALL Set2zero(N,Z(1,istage))
550       
551!~~~>   Prepare the loop-independent part of the right-hand side
552       CALL Set2zero(N,G)
553       IF (istage > 1) THEN
554           DO j = 1, istage-1
555               ! Gj(:) = sum_j Theta(i,j)*Zj(:) = H * sum_j A(i,j)*Fun(Zj(:))
556               CALL WAXPY(N,rkTheta(istage,j),Z(1,j),1,G,1)
557               ! Zi(:) = sum_j Alpha(i,j)*Zj(:)
558               IF (StartNewton) THEN
559                 CALL WAXPY(N,rkAlpha(istage,j),Z(1,j),1,Z(1,istage),1)
560               END IF
561           END DO
562       END IF
563
564       !~~~>  Initializations for Newton iteration
565       NewtonDone = .FALSE.
566       Fac = 0.5d0 ! Step reduction factor if too many iterations
567           
568NewtonLoop:DO NewtonIter = 1, NewtonMaxit
569
570!~~~>   Prepare the loop-dependent part of the right-hand side
571            CALL WADD(N,Y,Z(1,istage),TMP)              ! TMP <- Y + Zi
572            CALL FUN_CHEM(T+rkC(istage)*H,TMP,DZ)       ! DZ <- Fun(Y+Zi)
573            ISTATUS(Nfun) = ISTATUS(Nfun) + 1
574!            DZ(1:N) = G(1:N) - Z(1:N,istage) + (H*rkGamma)*DZ(1:N)
575            CALL WSCAL(N, H*rkGamma, DZ, 1)
576            CALL WAXPY (N, -ONE, Z(1,istage), 1, DZ, 1)
577            CALL WAXPY (N, ONE, G,1, DZ,1)
578
579!~~~>   Solve the linear system
580            CALL SDIRK_Solve ( H, N, E, IP, IER, DZ )
581           
582!~~~>   Check convergence of Newton iterations
583            CALL SDIRK_ErrorNorm(N, DZ, SCAL, NewtonIncrement)
584            IF ( NewtonIter == 1 ) THEN
585                Theta      = ABS(ThetaMin)
586                NewtonRate = 2.0d0 
587            ELSE
588                Theta = NewtonIncrement/NewtonIncrementOld
589                IF (Theta < 0.99d0) THEN
590                    NewtonRate = Theta/(ONE-Theta)
591                    ! Predict error at the end of Newton process
592                    NewtonPredictedErr = NewtonIncrement &
593                               *Theta**(NewtonMaxit-NewtonIter)/(ONE-Theta)
594                    IF (NewtonPredictedErr >= NewtonTol) THEN
595                      ! Non-convergence of Newton: predicted error too large
596                      Qnewton = MIN(10.0d0,NewtonPredictedErr/NewtonTol)
597                      Fac = 0.8d0*Qnewton**(-ONE/(1+NewtonMaxit-NewtonIter))
598                      EXIT NewtonLoop
599                    END IF
600                ELSE ! Non-convergence of Newton: Theta too large
601                    EXIT NewtonLoop
602                END IF
603            END IF
604            NewtonIncrementOld = NewtonIncrement
605            ! Update solution: Z(:) <-- Z(:)+DZ(:)
606            CALL WAXPY(N,ONE,DZ,1,Z(1,istage),1) 
607           
608            ! Check error in Newton iterations
609            NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol)
610            IF (NewtonDone) THEN
611               ! Tune error in TLM variables by defining the minimal number of Newton iterations.
612               saveNiter = NewtonIter+1
613               EXIT NewtonLoop
614            END IF
615           
616            END DO NewtonLoop
617           
618            IF (.NOT.NewtonDone) THEN
619                 !CALL RK_ErrorMsg(-12,T,H,Ierr);
620                 H = Fac*H; Reject=.TRUE.
621                 SkipJac = .TRUE.; SkipLU = .FALSE.
622                 CYCLE Tloop
623            END IF
624
625!~~~>  End of simplified Newton iterations for forward variables
626
627!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
628!~~~>  Solve for TLM variables
629!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
630
631!~~~>  Direct solution for TLM variables
632DirTLM:IF (DirectTLM) THEN
633
634         TMP(1:N) = Y(1:N) + Z(1:N,istage)
635         SkipJac = .FALSE.
636         CALL SDIRK_PrepareMatrix ( H, T+rkC(istage)*H, TMP, Jac, &
637                   SkipJac, SkipLU, E_tlm, IP_tlm, Reject, IER )
638         IF (IER /= 0) CYCLE TLoop
639
640TlmL:    DO itlm = 1, NTLM
641            G(1:N) = Y_tlm(1:N,itlm)
642            IF (istage > 1) THEN
643              ! Gj(:) = sum_j Theta(i,j)*Zj_tlm(:)
644              !       = H * sum_j A(i,j)*Jac(Zj(:))*(Yj_tlm+Zj_tlm)
645              DO j = 1, istage-1
646                  CALL WAXPY(N,rkTheta(istage,j),Z_tlm(1,j,itlm),1,G,1)
647              END DO
648            END IF
649            CALL SDIRK_Solve ( H, N, E_tlm, IP_tlm, IER, G )
650            Z_tlm(1:N,istage,itlm) = G(1:N) - Y_tlm(1:N,itlm)
651         END DO TlmL
652
653       ELSE DirTLM
654
655!~~~>  Jacobian of the current stage solution
656       TMP(1:N) = Y(1:N) + Z(1:N,istage)
657       CALL JAC_CHEM(T+rkC(istage)*H,TMP,Jac)
658       ISTATUS(Njac) = ISTATUS(Njac) + 1
659
660!~~~>  Simplified Newton iterations for TLM variables
661TlmLoop:DO itlm = 1,NTLM
662       NewtonRate = MAX(NewtonRate,Roundoff)**0.8d0
663
664!~~~>  Starting values for Newton iterations
665       CALL Set2zero(N,Z_tlm(1,istage,itlm))
666       
667!~~~>   Prepare the loop-independent part of the right-hand side
668#ifdef FULL_ALGEBRA 
669       DZ = MATMUL(Jac,Y_tlm(1,itlm))    ! DZ <- Jac(Y+Z)*Y_tlm
670#else     
671       CALL Jac_SP_Vec ( Jac, Y_tlm(1,itlm), DZ )   
672#endif     
673       G(1:N) = (H*rkGamma)*DZ(1:N)
674       IF (istage > 1) THEN
675           ! Gj(:) = sum_j Theta(i,j)*Zj_tlm(:)
676           !       = H * sum_j A(i,j)*Jac(Zj(:))*(Yj_tlm+Zj_tlm)
677           DO j = 1, istage-1
678               CALL WAXPY(N,rkTheta(istage,j),Z_tlm(1,j,itlm),1,G,1)
679           END DO
680       END IF
681       
682       
683       !~~~>  Initializations for Newton iteration
684       IF (TLMNewtonEst) THEN
685          NewtonDone = .FALSE.
686          Fac = 0.5d0 ! Step reduction factor if too many iterations
687
688          CALL SDIRK_ErrorScale(N,ITOL,AbsTol_tlm(1,itlm),RelTol_tlm(1,itlm), &
689               Y_tlm(1,itlm),SCAL_tlm)
690       END IF
691           
692NewtonLoopTLM:DO NewtonIterTLM = 1, NewtonMaxit
693
694!~~~>   Prepare the loop-dependent part of the right-hand side
695#ifdef FULL_ALGEBRA 
696            DZ = MATMUL(Jac,Z_tlm(1,istage,itlm))    ! DZ <- Jac(Y+Z)*Z_tlm
697#else     
698            CALL Jac_SP_Vec ( Jac, Z_tlm(1,istage,itlm), DZ )   
699#endif     
700            DZ(1:N) = (H*rkGamma)*DZ(1:N)+G(1:N)-Z_tlm(1:N,istage,itlm)
701
702            CALL SDIRK_Solve ( H, N, E, IP, IER, DZ )
703           
704            IF (TLMNewtonEst) THEN
705!~~~>   Check convergence of Newton iterations
706            CALL SDIRK_ErrorNorm(N, DZ, SCAL_tlm, NewtonIncrement)
707            IF ( NewtonIterTLM <= 1 ) THEN
708                ThetaTLM      = ABS(ThetaMin)
709                NewtonRate = 2.0d0 
710            ELSE
711                ThetaTLM = NewtonIncrement/NewtonIncrementOld
712                IF (ThetaTLM < 0.99d0) THEN
713                    NewtonRate = ThetaTLM/(ONE-ThetaTLM)
714                    ! Predict error at the end of Newton process
715                    NewtonPredictedErr = NewtonIncrement &
716                               *ThetaTLM**(NewtonMaxit-NewtonIterTLM)/(ONE-ThetaTLM)
717                    IF (NewtonPredictedErr >= NewtonTol) THEN
718                      ! Non-convergence of Newton: predicted error too large
719                      Qnewton = MIN(10.0d0,NewtonPredictedErr/NewtonTol)
720                      Fac = 0.8d0*Qnewton**(-ONE/(1+NewtonMaxit-NewtonIterTLM))
721                      EXIT NewtonLoopTLM
722                    END IF
723                ELSE ! Non-convergence of Newton: Theta too large
724                    EXIT NewtonLoopTLM
725                END IF
726            END IF
727            NewtonIncrementOld = NewtonIncrement
728            END IF !(TLMNewtonEst)
729           
730            ! Update solution: Z_tlm(:) <-- Z_tlm(:)+DZ(:)
731            CALL WAXPY(N,ONE,DZ,1,Z_tlm(1,istage,itlm),1) 
732           
733            ! Check error in Newton iterations
734            IF (TLMNewtonEst) THEN
735               NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol)
736               IF (NewtonDone) EXIT NewtonLoopTLM
737            ELSE
738               ! Minimum number of iterations same as FWD iterations
739               IF (NewtonIterTLM>=saveNiter) EXIT NewtonLoopTLM
740            END IF
741           
742            END DO NewtonLoopTLM
743           
744            IF ((TLMNewtonEst) .AND. (.NOT.NewtonDone)) THEN
745                 !CALL RK_ErrorMsg(-12,T,H,Ierr);
746                 H = Fac*H; Reject=.TRUE.
747                 SkipJac = .TRUE.; SkipLU = .FALSE.
748                 CYCLE Tloop
749            END IF
750
751      END DO TlmLoop
752!~~~> End of simplified Newton iterations for TLM
753
754    END IF DirTLM
755
756   END DO stages
757
758!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
759!~~~>  Error estimation
760!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
761      ISTATUS(Nstp) = ISTATUS(Nstp) + 1
762      CALL Set2zero(N,Yerr)
763      DO i = 1,rkS
764         IF (rkE(i)/=ZERO) CALL WAXPY(N,rkE(i),Z(1,i),1,Yerr,1)
765      END DO 
766
767      CALL SDIRK_Solve ( H, N, E, IP, IER, Yerr )
768      CALL SDIRK_ErrorNorm(N, Yerr, SCAL, Err)
769
770      IF (TLMtruncErr) THEN
771        CALL Set2zero(NVAR*NTLM,Yerr_tlm)
772        DO itlm=1,NTLM
773          DO j=1,rkS 
774            IF (rkE(j) /= ZERO) CALL WAXPY(N,rkE(j),Z_tlm(1,j,itlm),1,Yerr_tlm(1,itlm),1)
775          END DO
776          CALL SDIRK_Solve (H, N, E, IP, IER, Yerr_tlm(1,itlm))
777        END DO
778        CALL SDIRK_ErrorNorm_tlm(N,NTLM, Yerr_tlm, Err)
779      END IF
780
781!~~~> Computation of new step size Hnew
782      Fac  = FacSafe*(Err)**(-ONE/rkELO)
783      Fac  = MAX(FacMin,MIN(FacMax,Fac))
784      Hnew = H*Fac
785
786!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
787!~~~>  Accept/Reject step
788!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
789accept: IF ( Err < ONE ) THEN !~~~> Step is accepted
790
791         FirstStep=.FALSE.
792         ISTATUS(Nacc) = ISTATUS(Nacc) + 1
793
794!~~~> Update time and solution
795         T  =  T + H
796         ! Y(:) <-- Y(:) + Sum_j rkD(j)*Z_j(:)
797         DO i = 1,rkS 
798            IF (rkD(i)/=ZERO) THEN
799                 CALL WAXPY(N,rkD(i),Z(1,i),1,Y,1)
800                 DO itlm = 1, NTLM
801                   CALL WAXPY(N,rkD(i),Z_tlm(1,i,itlm),1,Y_tlm(1,itlm),1)
802                 END DO 
803            END IF     
804         END DO 
805       
806!~~~> Update scaling coefficients
807         CALL SDIRK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL)
808
809!~~~> Next time step
810         Hnew = Tdirection*MIN(ABS(Hnew),Hmax)
811         ! Last T and H
812         RSTATUS(Ntexit) = T
813         RSTATUS(Nhexit) = H
814         RSTATUS(Nhnew)  = Hnew
815         ! No step increase after a rejection
816         IF (Reject) Hnew = Tdirection*MIN(ABS(Hnew),ABS(H))
817         Reject = .FALSE.
818         IF ((T+Hnew/Qmin-Tfinal)*Tdirection > ZERO) THEN
819            H = Tfinal-T
820         ELSE
821            Hratio=Hnew/H
822            ! If step not changed too much keep Jacobian and reuse LU
823            SkipLU = ( (Theta <= ThetaMin) .AND. (Hratio >= Qmin) &
824                                     .AND. (Hratio <= Qmax) )
825            ! For TLM: do not skip LU (decrease TLM error)
826            SkipLU = .FALSE.
827            IF (.NOT.SkipLU) H = Hnew
828         END IF
829         ! If convergence is fast enough, do not update Jacobian
830!         SkipJac = (Theta <= ThetaMin)
831         SkipJac = .FALSE.
832
833      ELSE accept !~~~> Step is rejected
834
835         IF (FirstStep .OR. Reject) THEN
836             H = FacRej*H
837         ELSE
838             H = Hnew
839         END IF
840         Reject = .TRUE.
841         SkipJac = .TRUE.
842         SkipLU  = .FALSE. 
843         IF (ISTATUS(Nacc) >= 1) ISTATUS(Nrej) = ISTATUS(Nrej) + 1
844         
845      END IF accept
846     
847      END DO Tloop
848
849      ! Successful return
850      Ierr  = 1
851 
852      END SUBROUTINE SDIRK_IntegratorTLM
853
854
855!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
856      SUBROUTINE SDIRK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL)
857!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
858      IMPLICIT NONE
859      INTEGER :: N, i, ITOL
860      KPP_REAL :: AbsTol(NVAR), RelTol(NVAR), &
861                       Y(NVAR), SCAL(NVAR)
862      IF (ITOL == 0) THEN
863        DO i=1,N
864          SCAL(i) = ONE / ( AbsTol(1)+RelTol(1)*ABS(Y(i)) )
865        END DO
866      ELSE
867        DO i=1,N
868          SCAL(i) = ONE / ( AbsTol(i)+RelTol(i)*ABS(Y(i)) )
869        END DO
870      END IF
871      END SUBROUTINE SDIRK_ErrorScale
872     
873     
874!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
875      SUBROUTINE SDIRK_ErrorNorm(N, Y, SCAL, Err)
876!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
877!     
878      INTEGER :: i, N
879      KPP_REAL :: Y(N), SCAL(N), Err     
880      Err = ZERO
881      DO i=1,N
882           Err = Err+(Y(i)*SCAL(i))**2
883      END DO
884      Err = MAX( SQRT(Err/DBLE(N)), 1.0d-10 )
885!
886      END SUBROUTINE SDIRK_ErrorNorm
887
888!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
889      SUBROUTINE SDIRK_ErrorNorm_tlm(N,NTLM, Y_tlm, FWD_Err)
890!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
891!     
892      INTEGER :: itlm, NTLM, N
893      KPP_REAL :: Y_tlm(N,NTLM), SCAL_tlm(N), FWD_Err, Err 
894     
895      DO itlm=1,NTLM
896        CALL SDIRK_ErrorScale(N,ITOL,AbsTol_tlm(1,itlm),RelTol_tlm(1,itlm), &
897               Y_tlm(1,itlm),SCAL_tlm)
898        CALL SDIRK_ErrorNorm(N, Y_tlm(1,itlm), SCAL_tlm, Err)
899        FWD_Err = MAX(FWD_Err, Err)
900      END DO
901!
902      END SUBROUTINE SDIRK_ErrorNorm_tlm
903
904
905!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
906 SUBROUTINE SDIRK_ErrorMsg(Code,T,H,Ierr)
907!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
908!    Handles all error messages
909!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
910
911   KPP_REAL, INTENT(IN) :: T, H
912   INTEGER, INTENT(IN)  :: Code
913   INTEGER, INTENT(OUT) :: Ierr
914
915   Ierr = Code
916   PRINT * , &
917     'Forced exit from SDIRK due to the following error:'
918
919   SELECT CASE (Code)
920    CASE (-1)
921      PRINT * , '--> Improper value for maximal no of steps'
922    CASE (-2)
923      PRINT * , '--> Improper value for maximal no of Newton iterations'
924    CASE (-3)
925      PRINT * , '--> Hmin/Hmax/Hstart must be positive'
926    CASE (-4)
927      PRINT * , '--> FacMin/FacMax/FacRej must be positive'
928    CASE (-5)
929      PRINT * , '--> Improper tolerance values'
930    CASE (-6)
931      PRINT * , '--> No of steps exceeds maximum bound'
932    CASE (-7)
933      PRINT * , '--> Step size too small: T + 10*H = T', &
934            ' or H < Roundoff'
935    CASE (-8)
936      PRINT * , '--> Matrix is repeatedly singular'
937    CASE DEFAULT
938      PRINT *, 'Unknown Error code: ', Code
939   END SELECT
940
941   PRINT *, "T=", T, "and H=", H
942
943 END SUBROUTINE SDIRK_ErrorMsg
944     
945     
946!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
947      SUBROUTINE SDIRK_PrepareMatrix ( H, T, Y, FJAC, &
948                   SkipJac, SkipLU, E, IP, Reject, ISING )
949!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
950!~~~>  Compute the matrix E = I - 1/(H*Gamma)*Jac, and its decomposition
951!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
952     
953      IMPLICIT NONE
954     
955      KPP_REAL, INTENT(INOUT) :: H
956      KPP_REAL, INTENT(IN)    :: T, Y(NVAR)
957      LOGICAL, INTENT(INOUT)       :: SkipJac,SkipLU,Reject
958      INTEGER, INTENT(OUT)         :: ISING, IP(NVAR)
959#ifdef FULL_ALGEBRA
960      KPP_REAL, INTENT(INOUT) :: FJAC(NVAR,NVAR)
961      KPP_REAL, INTENT(OUT)   :: E(NVAR,NVAR)
962#else
963      KPP_REAL, INTENT(INOUT) :: FJAC(LU_NONZERO)
964      KPP_REAL, INTENT(OUT)   :: E(LU_NONZERO)
965#endif
966      KPP_REAL                :: HGammaInv
967      INTEGER                      :: i, j, ConsecutiveSng
968
969      ConsecutiveSng = 0
970      ISING = 1
971     
972Hloop: DO WHILE (ISING /= 0)
973     
974      HGammaInv = ONE/(H*rkGamma)
975
976!~~~>  Compute the Jacobian
977!      IF (SkipJac) THEN
978!          SkipJac = .FALSE.
979!      ELSE
980      IF (.NOT. SkipJac) THEN
981          CALL JAC_CHEM( T, Y, FJAC )
982          ISTATUS(Njac) = ISTATUS(Njac) + 1
983      END IF 
984     
985#ifdef FULL_ALGEBRA
986      DO j=1,NVAR
987         DO i=1,NVAR
988            E(i,j) = -FJAC(i,j)
989         END DO
990         E(j,j) = E(j,j)+HGammaInv
991      END DO
992      CALL DGETRF( NVAR, NVAR, E, NVAR, IP, ISING )
993#else
994      DO i = 1,LU_NONZERO
995            E(i) = -FJAC(i)
996      END DO
997      DO i = 1,NVAR
998          j = LU_DIAG(i); E(j) = E(j) + HGammaInv
999      END DO
1000      CALL KppDecomp ( E, ISING)
1001      IP(1) = 1
1002#endif
1003      ISTATUS(Ndec) = ISTATUS(Ndec) + 1
1004
1005      IF (ISING /= 0) THEN
1006          WRITE (6,*) ' MATRIX IS SINGULAR, ISING=',ISING,' T=',T,' H=',H
1007          ISTATUS(Nsng) = ISTATUS(Nsng) + 1; ConsecutiveSng = ConsecutiveSng + 1
1008          IF (ConsecutiveSng >= 6) RETURN ! Failure
1009          H = 0.5d0*H
1010          SkipJac = .TRUE.
1011          SkipLU  = .FALSE.
1012          Reject  = .TRUE.
1013      END IF
1014     
1015      END DO Hloop
1016
1017      END SUBROUTINE SDIRK_PrepareMatrix
1018
1019!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1020      SUBROUTINE SDIRK_Solve ( H, N, E, IP, ISING, RHS )
1021!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1022!~~~>  Solves the system (H*Gamma-Jac)*x = R
1023!      using the LU decomposition of E = I - 1/(H*Gamma)*Jac
1024!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1025      IMPLICIT NONE
1026      INTEGER, INTENT(IN)       :: N, IP(N), ISING
1027      KPP_REAL, INTENT(IN) :: H
1028#ifdef FULL_ALGEBRA
1029      KPP_REAL, INTENT(IN) :: E(NVAR,NVAR)
1030#else
1031      KPP_REAL, INTENT(IN) :: E(LU_NONZERO)
1032#endif
1033      KPP_REAL, INTENT(INOUT) :: RHS(N)
1034      KPP_REAL                :: HGammaInv
1035     
1036      HGammaInv = ONE/(H*rkGamma)
1037      CALL WSCAL(N,HGammaInv,RHS,1)
1038#ifdef FULL_ALGEBRA 
1039      CALL DGETRS( 'N', N, 1, E, N, IP, RHS, N, ISING )
1040#else
1041      CALL KppSolve(E, RHS)
1042#endif
1043      ISTATUS(Nsol) = ISTATUS(Nsol) + 1
1044 
1045      END SUBROUTINE SDIRK_Solve
1046
1047!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1048      SUBROUTINE Sdirk4a
1049!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1050      sdMethod = S4A
1051! Number of stages
1052      rkS = 5
1053
1054! Method coefficients
1055      rkGamma = .2666666666666666666666666666666667d0
1056
1057      rkA(1,1) = .2666666666666666666666666666666667d0
1058      rkA(2,1) = .5000000000000000000000000000000000d0
1059      rkA(2,2) = .2666666666666666666666666666666667d0
1060      rkA(3,1) = .3541539528432732316227461858529820d0
1061      rkA(3,2) = -.5415395284327323162274618585298197d-1
1062      rkA(3,3) = .2666666666666666666666666666666667d0
1063      rkA(4,1) = .8515494131138652076337791881433756d-1
1064      rkA(4,2) = -.6484332287891555171683963466229754d-1
1065      rkA(4,3) = .7915325296404206392428857585141242d-1
1066      rkA(4,4) = .2666666666666666666666666666666667d0
1067      rkA(5,1) = 2.100115700566932777970612055999074d0
1068      rkA(5,2) = -.7677800284445976813343102185062276d0
1069      rkA(5,3) = 2.399816361080026398094746205273880d0
1070      rkA(5,4) = -2.998818699869028161397714709433394d0
1071      rkA(5,5) = .2666666666666666666666666666666667d0
1072
1073      rkB(1)   = 2.100115700566932777970612055999074d0
1074      rkB(2)   = -.7677800284445976813343102185062276d0
1075      rkB(3)   = 2.399816361080026398094746205273880d0
1076      rkB(4)   = -2.998818699869028161397714709433394d0
1077      rkB(5)   = .2666666666666666666666666666666667d0
1078
1079      rkBhat(1)= 2.885264204387193942183851612883390d0
1080      rkBhat(2)= -.1458793482962771337341223443218041d0
1081      rkBhat(3)= 2.390008682465139866479830743628554d0
1082      rkBhat(4)= -4.129393538556056674929560012190140d0
1083      rkBhat(5)= 0.d0
1084
1085      rkC(1)   = .2666666666666666666666666666666667d0
1086      rkC(2)   = .7666666666666666666666666666666667d0
1087      rkC(3)   = .5666666666666666666666666666666667d0
1088      rkC(4)   = .3661315380631796996374935266701191d0
1089      rkC(5)   = 1.d0
1090
1091! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i}
1092      rkD(1)   = 0.d0
1093      rkD(2)   = 0.d0
1094      rkD(3)   = 0.d0
1095      rkD(4)   = 0.d0
1096      rkD(5)   = 1.d0
1097
1098! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i}
1099      rkE(1)   = -.6804000050475287124787034884002302d0
1100      rkE(2)   = 1.558961944525217193393931795738823d0
1101      rkE(3)   = -13.55893003128907927748632408763868d0
1102      rkE(4)   = 15.48522576958521253098585004571302d0
1103      rkE(5)   = 1.d0
1104
1105! Local order of Err estimate
1106      rkElo    = 4
1107
1108! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j}
1109      rkTheta(2,1) = 1.875000000000000000000000000000000d0
1110      rkTheta(3,1) = 1.708847304091539528432732316227462d0
1111      rkTheta(3,2) = -.2030773231622746185852981969486824d0
1112      rkTheta(4,1) = .2680325578937783958847157206823118d0
1113      rkTheta(4,2) = -.1828840955527181631794050728644549d0
1114      rkTheta(4,3) = .2968246986151577397160821594427966d0
1115      rkTheta(5,1) = .9096171815241460655379433581446771d0
1116      rkTheta(5,2) = -3.108254967778352416114774430509465d0
1117      rkTheta(5,3) = 12.33727431701306195581826123274001d0
1118      rkTheta(5,4) = -11.24557012450885560524143016037523d0
1119
1120! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j}
1121      rkAlpha(2,1) = 2.875000000000000000000000000000000d0
1122      rkAlpha(3,1) = .8500000000000000000000000000000000d0
1123      rkAlpha(3,2) = .4434782608695652173913043478260870d0
1124      rkAlpha(4,1) = .7352046091658870564637910527807370d0
1125      rkAlpha(4,2) = -.9525565003057343527941920657462074d-1
1126      rkAlpha(4,3) = .4290111305453813852259481840631738d0
1127      rkAlpha(5,1) = -16.10898993405067684831655675112808d0
1128      rkAlpha(5,2) = 6.559571569643355712998131800797873d0
1129      rkAlpha(5,3) = -15.90772144271326504260996815012482d0
1130      rkAlpha(5,4) = 25.34908987169226073668861694892683d0
1131               
1132!~~~> Coefficients for continuous solution
1133!          rkD(1,1)= 24.74416644927758d0
1134!          rkD(1,2)= -4.325375951824688d0
1135!          rkD(1,3)= 41.39683763286316d0
1136!          rkD(1,4)= -61.04144619901784d0
1137!          rkD(1,5)= -3.391332232917013d0
1138!          rkD(2,1)= -51.98245719616925d0
1139!          rkD(2,2)= 10.52501981094525d0
1140!          rkD(2,3)= -154.2067922191855d0
1141!          rkD(2,4)= 214.3082125319825d0
1142!          rkD(2,5)= 14.71166018088679d0
1143!          rkD(3,1)= 33.14347947522142d0
1144!          rkD(3,2)= -19.72986789558523d0
1145!          rkD(3,3)= 230.4878502285804d0
1146!          rkD(3,4)= -287.6629744338197d0
1147!          rkD(3,5)= -18.99932366302254d0
1148!          rkD(4,1)= -5.905188728329743d0
1149!          rkD(4,2)= 13.53022403646467d0
1150!          rkD(4,3)= -117.6778956422581d0
1151!          rkD(4,4)= 134.3962081008550d0
1152!          rkD(4,5)= 8.678995715052762d0
1153!
1154!         DO i=1,4  ! CONTi <-- Sum_j rkD(i,j)*Zj
1155!           CALL Set2zero(N,CONT(1,i))
1156!           DO j = 1,rkS
1157!             CALL WAXPY(N,rkD(i,j),Z(1,j),1,CONT(1,i),1)
1158!           END DO
1159!         END DO
1160         
1161          rkELO = 4.0d0
1162         
1163      END SUBROUTINE Sdirk4a
1164
1165!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1166      SUBROUTINE Sdirk4b
1167!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1168      sdMethod = S4B
1169! Number of stages
1170      rkS = 5
1171
1172! Method coefficients
1173      rkGamma = .25d0
1174
1175      rkA(1,1) = 0.25d0
1176      rkA(2,1) = 0.5d00
1177      rkA(2,2) = 0.25d0
1178      rkA(3,1) = 0.34d0
1179      rkA(3,2) =-0.40d-1
1180      rkA(3,3) = 0.25d0
1181      rkA(4,1) = 0.2727941176470588235294117647058824d0
1182      rkA(4,2) =-0.5036764705882352941176470588235294d-1
1183      rkA(4,3) = 0.2757352941176470588235294117647059d-1
1184      rkA(4,4) = 0.25d0
1185      rkA(5,1) = 1.041666666666666666666666666666667d0
1186      rkA(5,2) =-1.020833333333333333333333333333333d0
1187      rkA(5,3) = 7.812500000000000000000000000000000d0
1188      rkA(5,4) =-7.083333333333333333333333333333333d0
1189      rkA(5,5) = 0.25d0
1190
1191      rkB(1)   =  1.041666666666666666666666666666667d0
1192      rkB(2)   = -1.020833333333333333333333333333333d0
1193      rkB(3)   =  7.812500000000000000000000000000000d0
1194      rkB(4)   = -7.083333333333333333333333333333333d0
1195      rkB(5)   =  0.250000000000000000000000000000000d0
1196
1197      rkBhat(1)=  1.069791666666666666666666666666667d0
1198      rkBhat(2)= -0.894270833333333333333333333333333d0
1199      rkBhat(3)=  7.695312500000000000000000000000000d0
1200      rkBhat(4)= -7.083333333333333333333333333333333d0
1201      rkBhat(5)=  0.212500000000000000000000000000000d0
1202
1203      rkC(1)   = 0.25d0
1204      rkC(2)   = 0.75d0
1205      rkC(3)   = 0.55d0
1206      rkC(4)   = 0.50d0
1207      rkC(5)   = 1.00d0
1208
1209! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i}
1210      rkD(1)   = 0.0d0
1211      rkD(2)   = 0.0d0
1212      rkD(3)   = 0.0d0
1213      rkD(4)   = 0.0d0
1214      rkD(5)   = 1.0d0
1215
1216! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i}
1217      rkE(1)   =  0.5750d0
1218      rkE(2)   =  0.2125d0
1219      rkE(3)   = -4.6875d0
1220      rkE(4)   =  4.2500d0
1221      rkE(5)   =  0.1500d0
1222
1223! Local order of Err estimate
1224      rkElo    = 4
1225
1226! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j}
1227      rkTheta(2,1) = 2.d0
1228      rkTheta(3,1) = 1.680000000000000000000000000000000d0
1229      rkTheta(3,2) = -.1600000000000000000000000000000000d0
1230      rkTheta(4,1) = 1.308823529411764705882352941176471d0
1231      rkTheta(4,2) = -.1838235294117647058823529411764706d0
1232      rkTheta(4,3) = 0.1102941176470588235294117647058824d0
1233      rkTheta(5,1) = -3.083333333333333333333333333333333d0
1234      rkTheta(5,2) = -4.291666666666666666666666666666667d0
1235      rkTheta(5,3) =  34.37500000000000000000000000000000d0
1236      rkTheta(5,4) = -28.33333333333333333333333333333333d0
1237
1238! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j}
1239      rkAlpha(2,1) = 3.
1240      rkAlpha(3,1) = .8800000000000000000000000000000000d0
1241      rkAlpha(3,2) = .4400000000000000000000000000000000d0
1242      rkAlpha(4,1) = .1666666666666666666666666666666667d0
1243      rkAlpha(4,2) = -.8333333333333333333333333333333333d-1
1244      rkAlpha(4,3) = .9469696969696969696969696969696970d0
1245      rkAlpha(5,1) = -6.d0
1246      rkAlpha(5,2) = 9.d0
1247      rkAlpha(5,3) = -56.81818181818181818181818181818182d0
1248      rkAlpha(5,4) = 54.d0
1249     
1250      END SUBROUTINE Sdirk4b
1251
1252!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1253      SUBROUTINE Sdirk2a
1254!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1255      sdMethod = S2A
1256! Number of stages
1257      rkS = 2
1258
1259! Method coefficients
1260      rkGamma = .2928932188134524755991556378951510d0
1261
1262      rkA(1,1) = .2928932188134524755991556378951510d0
1263      rkA(2,1) = .7071067811865475244008443621048490d0
1264      rkA(2,2) = .2928932188134524755991556378951510d0
1265
1266      rkB(1)   = .7071067811865475244008443621048490d0
1267      rkB(2)   = .2928932188134524755991556378951510d0
1268
1269      rkBhat(1)= .6666666666666666666666666666666667d0
1270      rkBhat(2)= .3333333333333333333333333333333333d0
1271
1272      rkC(1)   = 0.292893218813452475599155637895151d0
1273      rkC(2)   = 1.0d0
1274
1275! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i}
1276      rkD(1)   = 0.0d0
1277      rkD(2)   = 1.0d0
1278
1279! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i}
1280      rkE(1)   =  0.4714045207910316829338962414032326d0
1281      rkE(2)   = -0.1380711874576983496005629080698993d0
1282
1283! Local order of Err estimate
1284      rkElo    = 2
1285
1286! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j}
1287      rkTheta(2,1) = 2.414213562373095048801688724209698d0
1288
1289! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j}
1290      rkAlpha(2,1) = 3.414213562373095048801688724209698d0
1291         
1292      END SUBROUTINE Sdirk2a
1293
1294!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1295      SUBROUTINE Sdirk2b
1296!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1297      sdMethod = S2B
1298! Number of stages
1299      rkS      = 2
1300
1301! Method coefficients
1302      rkGamma  = 1.707106781186547524400844362104849d0
1303
1304      rkA(1,1) = 1.707106781186547524400844362104849d0
1305      rkA(2,1) = -.707106781186547524400844362104849d0
1306      rkA(2,2) = 1.707106781186547524400844362104849d0
1307
1308      rkB(1)   = -.707106781186547524400844362104849d0
1309      rkB(2)   = 1.707106781186547524400844362104849d0
1310
1311      rkBhat(1)= .6666666666666666666666666666666667d0
1312      rkBhat(2)= .3333333333333333333333333333333333d0
1313
1314      rkC(1)   = 1.707106781186547524400844362104849d0
1315      rkC(2)   = 1.0d0
1316
1317! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i}
1318      rkD(1)   = 0.0d0
1319      rkD(2)   = 1.0d0
1320
1321! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i}
1322      rkE(1)   = -.4714045207910316829338962414032326d0
1323      rkE(2)   =  .8047378541243650162672295747365659d0
1324
1325! Local order of Err estimate
1326      rkElo    = 2
1327
1328! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j}
1329      rkTheta(2,1) = -.414213562373095048801688724209698d0
1330
1331! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j}
1332      rkAlpha(2,1) = .5857864376269049511983112757903019d0
1333     
1334      END SUBROUTINE Sdirk2b
1335
1336
1337!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1338      SUBROUTINE Sdirk3a
1339!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1340      sdMethod = S3A
1341! Number of stages
1342      rkS = 3
1343
1344! Method coefficients
1345      rkGamma = .2113248654051871177454256097490213d0
1346
1347      rkA(1,1) = .2113248654051871177454256097490213d0
1348      rkA(2,1) = .2113248654051871177454256097490213d0
1349      rkA(2,2) = .2113248654051871177454256097490213d0
1350      rkA(3,1) = .2113248654051871177454256097490213d0
1351      rkA(3,2) = .5773502691896257645091487805019573d0
1352      rkA(3,3) = .2113248654051871177454256097490213d0
1353
1354      rkB(1)   = .2113248654051871177454256097490213d0
1355      rkB(2)   = .5773502691896257645091487805019573d0
1356      rkB(3)   = .2113248654051871177454256097490213d0
1357
1358      rkBhat(1)= .2113248654051871177454256097490213d0
1359      rkBhat(2)= .6477918909913548037576239837516312d0
1360      rkBhat(3)= .1408832436034580784969504064993475d0
1361
1362      rkC(1)   = .2113248654051871177454256097490213d0
1363      rkC(2)   = .4226497308103742354908512194980427d0
1364      rkC(3)   = 1.d0
1365
1366! Ynew = Yold + h*Sum_i {rkB_i*k_i} = Yold + Sum_i {rkD_i*Z_i}
1367      rkD(1)   = 0.d0
1368      rkD(2)   = 0.d0
1369      rkD(3)   = 1.d0
1370
1371! Err = h * Sum_i {(rkB_i-rkBhat_i)*k_i} = Sum_i {rkE_i*Z_i}
1372      rkE(1)   =  0.9106836025229590978424821138352906d0
1373      rkE(2)   = -1.244016935856292431175815447168624d0
1374      rkE(3)   =  0.3333333333333333333333333333333333d0
1375
1376! Local order of Err estimate
1377      rkElo    = 2
1378
1379! h*Sum_j {rkA_ij*k_j} = Sum_j {rkTheta_ij*Z_j}
1380      rkTheta(2,1) =  1.0d0
1381      rkTheta(3,1) = -1.732050807568877293527446341505872d0
1382      rkTheta(3,2) =  2.732050807568877293527446341505872d0
1383
1384! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j}
1385      rkAlpha(2,1) =   2.0d0
1386      rkAlpha(3,1) = -12.92820323027550917410978536602349d0
1387      rkAlpha(3,2) =   8.83012701892219323381861585376468d0
1388
1389      END SUBROUTINE Sdirk3a
1390
1391!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1392   END SUBROUTINE SdirkTLM ! AND ALL ITS INTERNAL PROCEDURES
1393!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1394
1395
1396!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1397   SUBROUTINE FUN_CHEM( T, Y, P )
1398!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1399
1400      USE KPP_ROOT_Parameters, ONLY: NVAR
1401      USE KPP_ROOT_Global, ONLY: TIME, FIX, RCONST
1402      USE KPP_ROOT_Function, ONLY: Fun
1403      USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO
1404
1405      KPP_REAL :: T, Told
1406      KPP_REAL :: Y(NVAR), P(NVAR)
1407     
1408      Told = TIME
1409      TIME = T
1410      CALL Update_SUN()
1411      CALL Update_RCONST()
1412     
1413      CALL Fun( Y, FIX, RCONST, P )
1414     
1415      TIME = Told
1416     
1417    END SUBROUTINE FUN_CHEM
1418
1419
1420!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1421    SUBROUTINE JAC_CHEM( T, Y, JV )
1422!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1423
1424
1425      USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO
1426      USE KPP_ROOT_Global, ONLY: TIME, FIX, RCONST
1427      USE KPP_ROOT_Jacobian, ONLY: Jac_SP,LU_IROW,LU_ICOL
1428      USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO
1429 
1430      KPP_REAL :: T, Told
1431      KPP_REAL :: Y(NVAR)
1432#ifdef FULL_ALGEBRA
1433      KPP_REAL :: JS(LU_NONZERO), JV(NVAR,NVAR)
1434      INTEGER :: i, j
1435#else
1436      KPP_REAL :: JV(LU_NONZERO)
1437#endif
1438 
1439      Told = TIME
1440      TIME = T
1441      CALL Update_SUN()
1442      CALL Update_RCONST()
1443
1444#ifdef FULL_ALGEBRA
1445      CALL Jac_SP(Y, FIX, RCONST, JS)
1446      DO j=1,NVAR
1447         DO i=1,NVAR
1448           JV(i,j) = 0.0D0
1449         END DO
1450      END DO
1451      DO i=1,LU_NONZERO
1452         JV(LU_IROW(i),LU_ICOL(i)) = JS(i)
1453      END DO
1454#else
1455      CALL Jac_SP(Y, FIX, RCONST, JV)
1456#endif
1457      TIME = Told
1458
1459    END SUBROUTINE JAC_CHEM
1460
1461END MODULE KPP_ROOT_Integrator
Note: See TracBrowser for help on using the repository browser.