source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/runge_kutta.f90 @ 3523

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

Merge of branch palm4u into trunk

File size: 68.0 KB
Line 
1!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
2!  RungeKutta - Fully Implicit 3-stage Runge-Kutta methods based on:      !
3!          * Radau-2A   quadrature (order 5)                              !
4!          * Radau-1A   quadrature (order 5)                              !
5!          * Lobatto-3C quadrature (order 4)                              !
6!          * Gauss      quadrature (order 6)                              !
7!  By default the code employs the KPP sparse linear algebra routines     !
8!  Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK)        !
9!                                                                         !
10!    (C)  Adrian Sandu, August 2005                                       !
11!    Virginia Polytechnic Institute and State University                  !
12!    Contact: sandu@cs.vt.edu                                             !
13!    Revised by Philipp Miehe and Adrian Sandu, May 2006                  !
14!    This implementation is part of KPP - the Kinetic PreProcessor        !
15!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
16
17MODULE KPP_ROOT_Integrator
18
19  USE KPP_ROOT_Precision
20  USE KPP_ROOT_Parameters, ONLY: NVAR, NSPEC, NFIX, LU_NONZERO
21  USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME
22  USE KPP_ROOT_Jacobian,   ONLY: LU_DIAG
23  USE KPP_ROOT_LinearAlgebra
24
25  IMPLICIT NONE
26  PUBLIC
27  SAVE
28
29!~~~>  Statistics on the work performed by the Runge-Kutta method
30  INTEGER, PARAMETER :: Nfun=1, Njac=2, Nstp=3, Nacc=4, &
31    Nrej=5, Ndec=6, Nsol=7, Nsng=8, Ntexit=1, Nhacc=2, Nhnew=3
32 
33CONTAINS
34
35  ! **************************************************************************
36
37  SUBROUTINE INTEGRATE( TIN, TOUT, &
38    ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U, IERR_U )
39
40    USE KPP_ROOT_Parameters, ONLY: NVAR
41    USE KPP_ROOT_Global,     ONLY: ATOL,RTOL,VAR
42
43    IMPLICIT NONE
44
45    KPP_REAL :: TIN  ! TIN - Start Time
46    KPP_REAL :: TOUT ! TOUT - End Time
47    INTEGER,       INTENT(IN),  OPTIONAL :: ICNTRL_U(20)
48    KPP_REAL, INTENT(IN),  OPTIONAL :: RCNTRL_U(20)
49    INTEGER,       INTENT(OUT), OPTIONAL :: ISTATUS_U(20)
50    KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20)
51    INTEGER,       INTENT(OUT), OPTIONAL :: IERR_U
52
53    INTEGER :: IERR
54
55    KPP_REAL :: RCNTRL(20), RSTATUS(20), T1, T2
56    INTEGER :: ICNTRL(20), ISTATUS(20)
57    INTEGER, SAVE :: Ntotal = 0
58
59    RCNTRL(1:20) = 0.0_dp
60    ICNTRL(1:20) = 0
61
62    !~~~> fine-tune the integrator:
63    ICNTRL(2)  = 0   ! 0=vector tolerances, 1=scalar tolerances
64    ICNTRL(5)  = 8   ! Max no. of Newton iterations
65    ICNTRL(6)  = 0   ! Starting values for Newton are interpolated (0) or zero (1)
66    ICNTRL(10) = 1   ! 0 - classic or 1 - SDIRK error estimation
67    ICNTRL(11) = 0   ! Gustaffson (0) or classic(1) controller
68
69    !~~~> if optional parameters are given, and if they are >0,
70    !     then use them to 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   
79    T1 = TIN; T2 = TOUT
80    CALL RungeKutta(  NVAR, T1, T2, VAR, RTOL, ATOL, &
81                      RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR  )
82
83    Ntotal = Ntotal + ISTATUS(Nstp)
84    PRINT*,'NSTEPS=',ISTATUS(Nstp),' (',Ntotal,')','  O3=', VAR(ind_O3)
85
86    ! if optional parameters are given for output
87    ! use them to store information in them
88    IF (PRESENT(ISTATUS_U)) ISTATUS_U(:) = ISTATUS(:)
89    IF (PRESENT(RSTATUS_U)) RSTATUS_U(:) = RSTATUS(:)
90    IF (PRESENT(IERR_U)) IERR_U = IERR
91
92    IF (IERR < 0) THEN
93      PRINT *,'Runge-Kutta: Unsuccessful exit at T=', TIN,' (IERR=',IERR,')'
94    ENDIF
95
96  END SUBROUTINE INTEGRATE
97
98
99!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
100  SUBROUTINE RungeKutta( N,T,Tend,Y,RelTol,AbsTol,    &
101                         RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR )
102!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
103!
104!  This implementation is based on the book and the code Radau5:
105!
106!         E. HAIRER AND G. WANNER
107!         "SOLVING ORDINARY DIFFERENTIAL EQUATIONS II.
108!              STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS."
109!         SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS 14,
110!         SPRINGER-VERLAG (1991)
111!
112!         UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
113!         CH-1211 GENEVE 24, SWITZERLAND
114!         E-MAIL:  HAIRER@DIVSUN.UNIGE.CH,  WANNER@DIVSUN.UNIGE.CH
115!
116!   Methods:
117!          * Radau-2A   quadrature (order 5)                             
118!          * Radau-1A   quadrature (order 5)                             
119!          * Lobatto-3C quadrature (order 4)                             
120!          * Gauss      quadrature (order 6)                             
121!                                                                         
122!   (C)  Adrian Sandu, August 2005                                       
123!   Virginia Polytechnic Institute and State University                 
124!   Contact: sandu@cs.vt.edu                                             
125!   Revised by Philipp Miehe and Adrian Sandu, May 2006                 
126!   This implementation is part of KPP - the Kinetic PreProcessor       
127!
128!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
129!
130!~~~>   INPUT ARGUMENTS:
131!       ----------------
132!
133!    Note: For input parameters equal to zero the default values of the
134!          corresponding variables are used.
135!
136!     N           Dimension of the system
137!     T           Initial time value
138!
139!     Tend        Final T value (Tend-T may be positive or negative)
140!
141!     Y(N)        Initial values for Y
142!
143!     RelTol,AbsTol   Relative and absolute error tolerances.
144!          for ICNTRL(2) = 0: AbsTol, RelTol are N-dimensional vectors
145!                        = 1: AbsTol, RelTol are scalars
146!
147!~~~>  Integer input parameters:
148
149!    ICNTRL(1) = not used
150!
151!    ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors
152!              = 1: AbsTol, RelTol are scalars
153!
154!    ICNTRL(3) = RK method selection       
155!              = 1:  Radau-2A    (the default)
156!              = 2:  Lobatto-3C
157!              = 3:  Gauss
158!              = 4:  Radau-1A
159!              = 5:  Lobatto-3A (not yet implemented)
160!
161!    ICNTRL(4)  -> maximum number of integration steps
162!        For ICNTRL(4)=0 the default value of 10000 is used
163!
164!    ICNTRL(5)  -> maximum number of Newton iterations
165!        For ICNTRL(5)=0 the default value of 8 is used
166!
167!    ICNTRL(6)  -> starting values of Newton iterations:
168!        ICNTRL(6)=0 : starting values are obtained from
169!                      the extrapolated collocation solution
170!                      (the default)
171!        ICNTRL(6)=1 : starting values are zero
172!
173!    ICNTRL(10) -> switch for error estimation strategy
174!               ICNTRL(10) = 0: one additional stage at c=0,
175!                               see Hairer (default)
176!               ICNTRL(10) = 1: two additional stages at c=0
177!                               and SDIRK at c=1, stiffly accurate
178!
179!    ICNTRL(11) -> switch for step size strategy
180!              ICNTRL(11)=0:  mod. predictive controller (Gustafsson, default)
181!              ICNTRL(11)=1:  classical step size control
182!              the choice 1 seems to produce safer results;
183!              for simple problems, the choice 2 produces
184!              often slightly faster runs
185!
186!~~~>  Real input parameters:
187!
188!    RCNTRL(1)  -> Hmin, lower bound for the integration step size
189!                  (highly recommended to keep Hmin = ZERO, the default)
190!
191!    RCNTRL(2)  -> Hmax, upper bound for the integration step size
192!
193!    RCNTRL(3)  -> Hstart, the starting step size
194!
195!    RCNTRL(4)  -> FacMin, lower bound on step decrease factor (default=0.2)
196!
197!    RCNTRL(5)  -> FacMax, upper bound on step increase factor (default=6)
198!
199!    RCNTRL(6)  -> FacRej, step decrease factor after multiple rejections
200!                 (default=0.1)
201!
202!    RCNTRL(7)  -> FacSafe, by which the new step is slightly smaller
203!                  than the predicted value  (default=0.9)
204!
205!    RCNTRL(8)  -> ThetaMin. If Newton convergence rate smaller
206!                  than ThetaMin the Jacobian is not recomputed;
207!                  (default=0.001)
208!
209!    RCNTRL(9)  -> NewtonTol, stopping criterion for Newton's method
210!                  (default=0.03)
211!
212!    RCNTRL(10) -> Qmin
213!
214!    RCNTRL(11) -> Qmax. If Qmin < Hnew/Hold < Qmax, then the
215!                  step size is kept constant and the LU factorization
216!                  reused (default Qmin=1, Qmax=1.2)
217!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
218!
219!
220!    OUTPUT ARGUMENTS:
221!    -----------------
222!
223!    T           -> T value for which the solution has been computed
224!                     (after successful return T=Tend).
225!
226!    Y(N)        -> Numerical solution at T
227!
228!    IERR        -> Reports on successfulness upon return:
229!                    = 1 for success
230!                    < 0 for error (value equals error code)
231!
232!    ISTATUS(1)  -> No. of function calls
233!    ISTATUS(2)  -> No. of Jacobian calls
234!    ISTATUS(3)  -> No. of steps
235!    ISTATUS(4)  -> No. of accepted steps
236!    ISTATUS(5)  -> No. of rejected steps (except at very beginning)
237!    ISTATUS(6)  -> No. of LU decompositions
238!    ISTATUS(7)  -> No. of forward/backward substitutions
239!    ISTATUS(8)  -> No. of singular matrix decompositions
240!
241!    RSTATUS(1)  -> Texit, the time corresponding to the
242!                     computed Y upon return
243!    RSTATUS(2)  -> Hexit, last accepted step before exit
244!    RSTATUS(3)  -> Hnew, last predicted step (not yet taken)
245!                   For multiple restarts, use Hnew as Hstart
246!                     in the subsequent run
247!
248!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
249
250      IMPLICIT NONE
251     
252      INTEGER :: N
253      KPP_REAL :: Y(N),AbsTol(N),RelTol(N),RCNTRL(20),RSTATUS(20)
254      INTEGER :: ICNTRL(20), ISTATUS(20)
255      LOGICAL :: StartNewton, Gustafsson, SdirkError
256      INTEGER :: IERR, ITOL
257      KPP_REAL :: T,Tend
258
259      !~~~> Control arguments
260      INTEGER :: Max_no_steps, NewtonMaxit, rkMethod
261      KPP_REAL :: Hmin,Hmax,Hstart,Qmin,Qmax
262      KPP_REAL :: Roundoff, ThetaMin, NewtonTol
263      KPP_REAL :: FacSafe,FacMin,FacMax,FacRej
264      ! Runge-Kutta method parameters
265      INTEGER, PARAMETER :: R2A=1, R1A=2, L3C=3, GAU=4, L3A=5
266      KPP_REAL :: rkT(3,3), rkTinv(3,3), rkTinvAinv(3,3), rkAinvT(3,3), &
267                       rkA(0:3,0:3), rkB(0:3),  rkC(0:3), rkD(0:3), rkE(0:3), &
268                       rkBgam(0:4), rkBhat(0:4), rkTheta(0:3), rkF(0:4),      &
269                       rkGamma,  rkAlpha, rkBeta, rkELO
270       !~~~> Local variables
271      INTEGER :: i
272      KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
273   
274!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
275!        SETTING THE PARAMETERS
276!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
277      IERR = 0
278      ISTATUS(1:20) = 0
279      RSTATUS(1:20) = ZERO
280       
281!~~~> ICNTRL(1) - autonomous system - not used       
282!~~~> ITOL: 1 for vector and 0 for scalar AbsTol/RelTol
283      IF (ICNTRL(2) == 0) THEN
284         ITOL = 1
285      ELSE
286         ITOL = 0
287      END IF
288!~~~> Error control selection 
289      IF (ICNTRL(10) == 0) THEN
290         SdirkError = .FALSE.
291      ELSE
292         SdirkError = .TRUE.
293      END IF     
294!~~~> Method selection 
295      SELECT CASE (ICNTRL(3))     
296      CASE (0,1)
297         CALL Radau2A_Coefficients
298      CASE (2)
299         CALL Lobatto3C_Coefficients
300      CASE (3)
301         CALL Gauss_Coefficients
302      CASE (4)
303         CALL Radau1A_Coefficients
304      CASE (5)
305         CALL Lobatto3A_Coefficients
306      CASE DEFAULT
307         WRITE(6,*) 'ICNTRL(3)=',ICNTRL(3)
308         CALL RK_ErrorMsg(-13,T,ZERO,IERR)
309      END SELECT
310!~~~> Max_no_steps: the maximal number of time steps
311      IF (ICNTRL(4) == 0) THEN
312         Max_no_steps = 200000
313      ELSE
314         Max_no_steps=ICNTRL(4)
315         IF (Max_no_steps <= 0) THEN
316            WRITE(6,*) 'ICNTRL(4)=',ICNTRL(4)
317            CALL RK_ErrorMsg(-1,T,ZERO,IERR)
318         END IF
319      END IF
320!~~~> NewtonMaxit    maximal number of Newton iterations
321      IF (ICNTRL(5) == 0) THEN
322         NewtonMaxit = 8
323      ELSE
324         NewtonMaxit=ICNTRL(5)
325         IF (NewtonMaxit <= 0) THEN
326            WRITE(6,*) 'ICNTRL(5)=',ICNTRL(5)
327            CALL RK_ErrorMsg(-2,T,ZERO,IERR)
328          END IF
329      END IF
330!~~~> StartNewton:  Use extrapolation for starting values of Newton iterations
331      IF (ICNTRL(6) == 0) THEN
332         StartNewton = .TRUE.
333      ELSE
334         StartNewton = .FALSE.
335      END IF     
336!~~~> Gustafsson: step size controller
337      IF(ICNTRL(11) == 0)THEN
338         Gustafsson = .TRUE.
339      ELSE
340         Gustafsson = .FALSE.
341      END IF
342
343!~~~> Roundoff: smallest number s.t. 1.0 + Roundoff > 1.0
344      Roundoff=WLAMCH('E');
345
346!~~~> Hmin = minimal step size
347      IF (RCNTRL(1) == ZERO) THEN
348         Hmin = ZERO
349      ELSE
350         Hmin = MIN(ABS(RCNTRL(1)),ABS(Tend-T))
351      END IF
352!~~~> Hmax = maximal step size
353      IF (RCNTRL(2) == ZERO) THEN
354         Hmax = ABS(Tend-T)
355      ELSE
356         Hmax = MIN(ABS(RCNTRL(2)),ABS(Tend-T))
357      END IF
358!~~~> Hstart = starting step size
359      IF (RCNTRL(3) == ZERO) THEN
360         Hstart = ZERO
361      ELSE
362         Hstart = MIN(ABS(RCNTRL(3)),ABS(Tend-T))
363      END IF
364!~~~> FacMin: lower bound on step decrease factor
365      IF(RCNTRL(4) == ZERO)THEN
366         FacMin = 0.2d0
367      ELSE
368         FacMin = RCNTRL(4)
369      END IF
370!~~~> FacMax: upper bound on step increase factor
371      IF(RCNTRL(5) == ZERO)THEN
372         FacMax = 8.D0
373      ELSE
374         FacMax = RCNTRL(5)
375      END IF
376!~~~> FacRej: step decrease factor after 2 consecutive rejections
377      IF(RCNTRL(6) == ZERO)THEN
378         FacRej = 0.1d0
379      ELSE
380         FacRej = RCNTRL(6)
381      END IF
382!~~~> FacSafe:  by which the new step is slightly smaller
383!               than the predicted value
384      IF (RCNTRL(7) == ZERO) THEN
385         FacSafe=0.9d0
386      ELSE
387         FacSafe=RCNTRL(7)
388      END IF
389      IF ( (FacMax < ONE) .OR. (FacMin > ONE) .OR. &
390           (FacSafe <= 1.0d-3) .OR. (FacSafe >= ONE) ) THEN
391            WRITE(6,*)'RCNTRL(4:7)=',RCNTRL(4:7)
392            CALL RK_ErrorMsg(-4,T,ZERO,IERR)
393      END IF
394
395!~~~> ThetaMin:  decides whether the Jacobian should be recomputed
396      IF (RCNTRL(8) == ZERO) THEN
397         ThetaMin = 1.0d-3
398      ELSE
399         ThetaMin=RCNTRL(8)
400         IF (ThetaMin <= 0.0d0 .OR. ThetaMin >= 1.0d0) THEN
401            WRITE(6,*) 'RCNTRL(8)=', RCNTRL(8)
402            CALL RK_ErrorMsg(-5,T,ZERO,IERR)
403         END IF
404      END IF
405!~~~> NewtonTol:  stopping crierion for Newton's method
406      IF (RCNTRL(9) == ZERO) THEN
407         NewtonTol = 3.0d-2
408      ELSE
409         NewtonTol = RCNTRL(9)
410         IF (NewtonTol <= Roundoff) THEN
411            WRITE(6,*) 'RCNTRL(9)=',RCNTRL(9)
412            CALL RK_ErrorMsg(-6,T,ZERO,IERR)
413         END IF
414      END IF
415!~~~> Qmin AND Qmax: IF Qmin < Hnew/Hold < Qmax then step size = const.
416      IF (RCNTRL(10) == ZERO) THEN
417         Qmin=1.D0
418      ELSE
419         Qmin=RCNTRL(10)
420      END IF
421      IF (RCNTRL(11) == ZERO) THEN
422         Qmax=1.2D0
423      ELSE
424         Qmax=RCNTRL(11)
425      END IF
426      IF (Qmin > ONE .OR. Qmax < ONE) THEN
427         WRITE(6,*) 'RCNTRL(10:11)=',Qmin,Qmax
428         CALL RK_ErrorMsg(-7,T,ZERO,IERR)
429      END IF
430!~~~> Check if tolerances are reasonable
431      IF (ITOL == 0) THEN
432          IF (AbsTol(1) <= ZERO.OR.RelTol(1) <= 10.d0*Roundoff) THEN
433              WRITE (6,*) 'AbsTol/RelTol=',AbsTol,RelTol 
434              CALL RK_ErrorMsg(-8,T,ZERO,IERR)
435          END IF
436      ELSE
437          DO i=1,N
438          IF (AbsTol(i) <= ZERO.OR.RelTol(i) <= 10.d0*Roundoff) THEN
439              WRITE (6,*) 'AbsTol/RelTol(',i,')=',AbsTol(i),RelTol(i)
440              CALL RK_ErrorMsg(-8,T,ZERO,IERR)
441          END IF
442          END DO
443      END IF
444
445!~~~> Parameters are wrong
446      IF (IERR < 0) RETURN
447
448!~~~> Call the core method
449      CALL RK_Integrator( N,T,Tend,Y,IERR )
450
451   CONTAINS ! Internal procedures to RungeKutta
452
453!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
454   SUBROUTINE RK_Integrator( N,T,Tend,Y,IERR )
455!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
456
457      IMPLICIT NONE
458!~~~> Arguments
459      INTEGER,  INTENT(IN)         :: N
460      KPP_REAL, INTENT(IN)    :: Tend
461      KPP_REAL, INTENT(INOUT) :: T, Y(NVAR)
462      INTEGER,  INTENT(OUT)        :: IERR
463
464!~~~> Local variables
465#ifdef FULL_ALGEBRA
466      KPP_REAL    :: FJAC(NVAR,NVAR), E1(NVAR,NVAR)
467      COMPLEX(kind=dp) :: E2(NVAR,NVAR)   
468#else
469      KPP_REAL    :: FJAC(LU_NONZERO), E1(LU_NONZERO)
470      COMPLEX(kind=dp) :: E2(LU_NONZERO)   
471#endif               
472      KPP_REAL, DIMENSION(NVAR) :: Z1,Z2,Z3,Z4,SCAL,DZ1,DZ2,DZ3,DZ4, &
473                                G,TMP,F0
474      KPP_REAL  :: CONT(NVAR,3), Tdirection,  H, Hacc, Hnew, Hold, Fac, &
475                 FacGus, Theta, Err, ErrOld, NewtonRate, NewtonIncrement,    &
476                 Hratio, Qnewton, NewtonPredictedErr,NewtonIncrementOld, ThetaSD
477      INTEGER :: IP1(NVAR),IP2(NVAR),NewtonIter, ISING, Nconsecutive
478      LOGICAL :: Reject, FirstStep, SkipJac, NewtonDone, SkipLU
479     
480           
481!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
482!~~~>  Initial setting
483!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
484      Tdirection = SIGN(ONE,Tend-T)
485      H = MIN( MAX(ABS(Hmin),ABS(Hstart)) , Hmax )
486      IF (ABS(H) <= 10.d0*Roundoff) H = 1.0d-6
487      H = SIGN(H,Tdirection)
488      Hold      = H
489      Reject    = .FALSE.
490      FirstStep = .TRUE.
491      SkipJac   = .FALSE.
492      SkipLU    = .FALSE.
493      IF ((T+H*1.0001D0-Tend)*Tdirection >= ZERO) THEN
494         H = Tend-T
495      END IF
496      Nconsecutive = 0
497      CALL RK_ErrorScale(N,ITOL,AbsTol,RelTol,Y,SCAL)
498
499!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
500!~~~>  Time loop begins
501!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
502Tloop: DO WHILE ( (Tend-T)*Tdirection - Roundoff > ZERO )
503
504      !IF ( .NOT.Reject ) THEN     
505         CALL FUN_CHEM(T,Y,F0)
506         ISTATUS(Nfun) = ISTATUS(Nfun) + 1
507      !END IF   
508
509      IF ( .NOT.SkipLU ) THEN ! This time around skip the Jac update and LU
510        !~~~> Compute the Jacobian matrix
511        IF ( .NOT.SkipJac ) THEN
512          CALL JAC_CHEM(T,Y,FJAC)
513          ISTATUS(Njac) = ISTATUS(Njac) + 1
514        END IF
515        !~~~> Compute the matrices E1 and E2 and their decompositions
516        CALL RK_Decomp(N,H,FJAC,E1,IP1,E2,IP2,ISING)
517        IF (ISING /= 0) THEN
518          ISTATUS(Nsng) = ISTATUS(Nsng) + 1; Nconsecutive = Nconsecutive + 1
519          IF (Nconsecutive >= 5) THEN
520            CALL RK_ErrorMsg(-12,T,H,IERR); RETURN
521          END IF
522          H=H*0.5d0; Reject=.TRUE.; SkipJac = .TRUE.;  SkipLU = .FALSE.
523          CYCLE Tloop
524        ELSE
525          Nconsecutive = 0   
526        END IF   
527      END IF ! SkipLU
528   
529      ISTATUS(Nstp) = ISTATUS(Nstp) + 1
530      IF (ISTATUS(Nstp) > Max_no_steps) THEN
531        PRINT*,'Max number of time steps is ',Max_no_steps
532        CALL RK_ErrorMsg(-9,T,H,IERR); RETURN
533      END IF
534      IF (0.1D0*ABS(H) <= ABS(T)*Roundoff)  THEN
535        CALL RK_ErrorMsg(-10,T,H,IERR); RETURN
536      END IF
537     
538!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
539!~~~>  Loop for the simplified Newton iterations
540!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
541           
542      !~~~>  Starting values for Newton iteration
543      IF ( FirstStep .OR. (.NOT.StartNewton) ) THEN
544         CALL Set2zero(N,Z1)
545         CALL Set2zero(N,Z2)
546         CALL Set2zero(N,Z3)
547      ELSE
548         ! Evaluate quadratic polynomial
549         CALL RK_Interpolate('eval',N,H,Hold,Z1,Z2,Z3,CONT)
550      END IF
551     
552      !~~~>  Initializations for Newton iteration
553      NewtonDone = .FALSE.
554      Fac = 0.5d0 ! Step reduction if too many iterations
555     
556NewtonLoop:DO  NewtonIter = 1, NewtonMaxit
557 
558            !~~~> Prepare the right-hand side
559            CALL RK_PrepareRHS(N,T,H,Y,F0,Z1,Z2,Z3,DZ1,DZ2,DZ3)
560           
561            !~~~> Solve the linear systems
562            CALL RK_Solve( N,H,E1,IP1,E2,IP2,DZ1,DZ2,DZ3,ISING )
563           
564            NewtonIncrement = SQRT( ( RK_ErrorNorm(N,SCAL,DZ1)**2 + &
565                                RK_ErrorNorm(N,SCAL,DZ2)**2 + &
566                                RK_ErrorNorm(N,SCAL,DZ3)**2 )/3.0d0 )
567           
568            IF ( NewtonIter == 1 ) THEN
569                Theta      = ABS(ThetaMin)
570                NewtonRate = 2.0d0 
571            ELSE
572                Theta = NewtonIncrement/NewtonIncrementOld
573                IF (Theta < 0.99d0) THEN
574                    NewtonRate = Theta/(ONE-Theta)
575                ELSE ! Non-convergence of Newton: Theta too large
576                    EXIT NewtonLoop
577                END IF
578                IF ( NewtonIter < NewtonMaxit ) THEN 
579                  ! Predict error at the end of Newton process
580                  NewtonPredictedErr = NewtonIncrement &
581                               *Theta**(NewtonMaxit-NewtonIter)/(ONE-Theta)
582                  IF (NewtonPredictedErr >= NewtonTol) THEN
583                    ! Non-convergence of Newton: predicted error too large
584                    Qnewton = MIN(10.0d0,NewtonPredictedErr/NewtonTol)
585                    Fac=0.8d0*Qnewton**(-ONE/(1+NewtonMaxit-NewtonIter))
586                    EXIT NewtonLoop
587                  END IF
588                END IF
589            END IF
590
591            NewtonIncrementOld = MAX(NewtonIncrement,Roundoff) 
592            ! Update solution
593            CALL WAXPY(N,-ONE,DZ1,1,Z1,1) ! Z1 <- Z1 - DZ1
594            CALL WAXPY(N,-ONE,DZ2,1,Z2,1) ! Z2 <- Z2 - DZ2
595            CALL WAXPY(N,-ONE,DZ3,1,Z3,1) ! Z3 <- Z3 - DZ3
596           
597            ! Check error in Newton iterations
598            NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol)
599            IF (NewtonDone) EXIT NewtonLoop
600            IF (NewtonIter == NewtonMaxit) THEN
601                PRINT*, 'Slow or no convergence in Newton Iteration: Max no. of', &
602                        'Newton iterations reached'
603            END IF
604           
605      END DO NewtonLoop
606           
607      IF (.NOT.NewtonDone) THEN
608           !CALL RK_ErrorMsg(-12,T,H,IERR);
609           H = Fac*H; Reject=.TRUE.; SkipJac = .TRUE.;  SkipLU = .FALSE.
610           CYCLE Tloop
611      END IF
612
613!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
614!~~~> SDIRK Stage
615!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
616      IF (SdirkError) THEN
617
618!~~~>  Starting values for Newton iterations
619       Z4(1:N) = Z3(1:N)
620       
621!~~~>   Prepare the loop-independent part of the right-hand side
622!       G = H*rkBgam(0)*F0 + rkTheta(1)*Z1 + rkTheta(2)*Z2 + rkTheta(3)*Z3
623       CALL Set2Zero(N, G)
624       IF (rkMethod/=L3A) CALL WAXPY(N,rkBgam(0)*H, F0,1,G,1) 
625       CALL WAXPY(N,rkTheta(1),Z1,1,G,1)
626       CALL WAXPY(N,rkTheta(2),Z2,1,G,1)
627       CALL WAXPY(N,rkTheta(3),Z3,1,G,1)
628
629       !~~~>  Initializations for Newton iteration
630       NewtonDone = .FALSE.
631       Fac = 0.5d0 ! Step reduction factor if too many iterations
632           
633SDNewtonLoop:DO NewtonIter = 1, NewtonMaxit
634
635!~~~>   Prepare the loop-dependent part of the right-hand side
636            CALL WADD(N,Y,Z4,TMP)         ! TMP <- Y + Z4
637            CALL FUN_CHEM(T+H,TMP,DZ4)    ! DZ4 <- Fun(Y+Z4)         
638            ISTATUS(Nfun) = ISTATUS(Nfun) + 1
639!            DZ4(1:N) = (G(1:N)-Z4(1:N))*(rkGamma/H) + DZ4(1:N)
640            CALL WAXPY (N, -ONE*rkGamma/H, Z4, 1, DZ4, 1)
641            CALL WAXPY (N, rkGamma/H, G,1, DZ4,1)
642
643!~~~>   Solve the linear system
644#ifdef FULL_ALGEBRA 
645            CALL DGETRS( 'N', N, 1, E1, N, IP1, DZ4, N, ISING )
646#else
647            CALL KppSolve(E1, DZ4)
648#endif
649           
650!~~~>   Check convergence of Newton iterations
651            NewtonIncrement = RK_ErrorNorm(N,SCAL,DZ4)
652            IF ( NewtonIter == 1 ) THEN
653                ThetaSD      = ABS(ThetaMin)
654                NewtonRate = 2.0d0 
655            ELSE
656                ThetaSD = NewtonIncrement/NewtonIncrementOld
657                IF (ThetaSD < 0.99d0) THEN
658                    NewtonRate = ThetaSD/(ONE-ThetaSD)
659                    ! Predict error at the end of Newton process
660                    NewtonPredictedErr = NewtonIncrement &
661                               *ThetaSD**(NewtonMaxit-NewtonIter)/(ONE-ThetaSD)
662                    IF (NewtonPredictedErr >= NewtonTol) THEN
663                      ! Non-convergence of Newton: predicted error too large
664                      !PRINT*,'Error too large: ', NewtonPredictedErr
665                      Qnewton = MIN(10.0d0,NewtonPredictedErr/NewtonTol)
666                      Fac = 0.8d0*Qnewton**(-ONE/(1+NewtonMaxit-NewtonIter))
667                      EXIT SDNewtonLoop
668                    END IF
669                ELSE ! Non-convergence of Newton: Theta too large
670                    !PRINT*,'Theta too large: ',ThetaSD
671                    EXIT SDNewtonLoop
672                END IF
673            END IF
674            NewtonIncrementOld = NewtonIncrement
675            ! Update solution: Z4 <-- Z4 + DZ4
676            CALL WAXPY(N,ONE,DZ4,1,Z4,1) 
677           
678            ! Check error in Newton iterations
679            NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol)
680            IF (NewtonDone) EXIT SDNewtonLoop
681           
682            END DO SDNewtonLoop
683           
684            IF (.NOT.NewtonDone) THEN
685                H = Fac*H; Reject=.TRUE.; SkipJac = .TRUE.;  SkipLU = .FALSE.
686                CYCLE Tloop
687            END IF
688      END IF
689!~~~>  End of implified SDIRK Newton iterations
690
691!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
692!~~~> Error estimation
693!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
694      IF (SdirkError) THEN
695         CALL Set2Zero(N, DZ4)
696         IF (rkMethod==L3A) THEN
697           DZ4(1:N) = H*rkF(0)*F0(1:N)
698           IF (rkF(1) /= ZERO)  CALL WAXPY(N, rkF(1), Z1, 1, DZ4, 1)
699           IF (rkF(2) /= ZERO)  CALL WAXPY(N, rkF(2), Z2, 1, DZ4, 1)
700           IF (rkF(3) /= ZERO)  CALL WAXPY(N, rkF(3), Z3, 1, DZ4, 1)
701           TMP = Y + Z4
702           CALL FUN_CHEM(T+H,TMP,DZ1)
703           CALL WAXPY(N, H*rkBgam(4), DZ1, 1, DZ4, 1)
704         ELSE
705!         DZ4(1:N) =  rkD(1)*Z1 + rkD(2)*Z2 + rkD(3)*Z3 - Z4   
706           IF (rkD(1) /= ZERO)  CALL WAXPY(N, rkD(1), Z1, 1, DZ4, 1)
707           IF (rkD(2) /= ZERO)  CALL WAXPY(N, rkD(2), Z2, 1, DZ4, 1)
708           IF (rkD(3) /= ZERO)  CALL WAXPY(N, rkD(3), Z3, 1, DZ4, 1)
709           CALL WAXPY(N, -ONE, Z4, 1, DZ4, 1)
710         END IF
711         Err = RK_ErrorNorm(N,SCAL,DZ4)   
712      ELSE
713         CALL  RK_ErrorEstimate(N,H,T,Y,F0, &
714               E1,IP1,Z1,Z2,Z3,SCAL,Err,FirstStep,Reject)
715      END IF
716
717!~~~> Computation of new step size Hnew
718      Fac  = Err**(-ONE/rkELO)*   &
719             MIN(FacSafe,(ONE+2*NewtonMaxit)/(NewtonIter+2*NewtonMaxit))
720      Fac  = MIN(FacMax,MAX(FacMin,Fac))
721      Hnew = Fac*H
722
723!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
724!~~~> Accept/reject step
725!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
726accept:IF (Err < ONE) THEN !~~~> STEP IS ACCEPTED
727         FirstStep=.FALSE.
728         ISTATUS(Nacc) = ISTATUS(Nacc) + 1
729         IF (Gustafsson) THEN
730            !~~~> Predictive controller of Gustafsson
731            IF (ISTATUS(Nacc) > 1) THEN
732               FacGus=FacSafe*(H/Hacc)*(Err**2/ErrOld)**(-0.25d0)
733               FacGus=MIN(FacMax,MAX(FacMin,FacGus))
734               Fac=MIN(Fac,FacGus)
735               Hnew = Fac*H
736            END IF
737            Hacc=H
738            ErrOld=MAX(1.0d-2,Err)
739         END IF
740         Hold = H
741         T = T+H 
742         ! Update solution: Y <- Y + sum(d_i Z_i)
743         IF (rkD(1) /= ZERO)  CALL WAXPY(N,rkD(1),Z1,1,Y,1)
744         IF (rkD(2) /= ZERO)  CALL WAXPY(N,rkD(2),Z2,1,Y,1)
745         IF (rkD(3) /= ZERO)  CALL WAXPY(N,rkD(3),Z3,1,Y,1)
746         ! Construct the solution quadratic interpolant Q(c_i) = Z_i, i=1:3
747         IF (StartNewton) CALL RK_Interpolate('make',N,H,Hold,Z1,Z2,Z3,CONT)
748         CALL RK_ErrorScale(N,ITOL,AbsTol,RelTol,Y,SCAL)
749         RSTATUS(Ntexit) = T
750         RSTATUS(Nhnew)  = Hnew
751         RSTATUS(Nhacc)  = H
752         Hnew = Tdirection*MIN( MAX(ABS(Hnew),Hmin) , Hmax )
753         IF (Reject) Hnew = Tdirection*MIN(ABS(Hnew),ABS(H))
754         Reject = .FALSE.
755         IF ((T+Hnew/Qmin-Tend)*Tdirection >=  ZERO) THEN
756            H = Tend-T
757         ELSE
758            Hratio=Hnew/H
759            ! Reuse the LU decomposition
760            SkipLU = (Theta<=ThetaMin) .AND. (Hratio>=Qmin) .AND. (Hratio<=Qmax)
761            IF (.NOT.SkipLU) H=Hnew
762         END IF
763         ! If convergence is fast enough, do not update Jacobian
764!         SkipJac = (Theta <= ThetaMin)
765         SkipJac  = .FALSE.
766
767      ELSE accept !~~~> Step is rejected
768         IF (FirstStep .OR. Reject) THEN
769             H = FacRej*H
770         ELSE
771             H = Hnew
772         END IF
773         Reject   = .TRUE.
774         SkipJac  = .TRUE.      ! Skip if rejected - Jac is independent of H
775         SkipLU   = .FALSE. 
776         IF (ISTATUS(Nacc) >= 1) ISTATUS(Nrej) = ISTATUS(Nrej) + 1
777      END IF accept
778     
779    END DO Tloop
780   
781    ! Successful exit
782    IERR = 1 
783
784 END SUBROUTINE RK_Integrator
785
786
787!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
788 SUBROUTINE RK_ErrorMsg(Code,T,H,IERR)
789!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
790!    Handles all error messages
791!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
792
793   IMPLICIT NONE
794   KPP_REAL, INTENT(IN) :: T, H
795   INTEGER, INTENT(IN)  :: Code
796   INTEGER, INTENT(OUT) :: IERR
797
798   IERR = Code
799   PRINT * , &
800     'Forced exit from RungeKutta due to the following error:'
801
802
803   SELECT CASE (Code)
804    CASE (-1)
805      PRINT * , '--> Improper value for maximal no of steps'
806    CASE (-2)
807      PRINT * , '--> Improper value for maximal no of Newton iterations'
808    CASE (-3)
809      PRINT * , '--> Hmin/Hmax/Hstart must be positive'
810    CASE (-4)
811      PRINT * , '--> Improper values for FacMin/FacMax/FacSafe/FacRej'
812    CASE (-5)
813      PRINT * , '--> Improper value for ThetaMin'
814    CASE (-6)
815      PRINT * , '--> Newton stopping tolerance too small'
816    CASE (-7)
817      PRINT * , '--> Improper values for Qmin, Qmax'
818    CASE (-8)
819      PRINT * , '--> Tolerances are too small'
820    CASE (-9)
821      PRINT * , '--> No of steps exceeds maximum bound'
822    CASE (-10)
823      PRINT * , '--> Step size too small: T + 10*H = T', &
824            ' or H < Roundoff'
825    CASE (-11)
826      PRINT * , '--> Matrix is repeatedly singular'
827    CASE (-12)
828      PRINT * , '--> Non-convergence of Newton iterations'
829    CASE (-13)
830      PRINT * , '--> Requested RK method not implemented'
831    CASE DEFAULT
832      PRINT *, 'Unknown Error code: ', Code
833   END SELECT
834
835   WRITE(6,FMT="(5X,'T=',E12.5,'  H=',E12.5)") T, H 
836
837 END SUBROUTINE RK_ErrorMsg
838
839
840!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
841 SUBROUTINE RK_ErrorScale(N,ITOL,AbsTol,RelTol,Y,SCAL)
842!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
843!    Handles all error messages
844!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
845   IMPLICIT NONE
846   INTEGER, INTENT(IN)  :: N, ITOL
847   KPP_REAL, INTENT(IN) :: AbsTol(*), RelTol(*), Y(N)
848   KPP_REAL, INTENT(OUT) :: SCAL(N)
849   INTEGER :: i
850   
851   IF (ITOL==0) THEN
852       DO i=1,N
853          SCAL(i)= ONE/(AbsTol(1)+RelTol(1)*ABS(Y(i)))
854       END DO
855   ELSE
856       DO i=1,N
857          SCAL(i)=ONE/(AbsTol(i)+RelTol(i)*ABS(Y(i)))
858       END DO
859   END IF
860     
861 END SUBROUTINE RK_ErrorScale
862
863
864!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
865!  SUBROUTINE RK_Transform(N,Tr,Z1,Z2,Z3,W1,W2,W3)
866!!~~~>                 W <-- Tr x Z
867!!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
868!      IMPLICIT NONE
869!      INTEGER :: N, i
870!      KPP_REAL :: Tr(3,3),Z1(N),Z2(N),Z3(N),W1(N),W2(N),W3(N)
871!      KPP_REAL :: x1, x2, x3
872!      DO i=1,N
873!          x1 = Z1(i); x2 = Z2(i); x3 = Z3(i)
874!          W1(i) = Tr(1,1)*x1 + Tr(1,2)*x2 + Tr(1,3)*x3
875!          W2(i) = Tr(2,1)*x1 + Tr(2,2)*x2 + Tr(2,3)*x3
876!          W3(i) = Tr(3,1)*x1 + Tr(3,2)*x2 + Tr(3,3)*x3
877!      END DO
878!  END SUBROUTINE RK_Transform
879 
880!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
881  SUBROUTINE RK_Interpolate(action,N,H,Hold,Z1,Z2,Z3,CONT)
882!~~~>   Constructs or evaluates a quadratic polynomial
883!         that interpolates the Z solution at current step
884!         and provides starting values for the next step
885!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
886      INTEGER :: N, i
887      KPP_REAL :: H,Hold,Z1(N),Z2(N),Z3(N),CONT(N,3)
888      KPP_REAL :: r, x1, x2, x3, den
889      CHARACTER(LEN=4) :: action
890       
891      SELECT CASE (action) 
892      CASE ('make')
893         ! Construct the solution quadratic interpolant Q(c_i) = Z_i, i=1:3
894         den = (rkC(3)-rkC(2))*(rkC(2)-rkC(1))*(rkC(1)-rkC(3))
895         DO i=1,N
896             CONT(i,1)=(-rkC(3)**2*rkC(2)*Z1(i)+Z3(i)*rkC(2)*rkC(1)**2 &
897                        +rkC(2)**2*rkC(3)*Z1(i)-rkC(2)**2*rkC(1)*Z3(i) &
898                        +rkC(3)**2*rkC(1)*Z2(i)-Z2(i)*rkC(3)*rkC(1)**2)&
899                        /den-Z3(i)
900             CONT(i,2)= -( rkC(1)**2*(Z3(i)-Z2(i)) + rkC(2)**2*(Z1(i)  &
901                          -Z3(i)) +rkC(3)**2*(Z2(i)-Z1(i)) )/den
902             CONT(i,3)= ( rkC(1)*(Z3(i)-Z2(i)) + rkC(2)*(Z1(i)-Z3(i))  &
903                           +rkC(3)*(Z2(i)-Z1(i)) )/den
904         END DO
905      CASE ('eval')
906          ! Evaluate quadratic polynomial
907          r = H/Hold
908         x1 = ONE + rkC(1)*r
909         x2 = ONE + rkC(2)*r
910         x3 = ONE + rkC(3)*r
911         DO i=1,N
912            Z1(i) = CONT(i,1)+x1*(CONT(i,2)+x1*CONT(i,3))
913            Z2(i) = CONT(i,1)+x2*(CONT(i,2)+x2*CONT(i,3))
914            Z3(i) = CONT(i,1)+x3*(CONT(i,2)+x3*CONT(i,3))
915         END DO
916       END SELECT   
917  END SUBROUTINE RK_Interpolate
918
919
920!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
921   SUBROUTINE RK_PrepareRHS(N,T,H,Y,F0,Z1,Z2,Z3,R1,R2,R3)
922!~~~> Prepare the right-hand side for Newton iterations
923!     R = Z - hA x F
924!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
925      IMPLICIT NONE
926     
927      INTEGER :: N
928      KPP_REAL :: T, H
929      KPP_REAL, DIMENSION(N) :: Y,Z1,Z2,Z3,F0,F,R1,R2,R3,TMP
930
931      CALL WCOPY(N,Z1,1,R1,1) ! R1 <- Z1
932      CALL WCOPY(N,Z2,1,R2,1) ! R2 <- Z2
933      CALL WCOPY(N,Z3,1,R3,1) ! R3 <- Z3
934
935      IF (rkMethod==L3A) THEN
936         CALL WAXPY(N,-H*rkA(1,0),F0,1,R1,1) ! R1 <- R1 - h*A_10*F0
937         CALL WAXPY(N,-H*rkA(2,0),F0,1,R2,1) ! R2 <- R2 - h*A_20*F0
938         CALL WAXPY(N,-H*rkA(3,0),F0,1,R3,1) ! R3 <- R3 - h*A_30*F0
939      END IF
940
941      CALL WADD(N,Y,Z1,TMP)              ! TMP <- Y + Z1
942      CALL FUN_CHEM(T+rkC(1)*H,TMP,F)    ! F1 <- Fun(Y+Z1)         
943      CALL WAXPY(N,-H*rkA(1,1),F,1,R1,1) ! R1 <- R1 - h*A_11*F1
944      CALL WAXPY(N,-H*rkA(2,1),F,1,R2,1) ! R2 <- R2 - h*A_21*F1
945      CALL WAXPY(N,-H*rkA(3,1),F,1,R3,1) ! R3 <- R3 - h*A_31*F1
946
947      CALL WADD(N,Y,Z2,TMP)              ! TMP <- Y + Z2
948      CALL FUN_CHEM(T+rkC(2)*H,TMP,F)    ! F2 <- Fun(Y+Z2)       
949      CALL WAXPY(N,-H*rkA(1,2),F,1,R1,1) ! R1 <- R1 - h*A_12*F2
950      CALL WAXPY(N,-H*rkA(2,2),F,1,R2,1) ! R2 <- R2 - h*A_22*F2
951      CALL WAXPY(N,-H*rkA(3,2),F,1,R3,1) ! R3 <- R3 - h*A_32*F2
952
953      CALL WADD(N,Y,Z3,TMP)              ! TMP <- Y + Z3
954      CALL FUN_CHEM(T+rkC(3)*H,TMP,F)    ! F3 <- Fun(Y+Z3)     
955      CALL WAXPY(N,-H*rkA(1,3),F,1,R1,1) ! R1 <- R1 - h*A_13*F3
956      CALL WAXPY(N,-H*rkA(2,3),F,1,R2,1) ! R2 <- R2 - h*A_23*F3
957      CALL WAXPY(N,-H*rkA(3,3),F,1,R3,1) ! R3 <- R3 - h*A_33*F3
958           
959  END SUBROUTINE RK_PrepareRHS
960 
961
962!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
963   SUBROUTINE RK_Decomp(N,H,FJAC,E1,IP1,E2,IP2,ISING)
964   !~~~> Compute the matrices E1 and E2 and their decompositions
965!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
966      IMPLICIT NONE
967     
968      INTEGER :: N, ISING
969      KPP_REAL    :: H, Alpha, Beta, Gamma
970#ifdef FULL_ALGEBRA     
971      KPP_REAL    :: FJAC(NVAR,NVAR),E1(NVAR,NVAR)
972      COMPLEX(kind=dp) :: E2(N,N)
973#else     
974      KPP_REAL    :: FJAC(LU_NONZERO),E1(LU_NONZERO)
975      COMPLEX(kind=dp) :: E2(LU_NONZERO)
976#endif     
977      INTEGER :: IP1(N), IP2(N), i, j
978     
979      Gamma = rkGamma/H
980      Alpha = rkAlpha/H
981      Beta  = rkBeta /H
982
983#ifdef FULL_ALGEBRA     
984      DO j=1,N
985         DO  i=1,N
986            E1(i,j)=-FJAC(i,j)
987         END DO
988         E1(j,j)=E1(j,j)+Gamma
989      END DO
990      CALL DGETRF(N,N,E1,N,IP1,ISING) 
991#else     
992      DO i=1,LU_NONZERO
993         E1(i)=-FJAC(i)
994      END DO
995      DO i=1,NVAR
996         j=LU_DIAG(i); E1(j)=E1(j)+Gamma
997      END DO
998      CALL KppDecomp(E1,ISING)
999#endif     
1000     
1001      IF (ISING /= 0) THEN
1002         ISTATUS(Ndec) = ISTATUS(Ndec) + 1
1003         RETURN
1004      END IF
1005     
1006#ifdef FULL_ALGEBRA     
1007      DO j=1,N
1008        DO i=1,N
1009          E2(i,j) = DCMPLX( -FJAC(i,j), ZERO )
1010        END DO
1011        E2(j,j) = E2(j,j) + CMPLX( Alpha, Beta )
1012      END DO
1013      CALL ZGETRF(N,N,E2,N,IP2,ISING)     
1014#else 
1015      DO i=1,LU_NONZERO
1016         E2(i) = DCMPLX( -FJAC(i), ZERO )
1017      END DO
1018      DO i=1,NVAR
1019         j=LU_DIAG(i); E2(j)=E2(j) + CMPLX( Alpha, Beta )
1020      END DO
1021      CALL KppDecompCmplx(E2,ISING)   
1022#endif     
1023      ISTATUS(Ndec) = ISTATUS(Ndec) + 1
1024     
1025   END SUBROUTINE RK_Decomp
1026
1027
1028!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1029   SUBROUTINE RK_Solve(N,H,E1,IP1,E2,IP2,R1,R2,R3,ISING)
1030!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1031      IMPLICIT NONE
1032      INTEGER :: N,IP1(NVAR),IP2(NVAR),ISING
1033#ifdef FULL_ALGEBRA     
1034      KPP_REAL    :: E1(NVAR,NVAR)
1035      COMPLEX(kind=dp) :: E2(NVAR,NVAR)
1036      INTEGER :: ISING
1037#else     
1038      KPP_REAL    :: E1(LU_NONZERO)
1039      COMPLEX(kind=dp) :: E2(LU_NONZERO)
1040#endif     
1041      KPP_REAL    :: R1(N),R2(N),R3(N)
1042      KPP_REAL    :: H, x1, x2, x3
1043      COMPLEX(kind=dp) :: BC(N)
1044      INTEGER :: i
1045!     
1046     ! Z <- h^{-1) T^{-1) A^{-1) x Z
1047      DO i=1,N
1048          x1 = R1(i)/H; x2 = R2(i)/H; x3 = R3(i)/H
1049          R1(i) = rkTinvAinv(1,1)*x1 + rkTinvAinv(1,2)*x2 + rkTinvAinv(1,3)*x3
1050          R2(i) = rkTinvAinv(2,1)*x1 + rkTinvAinv(2,2)*x2 + rkTinvAinv(2,3)*x3
1051          R3(i) = rkTinvAinv(3,1)*x1 + rkTinvAinv(3,2)*x2 + rkTinvAinv(3,3)*x3
1052      END DO
1053
1054#ifdef FULL_ALGEBRA     
1055      CALL DGETRS ('N',N,1,E1,N,IP1,R1,N,ISING) 
1056#else     
1057      CALL KppSolve (E1,R1)
1058#endif     
1059!     
1060      DO i=1,N
1061        BC(i) = DCMPLX(R2(i),R3(i))
1062      END DO
1063#ifdef FULL_ALGEBRA     
1064      CALL ZGETRS ('N',N,1,E2,N,IP2,BC,N,ISING) 
1065#else     
1066      CALL KppSolveCmplx (E2,BC)
1067#endif     
1068      DO i=1,N
1069        R2(i) = DBLE( BC(i) )
1070        R3(i) = AIMAG( BC(i) )
1071      END DO
1072
1073      ! Z <- T x Z
1074      DO i=1,N
1075          x1 = R1(i); x2 = R2(i); x3 = R3(i)
1076          R1(i) = rkT(1,1)*x1 + rkT(1,2)*x2 + rkT(1,3)*x3
1077          R2(i) = rkT(2,1)*x1 + rkT(2,2)*x2 + rkT(2,3)*x3
1078          R3(i) = rkT(3,1)*x1 + rkT(3,2)*x2 + rkT(3,3)*x3
1079      END DO
1080
1081      ISTATUS(Nsol) = ISTATUS(Nsol) + 1
1082
1083   END SUBROUTINE RK_Solve
1084
1085
1086!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1087   SUBROUTINE RK_ErrorEstimate(N,H,T,Y,F0,   &
1088               E1,IP1,Z1,Z2,Z3,SCAL,Err,     &
1089               FirstStep,Reject)
1090!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1091      IMPLICIT NONE
1092     
1093      INTEGER :: N
1094#ifdef FULL_ALGEBRA     
1095      KPP_REAL :: E1(NVAR,NVAR)
1096      INTEGER  :: ISING
1097#else     
1098      KPP_REAL :: E1(LU_NONZERO)
1099#endif     
1100      KPP_REAL :: SCAL(N),Z1(N),Z2(N),Z3(N),F1(N),F2(N), &
1101                        F0(N),Y(N),TMP(N),T,H
1102      INTEGER :: IP1(N), i
1103      LOGICAL FirstStep,Reject
1104      KPP_REAL :: HrkE1,HrkE2,HrkE3,Err
1105
1106      HrkE1  = rkE(1)/H
1107      HrkE2  = rkE(2)/H
1108      HrkE3  = rkE(3)/H
1109
1110      DO  i=1,N
1111         F2(i)  = HrkE1*Z1(i)+HrkE2*Z2(i)+HrkE3*Z3(i)
1112         TMP(i) = rkE(0)*F0(i) + F2(i)
1113      END DO
1114
1115
1116#ifdef FULL_ALGEBRA     
1117      CALL DGETRS ('N',N,1,E1,N,IP1,TMP,N,ISING) 
1118      IF ((rkMethod==R1A).OR.(rkMethod==GAU).OR.(rkMethod==L3A)) THEN
1119           CALL DGETRS ('N',N,1,E1,N,IP1,TMP,N,ISING)
1120      END IF       
1121      IF (rkMethod==GAU) THEN
1122           CALL DGETRS ('N',N,1,E1,N,IP1,TMP,N,ISING)
1123      END IF       
1124#else     
1125      CALL KppSolve (E1, TMP)
1126      IF ((rkMethod==R1A).OR.(rkMethod==GAU).OR.(rkMethod==L3A)) THEN
1127            CALL KppSolve (E1,TMP)
1128      END IF       
1129      IF (rkMethod==GAU) THEN
1130            CALL KppSolve (E1,TMP)
1131      END IF       
1132#endif     
1133
1134      Err = RK_ErrorNorm(N,SCAL,TMP)
1135!
1136      IF (Err < ONE) RETURN
1137firej:IF (FirstStep.OR.Reject) THEN
1138          DO i=1,N
1139             TMP(i)=Y(i)+TMP(i)
1140          END DO
1141          CALL FUN_CHEM(T,TMP,F1)   
1142          ISTATUS(Nfun) = ISTATUS(Nfun) + 1
1143          DO i=1,N
1144             TMP(i)=F1(i)+F2(i)
1145          END DO
1146
1147#ifdef FULL_ALGEBRA     
1148          CALL DGETRS ('N',N,1,E1,N,IP1,TMP,N,ISING) 
1149#else     
1150          CALL KppSolve (E1, TMP)
1151#endif     
1152          Err = RK_ErrorNorm(N,SCAL,TMP)
1153       END IF firej
1154 
1155   END SUBROUTINE RK_ErrorEstimate
1156
1157
1158!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1159   KPP_REAL FUNCTION RK_ErrorNorm(N,SCAL,DY)
1160!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1161      IMPLICIT NONE
1162     
1163      INTEGER :: N
1164      KPP_REAL :: SCAL(N),DY(N)
1165      INTEGER :: i
1166
1167      RK_ErrorNorm = ZERO
1168        DO i=1,N
1169          RK_ErrorNorm = RK_ErrorNorm + (DY(i)*SCAL(i))**2
1170        END DO
1171      RK_ErrorNorm = MAX( SQRT(RK_ErrorNorm/N), 1.0d-10 )
1172 
1173   END FUNCTION RK_ErrorNorm
1174 
1175!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1176  SUBROUTINE Radau2A_Coefficients
1177!    The coefficients of the 3-stage Radau-2A method
1178!    (given to ~30 accurate digits)
1179!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
1180      IMPLICIT NONE
1181! The coefficients of the Radau2A method
1182      KPP_REAL :: b0
1183
1184!      b0 = 1.0d0
1185      IF (SdirkError) THEN
1186        b0 = 0.2d-1
1187      ELSE
1188        b0 = 0.5d-1
1189      END IF
1190
1191! The coefficients of the Radau2A method
1192      rkMethod = R2A
1193
1194      rkA(1,1) =  1.968154772236604258683861429918299d-1
1195      rkA(1,2) = -6.55354258501983881085227825696087d-2
1196      rkA(1,3) =  2.377097434822015242040823210718965d-2
1197      rkA(2,1) =  3.944243147390872769974116714584975d-1
1198      rkA(2,2) =  2.920734116652284630205027458970589d-1
1199      rkA(2,3) = -4.154875212599793019818600988496743d-2
1200      rkA(3,1) =  3.764030627004672750500754423692808d-1
1201      rkA(3,2) =  5.124858261884216138388134465196080d-1
1202      rkA(3,3) =  1.111111111111111111111111111111111d-1
1203
1204      rkB(1) = 3.764030627004672750500754423692808d-1
1205      rkB(2) = 5.124858261884216138388134465196080d-1
1206      rkB(3) = 1.111111111111111111111111111111111d-1
1207
1208      rkC(1) = 1.550510257216821901802715925294109d-1
1209      rkC(2) = 6.449489742783178098197284074705891d-1
1210      rkC(3) = 1.0d0
1211     
1212      ! New solution: H* Sum B_j*f(Z_j) = Sum D_j*Z_j
1213      rkD(1) = 0.0d0
1214      rkD(2) = 0.0d0
1215      rkD(3) = 1.0d0
1216
1217      ! Classical error estimator:
1218      ! H* Sum (B_j-Bhat_j)*f(Z_j) = H*E(0)*f(0) + Sum E_j*Z_j
1219      rkE(0) = 1.0d0*b0
1220      rkE(1) = -10.04880939982741556246032950764708d0*b0
1221      rkE(2) = 1.382142733160748895793662840980412d0*b0
1222      rkE(3) = -.3333333333333333333333333333333333d0*b0
1223
1224      ! Sdirk error estimator
1225      rkBgam(0) = b0
1226      rkBgam(1) = .3764030627004672750500754423692807d0-1.558078204724922382431975370686279d0*b0
1227      rkBgam(2) = .8914115380582557157653087040196118d0*b0+.5124858261884216138388134465196077d0
1228      rkBgam(3) = -.1637777184845662566367174924883037d0-.3333333333333333333333333333333333d0*b0
1229      rkBgam(4) = .2748888295956773677478286035994148d0
1230
1231      ! H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j
1232      rkTheta(1) = -1.520677486405081647234271944611547d0-10.04880939982741556246032950764708d0*b0
1233      rkTheta(2) = 2.070455145596436382729929151810376d0+1.382142733160748895793662840980413d0*b0
1234      rkTheta(3) = -.3333333333333333333333333333333333d0*b0-.3744441479783868387391430179970741d0
1235
1236      ! Local order of error estimator
1237      IF (b0==0.0d0) THEN
1238        rkELO  = 6.0d0
1239      ELSE     
1240        rkELO  = 4.0d0
1241      END IF   
1242
1243      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1244      !~~~> Diagonalize the RK matrix:               
1245      ! rkTinv * inv(rkA) * rkT =         
1246      !           |  rkGamma      0           0     |
1247      !           |      0      rkAlpha   -rkBeta   |
1248      !           |      0      rkBeta     rkAlpha  |
1249      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1250
1251      rkGamma = 3.637834252744495732208418513577775d0
1252      rkAlpha = 2.681082873627752133895790743211112d0
1253      rkBeta  = 3.050430199247410569426377624787569d0
1254
1255      rkT(1,1) =  9.443876248897524148749007950641664d-2
1256      rkT(1,2) = -1.412552950209542084279903838077973d-1
1257      rkT(1,3) = -3.00291941051474244918611170890539d-2
1258      rkT(2,1) =  2.502131229653333113765090675125018d-1
1259      rkT(2,2) =  2.041293522937999319959908102983381d-1
1260      rkT(2,3) =  3.829421127572619377954382335998733d-1
1261      rkT(3,1) =  1.0d0
1262      rkT(3,2) =  1.0d0
1263      rkT(3,3) =  0.0d0
1264
1265      rkTinv(1,1) =  4.178718591551904727346462658512057d0
1266      rkTinv(1,2) =  3.27682820761062387082533272429617d-1
1267      rkTinv(1,3) =  5.233764454994495480399309159089876d-1
1268      rkTinv(2,1) = -4.178718591551904727346462658512057d0
1269      rkTinv(2,2) = -3.27682820761062387082533272429617d-1
1270      rkTinv(2,3) =  4.766235545005504519600690840910124d-1
1271      rkTinv(3,1) = -5.02872634945786875951247343139544d-1
1272      rkTinv(3,2) =  2.571926949855605429186785353601676d0
1273      rkTinv(3,3) = -5.960392048282249249688219110993024d-1
1274
1275      rkTinvAinv(1,1) =  1.520148562492775501049204957366528d+1
1276      rkTinvAinv(1,2) =  1.192055789400527921212348994770778d0
1277      rkTinvAinv(1,3) =  1.903956760517560343018332287285119d0
1278      rkTinvAinv(2,1) = -9.669512977505946748632625374449567d0
1279      rkTinvAinv(2,2) = -8.724028436822336183071773193986487d0
1280      rkTinvAinv(2,3) =  3.096043239482439656981667712714881d0
1281      rkTinvAinv(3,1) = -1.409513259499574544876303981551774d+1
1282      rkTinvAinv(3,2) =  5.895975725255405108079130152868952d0
1283      rkTinvAinv(3,3) = -1.441236197545344702389881889085515d-1
1284
1285      rkAinvT(1,1) = .3435525649691961614912493915818282d0
1286      rkAinvT(1,2) = -.4703191128473198422370558694426832d0
1287      rkAinvT(1,3) = .3503786597113668965366406634269080d0
1288      rkAinvT(2,1) = .9102338692094599309122768354288852d0
1289      rkAinvT(2,2) = 1.715425895757991796035292755937326d0
1290      rkAinvT(2,3) = .4040171993145015239277111187301784d0
1291      rkAinvT(3,1) = 3.637834252744495732208418513577775d0
1292      rkAinvT(3,2) = 2.681082873627752133895790743211112d0
1293      rkAinvT(3,3) = -3.050430199247410569426377624787569d0
1294
1295  END SUBROUTINE Radau2A_Coefficients
1296
1297   
1298
1299   
1300!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1301  SUBROUTINE Lobatto3C_Coefficients
1302!    The coefficients of the 3-stage Lobatto-3C method
1303!    (given to ~30 accurate digits)
1304!    The parameter b0 can be chosen to tune the error estimator
1305!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
1306      IMPLICIT NONE
1307      KPP_REAL :: b0
1308
1309      rkMethod = L3C
1310
1311!      b0 = 1.0d0
1312      IF (SdirkError) THEN
1313        b0 = 0.2d0
1314      ELSE
1315        b0 = 0.5d0
1316      END IF
1317! The coefficients of the Lobatto3C method
1318
1319      rkA(1,1) =  .1666666666666666666666666666666667d0
1320      rkA(1,2) = -.3333333333333333333333333333333333d0
1321      rkA(1,3) =  .1666666666666666666666666666666667d0
1322      rkA(2,1) =  .1666666666666666666666666666666667d0
1323      rkA(2,2) =  .4166666666666666666666666666666667d0
1324      rkA(2,3) = -.8333333333333333333333333333333333d-1
1325      rkA(3,1) =  .1666666666666666666666666666666667d0
1326      rkA(3,2) =  .6666666666666666666666666666666667d0
1327      rkA(3,3) =  .1666666666666666666666666666666667d0
1328
1329      rkB(1) = .1666666666666666666666666666666667d0
1330      rkB(2) = .6666666666666666666666666666666667d0
1331      rkB(3) = .1666666666666666666666666666666667d0
1332
1333      rkC(1) = 0.0d0
1334      rkC(2) = 0.5d0
1335      rkC(3) = 1.0d0
1336
1337      ! Classical error estimator, embedded solution:
1338      rkBhat(0) = b0
1339      rkBhat(1) = .16666666666666666666666666666666667d0-b0
1340      rkBhat(2) = .66666666666666666666666666666666667d0
1341      rkBhat(3) = .16666666666666666666666666666666667d0
1342     
1343      ! New solution: h Sum_j b_j f(Z_j) = sum d_j Z_j
1344      rkD(1) = 0.0d0
1345      rkD(2) = 0.0d0
1346      rkD(3) = 1.0d0
1347
1348      ! Classical error estimator:
1349      !   H* Sum (B_j-Bhat_j)*f(Z_j) = H*E(0)*f(0) + Sum E_j*Z_j
1350      rkE(0) =   .3808338772072650364017425226487022*b0
1351      rkE(1) = -1.142501631621795109205227567946107*b0
1352      rkE(2) = -1.523335508829060145606970090594809*b0
1353      rkE(3) =   .3808338772072650364017425226487022*b0
1354
1355      ! Sdirk error estimator
1356      rkBgam(0) = b0
1357      rkBgam(1) = .1666666666666666666666666666666667d0-1.d0*b0
1358      rkBgam(2) = .6666666666666666666666666666666667d0
1359      rkBgam(3) = -.2141672105405983697350758559820354d0
1360      rkBgam(4) = .3808338772072650364017425226487021d0
1361
1362      ! H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j
1363      rkTheta(1) = -3.d0*b0-.3808338772072650364017425226487021d0
1364      rkTheta(2) = -4.d0*b0+1.523335508829060145606970090594808d0
1365      rkTheta(3) = -.142501631621795109205227567946106d0+b0
1366
1367      ! Local order of error estimator
1368      IF (b0==0.0d0) THEN
1369        rkELO  = 5.0d0
1370      ELSE     
1371        rkELO  = 4.0d0
1372      END IF   
1373
1374      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1375      !~~~> Diagonalize the RK matrix:               
1376      ! rkTinv * inv(rkA) * rkT =         
1377      !           |  rkGamma      0           0     |
1378      !           |      0      rkAlpha   -rkBeta   |
1379      !           |      0      rkBeta     rkAlpha  |
1380      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1381
1382      rkGamma = 2.625816818958466716011888933765284d0
1383      rkAlpha = 1.687091590520766641994055533117359d0
1384      rkBeta  = 2.508731754924880510838743672432351d0
1385
1386      rkT(1,1) = 1.d0
1387      rkT(1,2) = 1.d0
1388      rkT(1,3) = 0.d0
1389      rkT(2,1) = .4554100411010284672111720348287483d0
1390      rkT(2,2) = -.6027050205505142336055860174143743d0
1391      rkT(2,3) = -.4309321229203225731070721341350346d0
1392      rkT(3,1) = 2.195823345445647152832799205549709d0
1393      rkT(3,2) = -1.097911672722823576416399602774855d0
1394      rkT(3,3) = .7850032632435902184104551358922130d0
1395
1396      rkTinv(1,1) = .4205559181381766909344950150991349d0
1397      rkTinv(1,2) = .3488903392193734304046467270632057d0
1398      rkTinv(1,3) = .1915253879645878102698098373933487d0
1399      rkTinv(2,1) = .5794440818618233090655049849008650d0
1400      rkTinv(2,2) = -.3488903392193734304046467270632057d0
1401      rkTinv(2,3) = -.1915253879645878102698098373933487d0
1402      rkTinv(3,1) = -.3659705575742745254721332009249516d0
1403      rkTinv(3,2) = -1.463882230297098101888532803699806d0
1404      rkTinv(3,3) = .4702733607340189781407813565524989d0
1405
1406      rkTinvAinv(1,1) = 1.104302803159744452668648155627548d0
1407      rkTinvAinv(1,2) = .916122120694355522658740710823143d0
1408      rkTinvAinv(1,3) = .5029105849749601702795812241441172d0
1409      rkTinvAinv(2,1) = 1.895697196840255547331351844372453d0
1410      rkTinvAinv(2,2) = 3.083877879305644477341259289176857d0
1411      rkTinvAinv(2,3) = -1.502910584974960170279581224144117d0
1412      rkTinvAinv(3,1) = .8362439183082935036129145574774502d0
1413      rkTinvAinv(3,2) = -3.344975673233174014451658229909802d0
1414      rkTinvAinv(3,3) = .312908409479233358005944466882642d0
1415
1416      rkAinvT(1,1) = 2.625816818958466716011888933765282d0
1417      rkAinvT(1,2) = 1.687091590520766641994055533117358d0
1418      rkAinvT(1,3) = -2.508731754924880510838743672432351d0
1419      rkAinvT(2,1) = 1.195823345445647152832799205549710d0
1420      rkAinvT(2,2) = -2.097911672722823576416399602774855d0
1421      rkAinvT(2,3) = .7850032632435902184104551358922130d0
1422      rkAinvT(3,1) = 5.765829871932827589653709477334136d0
1423      rkAinvT(3,2) = .1170850640335862051731452613329320d0
1424      rkAinvT(3,3) = 4.078738281412060947659653944216779d0
1425
1426  END SUBROUTINE Lobatto3C_Coefficients
1427
1428!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1429  SUBROUTINE Gauss_Coefficients
1430!    The coefficients of the 3-stage Gauss method
1431!    (given to ~30 accurate digits)
1432!    The parameter b3 can be chosen by the user
1433!    to tune the error estimator
1434!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
1435      IMPLICIT NONE
1436      KPP_REAL :: b0
1437! The coefficients of the Gauss method
1438
1439
1440      rkMethod = GAU
1441     
1442!      b0 = 4.0d0
1443      b0 = 0.1d0
1444     
1445! The coefficients of the Gauss method
1446
1447      rkA(1,1) =  .1388888888888888888888888888888889d0
1448      rkA(1,2) = -.359766675249389034563954710966045d-1
1449      rkA(1,3) =  .97894440153083260495800422294756d-2
1450      rkA(2,1) =  .3002631949808645924380249472131556d0
1451      rkA(2,2) =  .2222222222222222222222222222222222d0
1452      rkA(2,3) = -.224854172030868146602471694353778d-1
1453      rkA(3,1) =  .2679883337624694517281977355483022d0
1454      rkA(3,2) =  .4804211119693833479008399155410489d0
1455      rkA(3,3) =  .1388888888888888888888888888888889d0
1456
1457      rkB(1) = .2777777777777777777777777777777778d0
1458      rkB(2) = .4444444444444444444444444444444444d0
1459      rkB(3) = .2777777777777777777777777777777778d0
1460
1461      rkC(1) = .1127016653792583114820734600217600d0
1462      rkC(2) = .5000000000000000000000000000000000d0
1463      rkC(3) = .8872983346207416885179265399782400d0
1464
1465      ! Classical error estimator, embedded solution:
1466      rkBhat(0) = b0
1467      rkBhat(1) =-1.4788305577012361475298775666303999d0*b0 &
1468                  +.27777777777777777777777777777777778d0
1469      rkBhat(2) =  .44444444444444444444444444444444444d0 &
1470                  +.66666666666666666666666666666666667d0*b0
1471      rkBhat(3) = -.18783610896543051913678910003626672d0*b0 &
1472                  +.27777777777777777777777777777777778d0
1473
1474      ! New solution: h Sum_j b_j f(Z_j) = sum d_j Z_j
1475      rkD(1) = .1666666666666666666666666666666667d1
1476      rkD(2) = -.1333333333333333333333333333333333d1
1477      rkD(3) = .1666666666666666666666666666666667d1
1478
1479      ! Classical error estimator:
1480      !   H* Sum (B_j-Bhat_j)*f(Z_j) = H*E(0)*f(0) + Sum E_j*Z_j
1481      rkE(0) = .2153144231161121782447335303806954d0*b0
1482      rkE(1) = -2.825278112319014084275808340593191d0*b0
1483      rkE(2) = .2870858974881495709929780405075939d0*b0
1484      rkE(3) = -.4558086256248162565397206448274867d-1*b0
1485
1486      ! Sdirk error estimator
1487      rkBgam(0) = 0.d0
1488      rkBgam(1) = .2373339543355109188382583162660537d0
1489      rkBgam(2) = .5879873931885192299409334646982414d0
1490      rkBgam(3) = -.4063577064014232702392531134499046d-1
1491      rkBgam(4) = .2153144231161121782447335303806955d0
1492
1493      ! H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j
1494      rkTheta(1) = -2.594040933093095272574031876464493d0
1495      rkTheta(2) = 1.824611539036311947589425112250199d0
1496      rkTheta(3) = .1856563166634371860478043996459493d0
1497
1498      ! ELO = local order of classical error estimator
1499      rkELO = 4.0d0
1500
1501      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1502      !~~~> Diagonalize the RK matrix:               
1503      ! rkTinv * inv(rkA) * rkT =         
1504      !           |  rkGamma      0           0     |
1505      !           |      0      rkAlpha   -rkBeta   |
1506      !           |      0      rkBeta     rkAlpha  |
1507      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1508
1509      rkGamma = 4.644370709252171185822941421408064d0
1510      rkAlpha = 3.677814645373914407088529289295970d0
1511      rkBeta  = 3.508761919567443321903661209182446d0
1512
1513      rkT(1,1) =  .7215185205520017032081769924397664d-1
1514      rkT(1,2) = -.8224123057363067064866206597516454d-1
1515      rkT(1,3) = -.6012073861930850173085948921439054d-1
1516      rkT(2,1) =  .1188325787412778070708888193730294d0
1517      rkT(2,2) =  .5306509074206139504614411373957448d-1
1518      rkT(2,3) =  .3162050511322915732224862926182701d0
1519      rkT(3,1) = 1.d0
1520      rkT(3,2) = 1.d0
1521      rkT(3,3) = 0.d0
1522
1523      rkTinv(1,1) =  5.991698084937800775649580743981285d0
1524      rkTinv(1,2) =  1.139214295155735444567002236934009d0
1525      rkTinv(1,3) =   .4323121137838583855696375901180497d0
1526      rkTinv(2,1) = -5.991698084937800775649580743981285d0
1527      rkTinv(2,2) = -1.139214295155735444567002236934009d0
1528      rkTinv(2,3) =   .5676878862161416144303624098819503d0
1529      rkTinv(3,1) = -1.246213273586231410815571640493082d0
1530      rkTinv(3,2) =  2.925559646192313662599230367054972d0
1531      rkTinv(3,3) =  -.2577352012734324923468722836888244d0
1532
1533      rkTinvAinv(1,1) =  27.82766708436744962047620566703329d0
1534      rkTinvAinv(1,2) =   5.290933503982655311815946575100597d0
1535      rkTinvAinv(1,3) =   2.007817718512643701322151051660114d0
1536      rkTinvAinv(2,1) = -17.66368928942422710690385180065675d0
1537      rkTinvAinv(2,2) = -14.45491129892587782538830044147713d0
1538      rkTinvAinv(2,3) =   2.992182281487356298677848948339886d0
1539      rkTinvAinv(3,1) = -25.60678350282974256072419392007303d0
1540      rkTinvAinv(3,2) =   6.762434375611708328910623303779923d0
1541      rkTinvAinv(3,3) =   1.043979339483109825041215970036771d0
1542     
1543      rkAinvT(1,1) = .3350999483034677402618981153470483d0
1544      rkAinvT(1,2) = -.5134173605009692329246186488441294d0
1545      rkAinvT(1,3) = .6745196507033116204327635673208923d-1
1546      rkAinvT(2,1) = .5519025480108928886873752035738885d0
1547      rkAinvT(2,2) = 1.304651810077110066076640761092008d0
1548      rkAinvT(2,3) = .9767507983414134987545585703726984d0
1549      rkAinvT(3,1) = 4.644370709252171185822941421408064d0
1550      rkAinvT(3,2) = 3.677814645373914407088529289295970d0
1551      rkAinvT(3,3) = -3.508761919567443321903661209182446d0
1552     
1553  END SUBROUTINE Gauss_Coefficients
1554
1555
1556
1557!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1558  SUBROUTINE Radau1A_Coefficients
1559!    The coefficients of the 3-stage Gauss method
1560!    (given to ~30 accurate digits)
1561!    The parameter b3 can be chosen by the user
1562!    to tune the error estimator
1563!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
1564      IMPLICIT NONE
1565!      KPP_REAL :: b0 = 0.3d0
1566      KPP_REAL :: b0 = 0.1d0
1567
1568! The coefficients of the Radau1A method
1569
1570      rkMethod = R1A
1571
1572      rkA(1,1) =  .1111111111111111111111111111111111d0
1573      rkA(1,2) = -.1916383190435098943442935597058829d0
1574      rkA(1,3) =  .8052720793239878323318244859477174d-1
1575      rkA(2,1) =  .1111111111111111111111111111111111d0
1576      rkA(2,2) =  .2920734116652284630205027458970589d0
1577      rkA(2,3) = -.481334970546573839513422644787591d-1
1578      rkA(3,1) =  .1111111111111111111111111111111111d0
1579      rkA(3,2) =  .5370223859435462728402311533676479d0
1580      rkA(3,3) =  .1968154772236604258683861429918299d0
1581
1582      rkB(1) = .1111111111111111111111111111111111d0
1583      rkB(2) = .5124858261884216138388134465196080d0
1584      rkB(3) = .3764030627004672750500754423692808d0
1585
1586      rkC(1) = 0.d0
1587      rkC(2) = .3550510257216821901802715925294109d0
1588      rkC(3) = .8449489742783178098197284074705891d0
1589
1590      ! Classical error estimator, embedded solution:
1591      rkBhat(0) = b0
1592      rkBhat(1) = .11111111111111111111111111111111111d0-b0
1593      rkBhat(2) = .51248582618842161383881344651960810d0
1594      rkBhat(3) = .37640306270046727505007544236928079d0
1595
1596      ! New solution: H* Sum B_j*f(Z_j) = Sum D_j*Z_j
1597      rkD(1) = .3333333333333333333333333333333333d0
1598      rkD(2) = -.8914115380582557157653087040196127d0
1599      rkD(3) = .1558078204724922382431975370686279d1
1600
1601      ! Classical error estimator:
1602      ! H* Sum (b_j-bhat_j) f(Z_j) = H*E(0)*F(0) + Sum E_j Z_j
1603      rkE(0) =   .2748888295956773677478286035994148d0*b0
1604      rkE(1) = -1.374444147978386838739143017997074d0*b0
1605      rkE(2) = -1.335337922441686804550326197041126d0*b0
1606      rkE(3) =   .235782604058977333559011782643466d0*b0
1607
1608      ! Sdirk error estimator
1609      rkBgam(0) = 0.0d0
1610      rkBgam(1) = .1948150124588532186183490991130616d-1
1611      rkBgam(2) = .7575249005733381398986810981093584d0
1612      rkBgam(3) = -.518952314149008295083446116200793d-1
1613      rkBgam(4) = .2748888295956773677478286035994148d0
1614
1615      ! H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j
1616      rkTheta(1) = -1.224370034375505083904362087063351d0
1617      rkTheta(2) = .9340045331532641409047527962010133d0
1618      rkTheta(3) = .4656990124352088397561234800640929d0
1619
1620      ! ELO = local order of classical error estimator
1621      rkELO = 4.0d0
1622
1623      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1624      !~~~> Diagonalize the RK matrix:               
1625      ! rkTinv * inv(rkA) * rkT =         
1626      !           |  rkGamma      0           0     |
1627      !           |      0      rkAlpha   -rkBeta   |
1628      !           |      0      rkBeta     rkAlpha  |
1629      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1630
1631      rkGamma = 3.637834252744495732208418513577775d0
1632      rkAlpha = 2.681082873627752133895790743211112d0
1633      rkBeta  = 3.050430199247410569426377624787569d0
1634
1635      rkT(1,1) =  .424293819848497965354371036408369d0
1636      rkT(1,2) = -.3235571519651980681202894497035503d0
1637      rkT(1,3) = -.522137786846287839586599927945048d0
1638      rkT(2,1) =  .57594609499806128896291585429339d-1
1639      rkT(2,2) =  .3148663231849760131614374283783d-2
1640      rkT(2,3) =  .452429247674359778577728510381731d0
1641      rkT(3,1) = 1.d0
1642      rkT(3,2) = 1.d0
1643      rkT(3,3) = 0.d0
1644
1645      rkTinv(1,1) = 1.233523612685027760114769983066164d0
1646      rkTinv(1,2) = 1.423580134265707095505388133369554d0
1647      rkTinv(1,3) = .3946330125758354736049045150429624d0
1648      rkTinv(2,1) = -1.233523612685027760114769983066164d0
1649      rkTinv(2,2) = -1.423580134265707095505388133369554d0
1650      rkTinv(2,3) = .6053669874241645263950954849570376d0
1651      rkTinv(3,1) = -.1484438963257383124456490049673414d0
1652      rkTinv(3,2) = 2.038974794939896109682070471785315d0
1653      rkTinv(3,3) = -.544501292892686735299355831692542d-1
1654
1655      rkTinvAinv(1,1) =  4.487354449794728738538663081025420d0
1656      rkTinvAinv(1,2) =  5.178748573958397475446442544234494d0
1657      rkTinvAinv(1,3) =  1.435609490412123627047824222335563d0
1658      rkTinvAinv(2,1) = -2.854361287939276673073807031221493d0
1659      rkTinvAinv(2,2) = -1.003648660720543859000994063139137d+1
1660      rkTinvAinv(2,3) =  1.789135380979465422050817815017383d0
1661      rkTinvAinv(3,1) = -4.160768067752685525282947313530352d0
1662      rkTinvAinv(3,2) =  1.124128569859216916690209918405860d0
1663      rkTinvAinv(3,3) =  1.700644430961823796581896350418417d0
1664
1665      rkAinvT(1,1) = 1.543510591072668287198054583233180d0
1666      rkAinvT(1,2) = -2.460228411937788329157493833295004d0
1667      rkAinvT(1,3) = -.412906170450356277003910443520499d0
1668      rkAinvT(2,1) = .209519643211838264029272585946993d0
1669      rkAinvT(2,2) = 1.388545667194387164417459732995766d0
1670      rkAinvT(2,3) = 1.20339553005832004974976023130002d0
1671      rkAinvT(3,1) = 3.637834252744495732208418513577775d0
1672      rkAinvT(3,2) = 2.681082873627752133895790743211112d0
1673      rkAinvT(3,3) = -3.050430199247410569426377624787569d0
1674
1675  END SUBROUTINE Radau1A_Coefficients
1676
1677
1678!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1679  SUBROUTINE Lobatto3A_Coefficients
1680!    The coefficients of the 4-stage Lobatto-3A method
1681!    (given to ~30 accurate digits)
1682!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~     
1683      IMPLICIT NONE
1684
1685! The coefficients of the Lobatto-3A method
1686
1687      rkMethod = L3A
1688
1689      rkA(0,0) = 0.0d0
1690      rkA(0,1) = 0.0d0
1691      rkA(0,2) = 0.0d0
1692      rkA(0,3) = 0.0d0
1693      rkA(1,0) = .11030056647916491413674311390609397d0
1694      rkA(1,1) = .1896994335208350858632568860939060d0
1695      rkA(1,2) = -.339073642291438837776604807792215d-1
1696      rkA(1,3) = .1030056647916491413674311390609397d-1
1697      rkA(2,0) = .73032766854168419196590219427239365d-1
1698      rkA(2,1) = .4505740308958105504443271474458881d0
1699      rkA(2,2) = .2269672331458315808034097805727606d0
1700      rkA(2,3) = -.2696723314583158080340978057276063d-1
1701      rkA(3,0) = .83333333333333333333333333333333333d-1
1702      rkA(3,1) = .4166666666666666666666666666666667d0
1703      rkA(3,2) = .4166666666666666666666666666666667d0
1704      rkA(3,3) = .8333333333333333333333333333333333d-1
1705     
1706      rkB(0) = .83333333333333333333333333333333333d-1
1707      rkB(1) = .4166666666666666666666666666666667d0
1708      rkB(2) = .4166666666666666666666666666666667d0
1709      rkB(3) = .8333333333333333333333333333333333d-1
1710
1711      rkC(0) = 0.0d0
1712      rkC(1) = .2763932022500210303590826331268724d0
1713      rkC(2) = .7236067977499789696409173668731276d0
1714      rkC(3) = 1.0d0
1715
1716      ! New solution: H*Sum B_j*f(Z_j) = Sum D_j*Z_j
1717      rkD(0) = 0.0d0
1718      rkD(1) = 0.0d0
1719      rkD(2) = 0.0d0
1720      rkD(3) = 1.0d0
1721     
1722      ! Classical error estimator, embedded solution:
1723      rkBhat(0) = .90909090909090909090909090909090909d-1
1724      rkBhat(1) = .39972675774621371442114262372173276d0
1725      rkBhat(2) = .43360657558711961891219070961160058d0
1726      rkBhat(3) = .15151515151515151515151515151515152d-1
1727
1728      ! Classical error estimator:
1729      ! H* Sum (B_j-Bhat_j)*f(Z_j) = H*E(0)*f(0) + Sum E_j*Z_j
1730
1731      rkE(0) =  .1957403846510110711315759367097231d-1
1732      rkE(1) = -.1986820345632580910316020806676438d0
1733      rkE(2) =  .1660586371214229125096727578826900d0
1734      rkE(3) = -.9787019232550553556578796835486154d-1
1735
1736      ! Sdirk error estimator:
1737      rkF(0) =  0.0d0
1738      rkF(1) = -.66535815876916686607437314126436349d0
1739      rkF(2) =  1.7419302743497277572980407931678409d0
1740      rkF(3) = -1.2918865386966730694684011822841728d0
1741       
1742      ! ELO = local order of classical error estimator
1743      rkELO = 4.0d0
1744     
1745      ! Sdirk error estimator:
1746      rkBgam(0) =  .2950472755430528877214995073815946d-1
1747      rkBgam(1) =  .5370310883226113978352873633882769d0
1748      rkBgam(2) =  .2963022450107219354980459699450564d0
1749      rkBgam(3) = -.7815248400375080035021681445218837d-1
1750      rkBgam(4) =  .2153144231161121782447335303806956d0
1751     
1752      ! H* Sum Bgam_j*f(Z_j) = H*Bgam(0)*f(0) + Sum Theta_j*Z_j
1753      rkTheta(0) = 0.0d0
1754      rkTheta(1) = -.6653581587691668660743731412643631d0
1755      rkTheta(2) = 1.741930274349727757298040793167842d0
1756      rkTheta(3) = -.291886538696673069468401182284174d0
1757
1758
1759      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1760      !~~~> Diagonalize the RK matrix:               
1761      ! rkTinv * inv(rkA) * rkT =         
1762      !           |  rkGamma      0           0     |
1763      !           |      0      rkAlpha   -rkBeta   |
1764      !           |      0      rkBeta     rkAlpha  |
1765      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1766
1767      rkGamma = 4.644370709252171185822941421408063d0
1768      rkAlpha = 3.677814645373914407088529289295968d0
1769      rkBeta  = 3.508761919567443321903661209182446d0
1770
1771      rkT(1,1) = .5303036326129938105898786144870856d-1
1772      rkT(1,2) = -.7776129960563076320631956091016914d-1
1773      rkT(1,3) = .6043307469475508514468017399717112d-2
1774      rkT(2,1) = .2637242522173698467283726114649606d0
1775      rkT(2,2) = .2193839918662961493126393244533346d0
1776      rkT(2,3) = .3198765142300936188514264752235344d0
1777      rkT(3,1) = 1.d0
1778      rkT(3,2) = 1.d0
1779      rkT(3,3) = 0.d0
1780
1781      rkTinv(1,1) = 7.695032983257654470769069079238553d0
1782      rkTinv(1,2) = -.1453793830957233720334601186354032d0
1783      rkTinv(1,3) = .6302696746849084900422461036874826d0
1784      rkTinv(2,1) = -7.695032983257654470769069079238553d0
1785      rkTinv(2,2) = .1453793830957233720334601186354032d0
1786      rkTinv(2,3) = .3697303253150915099577538963125174d0
1787      rkTinv(3,1) = -1.066660885401270392058552736086173d0
1788      rkTinv(3,2) = 3.146358406832537460764521760668932d0
1789      rkTinv(3,3) = -.7732056038202974770406168510664738d0
1790
1791      rkTinvAinv(1,1) = 35.73858579417120341641749040405149d0
1792      rkTinvAinv(1,2) = -.675195748578927863668368190236025d0
1793      rkTinvAinv(1,3) = 2.927206016036483646751158874041632d0
1794      rkTinvAinv(2,1) = -24.55824590667225493437162206039511d0
1795      rkTinvAinv(2,2) = -10.50514413892002061837750015342036
1796      rkTinvAinv(2,3) = 4.072793983963516353248841125958369d0
1797      rkTinvAinv(3,1) = -30.92301972744621647251975054630589d0
1798      rkTinvAinv(3,2) = 12.08182467154052413351908559269928d0
1799      rkTinvAinv(3,3) = -1.546411207640594954081233702132946d0
1800
1801      rkAinvT(1,1) = .2462926658317812882584158369803835d0
1802      rkAinvT(1,2) = -.2647871194157644619747121197289574d0
1803      rkAinvT(1,3) = .2950720515900466654896406799284586d0
1804      rkAinvT(2,1) = 1.224833192317784474576995878738004d0
1805      rkAinvT(2,2) = 1.929224190340981580557006261869763d0
1806      rkAinvT(2,3) = .4066803323234419988910915619080306d0
1807      rkAinvT(3,1) = 4.644370709252171185822941421408064d0
1808      rkAinvT(3,2) = 3.677814645373914407088529289295968d0
1809      rkAinvT(3,3) = -3.508761919567443321903661209182446d0
1810
1811  END SUBROUTINE Lobatto3A_Coefficients
1812 
1813!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1814  END SUBROUTINE RungeKutta ! and all its internal procedures
1815!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1816
1817!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1818  SUBROUTINE FUN_CHEM(T, V, FCT)
1819!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1820
1821    USE KPP_ROOT_Parameters
1822    USE KPP_ROOT_Global
1823    USE KPP_ROOT_Function, ONLY: Fun
1824    USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO
1825
1826    IMPLICIT NONE
1827
1828    KPP_REAL :: V(NVAR), FCT(NVAR)
1829    KPP_REAL :: T, Told
1830
1831    Told = TIME
1832    TIME = T
1833    CALL Update_SUN()
1834    CALL Update_RCONST()
1835    CALL Update_PHOTO()
1836    TIME = Told
1837   
1838    CALL Fun(V, FIX, RCONST, FCT)
1839   
1840  END SUBROUTINE FUN_CHEM
1841
1842
1843!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1844  SUBROUTINE JAC_CHEM (T, V, JF)
1845!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1846
1847    USE KPP_ROOT_Parameters
1848    USE KPP_ROOT_Global
1849    USE KPP_ROOT_JacobianSP
1850    USE KPP_ROOT_Jacobian, ONLY: Jac_SP
1851    USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO
1852
1853    IMPLICIT NONE
1854
1855    KPP_REAL :: V(NVAR), T , Told
1856#ifdef FULL_ALGEBRA   
1857    KPP_REAL :: JV(LU_NONZERO), JF(NVAR,NVAR)
1858    INTEGER :: i, j 
1859#else
1860    KPP_REAL :: JF(LU_NONZERO)
1861#endif   
1862
1863    Told = TIME
1864    TIME = T
1865    CALL Update_SUN()
1866    CALL Update_RCONST()
1867    CALL Update_PHOTO()
1868    TIME = Told
1869   
1870#ifdef FULL_ALGEBRA   
1871    CALL Jac_SP(V, FIX, RCONST, JV)
1872    DO j=1,NVAR
1873      DO i=1,NVAR
1874         JF(i,j) = 0.0d0
1875      END DO
1876    END DO
1877    DO i=1,LU_NONZERO
1878       JF(LU_IROW(i),LU_ICOL(i)) = JV(i)
1879    END DO
1880#else
1881    CALL Jac_SP(V, FIX, RCONST, JF) 
1882#endif   
1883
1884  END SUBROUTINE JAC_CHEM
1885
1886
1887!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1888
1889END MODULE KPP_ROOT_Integrator
1890
1891
Note: See TracBrowser for help on using the repository browser.