source: palm/trunk/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm/chem_gasphase_mod_Integrator.f90 @ 2696

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

Merge of branch palm4u into trunk

File size: 45.8 KB
Line 
1! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2!
3! Numerical Integrator (Time-Stepping) File
4!
5! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
6!       (http://www.cs.vt.edu/~asandu/Software/KPP)
7! KPP is distributed under GPL, the general public licence
8!       (http://www.gnu.org/copyleft/gpl.html)
9! (C) 1995-1997, V. Damian & A. Sandu, CGRER, Univ. Iowa
10! (C) 1997-2005, A. Sandu, Michigan Tech, Virginia Tech
11!     With important contributions from:
12!        M. Damian, Villanova University, USA
13!        R. Sander, Max-Planck Institute for Chemistry, Mainz, Germany
14!
15! File                 : chem_gasphase_mod_Integrator.f90
16! Time                 : Fri Dec  1 18:10:53 2017
17! Working directory    : /data/kanani/branches/palm4u/GASPHASE_PREPROC/tmp_kpp4palm
18! Equation file        : chem_gasphase_mod.kpp
19! Output root filename : chem_gasphase_mod
20!
21! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
22
23
24
25! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
26!
27! INTEGRATE - Integrator routine
28!   Arguments :
29!      TIN       - Start Time for Integration
30!      TOUT      - End Time for Integration
31!
32! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
33
34!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
35!  Rosenbrock - Implementation of several Rosenbrock methods:             !
36!               * Ros2                                                    !
37!               * Ros3                                                    !
38!               * Ros4                                                    !
39!               * Rodas3                                                  !
40!               * Rodas4                                                  !
41!  By default the code employs the KPP sparse linear algebra routines     !
42!  Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK)        !
43!                                                                         !
44!    (C)  Adrian Sandu, August 2004                                       !
45!    Virginia Polytechnic Institute and State University                  !
46!    Contact: sandu@cs.vt.edu                                             !
47!    Revised by Philipp Miehe and Adrian Sandu, May 2006                  !                               !
48!    This implementation is part of KPP - the Kinetic PreProcessor        !
49!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
50
51MODULE chem_gasphase_mod_Integrator
52
53  USE chem_gasphase_mod_Parameters, ONLY: NVAR, NFIX, NSPEC, LU_NONZERO
54  USE chem_gasphase_mod_Global
55  IMPLICIT NONE
56  PUBLIC
57  SAVE
58 
59!~~~>  Statistics on the work performed by the Rosenbrock method
60  INTEGER, PARAMETER :: Nfun=1, Njac=2, Nstp=3, Nacc=4, &
61                        Nrej=5, Ndec=6, Nsol=7, Nsng=8, &
62                        Ntexit=1, Nhexit=2, Nhnew = 3
63
64CONTAINS
65
66SUBROUTINE INTEGRATE( TIN, TOUT, &
67  ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U, IERR_U )
68
69   IMPLICIT NONE
70
71   REAL(kind=dp), INTENT(IN) :: TIN  ! Start Time
72   REAL(kind=dp), INTENT(IN) :: TOUT ! End Time
73   ! Optional input parameters and statistics
74   INTEGER,       INTENT(IN),  OPTIONAL :: ICNTRL_U(20)
75   REAL(kind=dp), INTENT(IN),  OPTIONAL :: RCNTRL_U(20)
76   INTEGER,       INTENT(OUT), OPTIONAL :: ISTATUS_U(20)
77   REAL(kind=dp), INTENT(OUT), OPTIONAL :: RSTATUS_U(20)
78   INTEGER,       INTENT(OUT), OPTIONAL :: IERR_U
79
80   REAL(kind=dp) :: RCNTRL(20), RSTATUS(20)
81   INTEGER       :: ICNTRL(20), ISTATUS(20), IERR
82
83   INTEGER, SAVE :: Ntotal = 0
84
85   ICNTRL(:)  = 0
86   RCNTRL(:)  = 0.0_dp
87   ISTATUS(:) = 0
88   RSTATUS(:) = 0.0_dp
89
90    !~~~> fine-tune the integrator:
91   ICNTRL(1) = 0        ! 0 - non-autonomous, 1 - autonomous
92   ICNTRL(2) = 0        ! 0 - vector tolerances, 1 - scalars
93
94   ! If optional parameters are given, and if they are >0,
95   ! then they overwrite default settings.
96   IF (PRESENT(ICNTRL_U)) THEN
97     WHERE(ICNTRL_U(:) > 0) ICNTRL(:) = ICNTRL_U(:)
98   END IF
99   IF (PRESENT(RCNTRL_U)) THEN
100     WHERE(RCNTRL_U(:) > 0) RCNTRL(:) = RCNTRL_U(:)
101   END IF
102
103
104   CALL Rosenbrock(NVAR,VAR,TIN,TOUT,   &
105         ATOL,RTOL,                &
106         RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR)
107
108   !~~~> Debug option: show no of steps
109   ! Ntotal = Ntotal + ISTATUS(Nstp)
110   ! PRINT*,'NSTEPS=',ISTATUS(Nstp),' (',Ntotal,')','  O3=', VAR(ind_O3)
111
112   STEPMIN = RSTATUS(Nhexit)
113   ! if optional parameters are given for output they
114   ! are updated with the return information
115   IF (PRESENT(ISTATUS_U)) ISTATUS_U(:) = ISTATUS(:)
116   IF (PRESENT(RSTATUS_U)) RSTATUS_U(:) = RSTATUS(:)
117   IF (PRESENT(IERR_U))    IERR_U       = IERR
118
119END SUBROUTINE INTEGRATE
120
121!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
122SUBROUTINE Rosenbrock(N,Y,Tstart,Tend, &
123           AbsTol,RelTol,              &
124           RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR)
125!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
126!
127!    Solves the system y'=F(t,y) using a Rosenbrock method defined by:
128!
129!     G = 1/(H*gamma(1)) - Jac(t0,Y0)
130!     T_i = t0 + Alpha(i)*H
131!     Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
132!     G * K_i = Fun( T_i, Y_i ) + \sum_{j=1}^S C(i,j)/H * K_j +
133!         gamma(i)*dF/dT(t0, Y0)
134!     Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
135!
136!    For details on Rosenbrock methods and their implementation consult:
137!      E. Hairer and G. Wanner
138!      "Solving ODEs II. Stiff and differential-algebraic problems".
139!      Springer series in computational mathematics, Springer-Verlag, 1996.
140!    The codes contained in the book inspired this implementation.
141!
142!    (C)  Adrian Sandu, August 2004
143!    Virginia Polytechnic Institute and State University
144!    Contact: sandu@cs.vt.edu
145!    Revised by Philipp Miehe and Adrian Sandu, May 2006                 
146!    This implementation is part of KPP - the Kinetic PreProcessor
147!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
148!
149!~~~>   INPUT ARGUMENTS:
150!
151!-     Y(N)    = vector of initial conditions (at T=Tstart)
152!-    [Tstart,Tend]  = time range of integration
153!     (if Tstart>Tend the integration is performed backwards in time)
154!-    RelTol, AbsTol = user precribed accuracy
155!- SUBROUTINE  Fun( T, Y, Ydot ) = ODE function,
156!                       returns Ydot = Y' = F(T,Y)
157!- SUBROUTINE  Jac( T, Y, Jcb ) = Jacobian of the ODE function,
158!                       returns Jcb = dFun/dY
159!-    ICNTRL(1:20)    = integer inputs parameters
160!-    RCNTRL(1:20)    = real inputs parameters
161!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
162!
163!~~~>     OUTPUT ARGUMENTS:
164!
165!-    Y(N)    -> vector of final states (at T->Tend)
166!-    ISTATUS(1:20)   -> integer output parameters
167!-    RSTATUS(1:20)   -> real output parameters
168!-    IERR            -> job status upon return
169!                        success (positive value) or
170!                        failure (negative value)
171!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
172!
173!~~~>     INPUT PARAMETERS:
174!
175!    Note: For input parameters equal to zero the default values of the
176!       corresponding variables are used.
177!
178!    ICNTRL(1) = 1: F = F(y)   Independent of T (AUTONOMOUS)
179!              = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
180!
181!    ICNTRL(2) = 0: AbsTol, RelTol are N-dimensional vectors
182!              = 1: AbsTol, RelTol are scalars
183!
184!    ICNTRL(3)  -> selection of a particular Rosenbrock method
185!        = 0 :    Rodas3 (default)
186!        = 1 :    Ros2
187!        = 2 :    Ros3
188!        = 3 :    Ros4
189!        = 4 :    Rodas3
190!        = 5 :    Rodas4
191!
192!    ICNTRL(4)  -> maximum number of integration steps
193!        For ICNTRL(4)=0) the default value of 100000 is used
194!
195!    RCNTRL(1)  -> Hmin, lower bound for the integration step size
196!          It is strongly recommended to keep Hmin = ZERO
197!    RCNTRL(2)  -> Hmax, upper bound for the integration step size
198!    RCNTRL(3)  -> Hstart, starting value for the integration step size
199!
200!    RCNTRL(4)  -> FacMin, lower bound on step decrease factor (default=0.2)
201!    RCNTRL(5)  -> FacMax, upper bound on step increase factor (default=6)
202!    RCNTRL(6)  -> FacRej, step decrease factor after multiple rejections
203!                          (default=0.1)
204!    RCNTRL(7)  -> FacSafe, by which the new step is slightly smaller
205!         than the predicted value  (default=0.9)
206!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
207!
208!
209!    OUTPUT ARGUMENTS:
210!    -----------------
211!
212!    T           -> T value for which the solution has been computed
213!                     (after successful return T=Tend).
214!
215!    Y(N)        -> Numerical solution at T
216!
217!    IDID        -> Reports on successfulness upon return:
218!                    = 1 for success
219!                    < 0 for error (value equals error code)
220!
221!    ISTATUS(1)  -> No. of function calls
222!    ISTATUS(2)  -> No. of jacobian calls
223!    ISTATUS(3)  -> No. of steps
224!    ISTATUS(4)  -> No. of accepted steps
225!    ISTATUS(5)  -> No. of rejected steps (except at very beginning)
226!    ISTATUS(6)  -> No. of LU decompositions
227!    ISTATUS(7)  -> No. of forward/backward substitutions
228!    ISTATUS(8)  -> No. of singular matrix decompositions
229!
230!    RSTATUS(1)  -> Texit, the time corresponding to the
231!                     computed Y upon return
232!    RSTATUS(2)  -> Hexit, last accepted step before exit
233!    RSTATUS(3)  -> Hnew, last predicted step (not yet taken)
234!                   For multiple restarts, use Hnew as Hstart
235!                     in the subsequent run
236!
237!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
238
239  USE chem_gasphase_mod_Parameters
240  USE chem_gasphase_mod_LinearAlgebra
241  IMPLICIT NONE
242
243!~~~>  Arguments
244   INTEGER,       INTENT(IN)    :: N
245   REAL(kind=dp), INTENT(INOUT) :: Y(N)
246   REAL(kind=dp), INTENT(IN)    :: Tstart,Tend
247   REAL(kind=dp), INTENT(IN)    :: AbsTol(N),RelTol(N)
248   INTEGER,       INTENT(IN)    :: ICNTRL(20)
249   REAL(kind=dp), INTENT(IN)    :: RCNTRL(20)
250   INTEGER,       INTENT(INOUT) :: ISTATUS(20)
251   REAL(kind=dp), INTENT(INOUT) :: RSTATUS(20)
252   INTEGER, INTENT(OUT)   :: IERR
253!~~~>  Parameters of the Rosenbrock method, up to 6 stages
254   INTEGER ::  ros_S, rosMethod
255   INTEGER, PARAMETER :: RS2=1, RS3=2, RS4=3, RD3=4, RD4=5, RG3=6
256   REAL(kind=dp) :: ros_A(15), ros_C(15), ros_M(6), ros_E(6), &
257                    ros_Alpha(6), ros_Gamma(6), ros_ELO
258   LOGICAL :: ros_NewF(6)
259   CHARACTER(LEN=12) :: ros_Name
260!~~~>  Local variables
261   REAL(kind=dp) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
262   REAL(kind=dp) :: Hmin, Hmax, Hstart
263   REAL(kind=dp) :: Texit
264   INTEGER       :: i, UplimTol, Max_no_steps
265   LOGICAL       :: Autonomous, VectorTol
266!~~~>   Parameters
267   REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp, ONE  = 1.0_dp
268   REAL(kind=dp), PARAMETER :: DeltaMin = 1.0E-5_dp
269
270!~~~>  Initialize statistics
271   ISTATUS(1:8) = 0
272   RSTATUS(1:3) = ZERO
273
274!~~~>  Autonomous or time dependent ODE. Default is time dependent.
275   Autonomous = .NOT.(ICNTRL(1) == 0)
276
277!~~~>  For Scalar tolerances (ICNTRL(2).NE.0)  the code uses AbsTol(1) and RelTol(1)
278!   For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:N) and RelTol(1:N)
279   IF (ICNTRL(2) == 0) THEN
280      VectorTol = .TRUE.
281      UplimTol  = N
282   ELSE
283      VectorTol = .FALSE.
284      UplimTol  = 1
285   END IF
286
287!~~~>   Initialize the particular Rosenbrock method selected
288   SELECT CASE (ICNTRL(3))
289     CASE (1)
290       CALL Ros2
291     CASE (2)
292       CALL Ros3
293     CASE (3)
294       CALL Ros4
295     CASE (0,4)
296       CALL Rodas3
297     CASE (5)
298       CALL Rodas4
299     CASE (6)
300       CALL Rang3
301     CASE DEFAULT
302       PRINT * , 'Unknown Rosenbrock method: ICNTRL(3)=',ICNTRL(3) 
303       CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
304       RETURN
305   END SELECT
306
307!~~~>   The maximum number of steps admitted
308   IF (ICNTRL(4) == 0) THEN
309      Max_no_steps = 200000
310   ELSEIF (ICNTRL(4) > 0) THEN
311      Max_no_steps=ICNTRL(4)
312   ELSE
313      PRINT * ,'User-selected max no. of steps: ICNTRL(4)=',ICNTRL(4)
314      CALL ros_ErrorMsg(-1,Tstart,ZERO,IERR)
315      RETURN
316   END IF
317
318!~~~>  Unit roundoff (1+Roundoff>1)
319   Roundoff = WLAMCH('E')
320
321!~~~>  Lower bound on the step size: (positive value)
322   IF (RCNTRL(1) == ZERO) THEN
323      Hmin = ZERO
324   ELSEIF (RCNTRL(1) > ZERO) THEN
325      Hmin = RCNTRL(1)
326   ELSE
327      PRINT * , 'User-selected Hmin: RCNTRL(1)=', RCNTRL(1)
328      CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
329      RETURN
330   END IF
331!~~~>  Upper bound on the step size: (positive value)
332   IF (RCNTRL(2) == ZERO) THEN
333      Hmax = ABS(Tend-Tstart)
334   ELSEIF (RCNTRL(2) > ZERO) THEN
335      Hmax = MIN(ABS(RCNTRL(2)),ABS(Tend-Tstart))
336   ELSE
337      PRINT * , 'User-selected Hmax: RCNTRL(2)=', RCNTRL(2)
338      CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
339      RETURN
340   END IF
341!~~~>  Starting step size: (positive value)
342   IF (RCNTRL(3) == ZERO) THEN
343      Hstart = MAX(Hmin,DeltaMin)
344   ELSEIF (RCNTRL(3) > ZERO) THEN
345      Hstart = MIN(ABS(RCNTRL(3)),ABS(Tend-Tstart))
346   ELSE
347      PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3)
348      CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
349      RETURN
350   END IF
351!~~~>  Step size can be changed s.t.  FacMin < Hnew/Hold < FacMax
352   IF (RCNTRL(4) == ZERO) THEN
353      FacMin = 0.2_dp
354   ELSEIF (RCNTRL(4) > ZERO) THEN
355      FacMin = RCNTRL(4)
356   ELSE
357      PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4)
358      CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
359      RETURN
360   END IF
361   IF (RCNTRL(5) == ZERO) THEN
362      FacMax = 6.0_dp
363   ELSEIF (RCNTRL(5) > ZERO) THEN
364      FacMax = RCNTRL(5)
365   ELSE
366      PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5)
367      CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
368      RETURN
369   END IF
370!~~~>   FacRej: Factor to decrease step after 2 succesive rejections
371   IF (RCNTRL(6) == ZERO) THEN
372      FacRej = 0.1_dp
373   ELSEIF (RCNTRL(6) > ZERO) THEN
374      FacRej = RCNTRL(6)
375   ELSE
376      PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6)
377      CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
378      RETURN
379   END IF
380!~~~>   FacSafe: Safety Factor in the computation of new step size
381   IF (RCNTRL(7) == ZERO) THEN
382      FacSafe = 0.9_dp
383   ELSEIF (RCNTRL(7) > ZERO) THEN
384      FacSafe = RCNTRL(7)
385   ELSE
386      PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7)
387      CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
388      RETURN
389   END IF
390!~~~>  Check if tolerances are reasonable
391    DO i=1,UplimTol
392      IF ( (AbsTol(i) <= ZERO) .OR. (RelTol(i) <= 10.0_dp*Roundoff) &
393         .OR. (RelTol(i) >= 1.0_dp) ) THEN
394        PRINT * , ' AbsTol(',i,') = ',AbsTol(i)
395        PRINT * , ' RelTol(',i,') = ',RelTol(i)
396        CALL ros_ErrorMsg(-5,Tstart,ZERO,IERR)
397        RETURN
398      END IF
399    END DO
400
401
402!~~~>  CALL Rosenbrock method
403   CALL ros_Integrator(Y, Tstart, Tend, Texit,   &
404        AbsTol, RelTol,                          &
405!  Integration parameters
406        Autonomous, VectorTol, Max_no_steps,     &
407        Roundoff, Hmin, Hmax, Hstart,            &
408        FacMin, FacMax, FacRej, FacSafe,         &
409!  Error indicator
410        IERR)
411
412!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
413CONTAINS !  SUBROUTINES internal to Rosenbrock
414!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
415   
416!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
417 SUBROUTINE ros_ErrorMsg(Code,T,H,IERR)
418!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
419!    Handles all error messages
420!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
421 
422   REAL(kind=dp), INTENT(IN) :: T, H
423   INTEGER, INTENT(IN)  :: Code
424   INTEGER, INTENT(OUT) :: IERR
425   
426   IERR = Code
427   PRINT * , &
428     'Forced exit from Rosenbrock due to the following error:' 
429     
430   SELECT CASE (Code)
431    CASE (-1)   
432      PRINT * , '--> Improper value for maximal no of steps'
433    CASE (-2)   
434      PRINT * , '--> Selected Rosenbrock method not implemented'
435    CASE (-3)   
436      PRINT * , '--> Hmin/Hmax/Hstart must be positive'
437    CASE (-4)   
438      PRINT * , '--> FacMin/FacMax/FacRej must be positive'
439    CASE (-5) 
440      PRINT * , '--> Improper tolerance values'
441    CASE (-6) 
442      PRINT * , '--> No of steps exceeds maximum bound'
443    CASE (-7) 
444      PRINT * , '--> Step size too small: T + 10*H = T', &
445            ' or H < Roundoff'
446    CASE (-8)   
447      PRINT * , '--> Matrix is repeatedly singular'
448    CASE DEFAULT
449      PRINT *, 'Unknown Error code: ', Code
450   END SELECT
451   
452   PRINT *, "T=", T, "and H=", H
453     
454 END SUBROUTINE ros_ErrorMsg
455
456!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
457 SUBROUTINE ros_Integrator (Y, Tstart, Tend, T,  &
458        AbsTol, RelTol,                          &
459!~~~> Integration parameters
460        Autonomous, VectorTol, Max_no_steps,     &
461        Roundoff, Hmin, Hmax, Hstart,            &
462        FacMin, FacMax, FacRej, FacSafe,         &
463!~~~> Error indicator
464        IERR )
465!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
466!   Template for the implementation of a generic Rosenbrock method
467!      defined by ros_S (no of stages)
468!      and its coefficients ros_{A,C,M,E,Alpha,Gamma}
469!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
470
471  IMPLICIT NONE
472
473!~~~> Input: the initial condition at Tstart; Output: the solution at T
474   REAL(kind=dp), INTENT(INOUT) :: Y(N)
475!~~~> Input: integration interval
476   REAL(kind=dp), INTENT(IN) :: Tstart,Tend
477!~~~> Output: time at which the solution is returned (T=Tend if success)
478   REAL(kind=dp), INTENT(OUT) ::  T
479!~~~> Input: tolerances
480   REAL(kind=dp), INTENT(IN) ::  AbsTol(N), RelTol(N)
481!~~~> Input: integration parameters
482   LOGICAL, INTENT(IN) :: Autonomous, VectorTol
483   REAL(kind=dp), INTENT(IN) :: Hstart, Hmin, Hmax
484   INTEGER, INTENT(IN) :: Max_no_steps
485   REAL(kind=dp), INTENT(IN) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
486!~~~> Output: Error indicator
487   INTEGER, INTENT(OUT) :: IERR
488! ~~~~ Local variables
489   REAL(kind=dp) :: Ynew(N), Fcn0(N), Fcn(N)
490   REAL(kind=dp) :: K(N*ros_S), dFdT(N)
491#ifdef FULL_ALGEBRA   
492   REAL(kind=dp) :: Jac0(N,N), Ghimj(N,N)
493#else
494   REAL(kind=dp) :: Jac0(LU_NONZERO), Ghimj(LU_NONZERO)
495#endif
496   REAL(kind=dp) :: H, Hnew, HC, HG, Fac, Tau
497   REAL(kind=dp) :: Err, Yerr(N)
498   INTEGER :: Pivot(N), Direction, ioffset, j, istage
499   LOGICAL :: RejectLastH, RejectMoreH, Singular
500!~~~>  Local parameters
501   REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp, ONE  = 1.0_dp
502   REAL(kind=dp), PARAMETER :: DeltaMin = 1.0E-5_dp
503!~~~>  Locally called functions
504!    REAL(kind=dp) WLAMCH
505!    EXTERNAL WLAMCH
506!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
507
508
509!~~~>  Initial preparations
510   T = Tstart
511   RSTATUS(Nhexit) = ZERO
512   H = MIN( MAX(ABS(Hmin),ABS(Hstart)) , ABS(Hmax) )
513   IF (ABS(H) <= 10.0_dp*Roundoff) H = DeltaMin
514
515   IF (Tend  >=  Tstart) THEN
516     Direction = +1
517   ELSE
518     Direction = -1
519   END IF
520   H = Direction*H
521
522   RejectLastH=.FALSE.
523   RejectMoreH=.FALSE.
524
525!~~~> Time loop begins below
526
527TimeLoop: DO WHILE ( (Direction > 0).AND.((T-Tend)+Roundoff <= ZERO) &
528       .OR. (Direction < 0).AND.((Tend-T)+Roundoff <= ZERO) )
529
530   IF ( ISTATUS(Nstp) > Max_no_steps ) THEN  ! Too many steps
531      CALL ros_ErrorMsg(-6,T,H,IERR)
532      RETURN
533   END IF
534   IF ( ((T+0.1_dp*H) == T).OR.(H <= Roundoff) ) THEN  ! Step size too small
535      CALL ros_ErrorMsg(-7,T,H,IERR)
536      RETURN
537   END IF
538
539!~~~>  Limit H if necessary to avoid going beyond Tend
540   H = MIN(H,ABS(Tend-T))
541
542!~~~>   Compute the function at current time
543   CALL FunTemplate(T,Y,Fcn0)
544   ISTATUS(Nfun) = ISTATUS(Nfun) + 1
545
546!~~~>  Compute the function derivative with respect to T
547   IF (.NOT.Autonomous) THEN
548      CALL ros_FunTimeDerivative ( T, Roundoff, Y, &
549                Fcn0, dFdT )
550   END IF
551
552!~~~>   Compute the Jacobian at current time
553   CALL JacTemplate(T,Y,Jac0)
554   ISTATUS(Njac) = ISTATUS(Njac) + 1
555
556!~~~>  Repeat step calculation until current step accepted
557UntilAccepted: DO
558
559   CALL ros_PrepareMatrix(H,Direction,ros_Gamma(1), &
560          Jac0,Ghimj,Pivot,Singular)
561   IF (Singular) THEN ! More than 5 consecutive failed decompositions
562       CALL ros_ErrorMsg(-8,T,H,IERR)
563       RETURN
564   END IF
565
566!~~~>   Compute the stages
567Stage: DO istage = 1, ros_S
568
569      ! Current istage offset. Current istage vector is K(ioffset+1:ioffset+N)
570       ioffset = N*(istage-1)
571
572      ! For the 1st istage the function has been computed previously
573       IF ( istage == 1 ) THEN
574         !slim: CALL WCOPY(N,Fcn0,1,Fcn,1)
575         Fcn(1:N) = Fcn0(1:N)
576      ! istage>1 and a new function evaluation is needed at the current istage
577       ELSEIF ( ros_NewF(istage) ) THEN
578         !slim: CALL WCOPY(N,Y,1,Ynew,1)
579         Ynew(1:N) = Y(1:N)
580         DO j = 1, istage-1
581           CALL WAXPY(N,ros_A((istage-1)*(istage-2)/2+j), &
582            K(N*(j-1)+1),1,Ynew,1)
583         END DO
584         Tau = T + ros_Alpha(istage)*Direction*H
585         CALL FunTemplate(Tau,Ynew,Fcn)
586         ISTATUS(Nfun) = ISTATUS(Nfun) + 1
587       END IF ! if istage == 1 elseif ros_NewF(istage)
588       !slim: CALL WCOPY(N,Fcn,1,K(ioffset+1),1)
589       K(ioffset+1:ioffset+N) = Fcn(1:N)
590       DO j = 1, istage-1
591         HC = ros_C((istage-1)*(istage-2)/2+j)/(Direction*H)
592         CALL WAXPY(N,HC,K(N*(j-1)+1),1,K(ioffset+1),1)
593       END DO
594       IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN
595         HG = Direction*H*ros_Gamma(istage)
596         CALL WAXPY(N,HG,dFdT,1,K(ioffset+1),1)
597       END IF
598       CALL ros_Solve(Ghimj, Pivot, K(ioffset+1))
599
600   END DO Stage
601
602
603!~~~>  Compute the new solution
604   !slim: CALL WCOPY(N,Y,1,Ynew,1)
605   Ynew(1:N) = Y(1:N)
606   DO j=1,ros_S
607         CALL WAXPY(N,ros_M(j),K(N*(j-1)+1),1,Ynew,1)
608   END DO
609
610!~~~>  Compute the error estimation
611   !slim: CALL WSCAL(N,ZERO,Yerr,1)
612   Yerr(1:N) = ZERO
613   DO j=1,ros_S
614        CALL WAXPY(N,ros_E(j),K(N*(j-1)+1),1,Yerr,1)
615   END DO
616   Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol )
617
618!~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax
619   Fac  = MIN(FacMax,MAX(FacMin,FacSafe/Err**(ONE/ros_ELO)))
620   Hnew = H*Fac
621
622!~~~>  Check the error magnitude and adjust step size
623   ISTATUS(Nstp) = ISTATUS(Nstp) + 1
624   IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN  !~~~> Accept step
625      ISTATUS(Nacc) = ISTATUS(Nacc) + 1
626      !slim: CALL WCOPY(N,Ynew,1,Y,1)
627      Y(1:N) = Ynew(1:N)
628      T = T + Direction*H
629      Hnew = MAX(Hmin,MIN(Hnew,Hmax))
630      IF (RejectLastH) THEN  ! No step size increase after a rejected step
631         Hnew = MIN(Hnew,H)
632      END IF
633      RSTATUS(Nhexit) = H
634      RSTATUS(Nhnew)  = Hnew
635      RSTATUS(Ntexit) = T
636      RejectLastH = .FALSE.
637      RejectMoreH = .FALSE.
638      H = Hnew
639      EXIT UntilAccepted ! EXIT THE LOOP: WHILE STEP NOT ACCEPTED
640   ELSE           !~~~> Reject step
641      IF (RejectMoreH) THEN
642         Hnew = H*FacRej
643      END IF
644      RejectMoreH = RejectLastH
645      RejectLastH = .TRUE.
646      H = Hnew
647      IF (ISTATUS(Nacc) >= 1)  ISTATUS(Nrej) = ISTATUS(Nrej) + 1
648   END IF ! Err <= 1
649
650   END DO UntilAccepted
651
652   END DO TimeLoop
653
654!~~~> Succesful exit
655   IERR = 1  !~~~> The integration was successful
656
657  END SUBROUTINE ros_Integrator
658
659
660!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
661  REAL(kind=dp) FUNCTION ros_ErrorNorm ( Y, Ynew, Yerr, &
662                               AbsTol, RelTol, VectorTol )
663!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
664!~~~> Computes the "scaled norm" of the error vector Yerr
665!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
666   IMPLICIT NONE
667
668! Input arguments
669   REAL(kind=dp), INTENT(IN) :: Y(N), Ynew(N), &
670          Yerr(N), AbsTol(N), RelTol(N)
671   LOGICAL, INTENT(IN) ::  VectorTol
672! Local variables
673   REAL(kind=dp) :: Err, Scale, Ymax
674   INTEGER  :: i
675   REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp
676
677   Err = ZERO
678   DO i=1,N
679     Ymax = MAX(ABS(Y(i)),ABS(Ynew(i)))
680     IF (VectorTol) THEN
681       Scale = AbsTol(i)+RelTol(i)*Ymax
682     ELSE
683       Scale = AbsTol(1)+RelTol(1)*Ymax
684     END IF
685     Err = Err+(Yerr(i)/Scale)**2
686   END DO
687   Err  = SQRT(Err/N)
688
689   ros_ErrorNorm = MAX(Err,1.0d-10)
690
691  END FUNCTION ros_ErrorNorm
692
693
694!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
695  SUBROUTINE ros_FunTimeDerivative ( T, Roundoff, Y, &
696                Fcn0, dFdT )
697!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
698!~~~> The time partial derivative of the function by finite differences
699!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
700   IMPLICIT NONE
701
702!~~~> Input arguments
703   REAL(kind=dp), INTENT(IN) :: T, Roundoff, Y(N), Fcn0(N)
704!~~~> Output arguments
705   REAL(kind=dp), INTENT(OUT) :: dFdT(N)
706!~~~> Local variables
707   REAL(kind=dp) :: Delta
708   REAL(kind=dp), PARAMETER :: ONE = 1.0_dp, DeltaMin = 1.0E-6_dp
709
710   Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T))
711   CALL FunTemplate(T+Delta,Y,dFdT)
712   ISTATUS(Nfun) = ISTATUS(Nfun) + 1
713   CALL WAXPY(N,(-ONE),Fcn0,1,dFdT,1)
714   CALL WSCAL(N,(ONE/Delta),dFdT,1)
715
716  END SUBROUTINE ros_FunTimeDerivative
717
718
719!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
720  SUBROUTINE ros_PrepareMatrix ( H, Direction, gam, &
721             Jac0, Ghimj, Pivot, Singular )
722! --- --- --- --- --- --- --- --- --- --- --- --- ---
723!  Prepares the LHS matrix for stage calculations
724!  1.  Construct Ghimj = 1/(H*ham) - Jac0
725!      "(Gamma H) Inverse Minus Jacobian"
726!  2.  Repeat LU decomposition of Ghimj until successful.
727!       -half the step size if LU decomposition fails and retry
728!       -exit after 5 consecutive fails
729! --- --- --- --- --- --- --- --- --- --- --- --- ---
730   IMPLICIT NONE
731
732!~~~> Input arguments
733#ifdef FULL_ALGEBRA   
734   REAL(kind=dp), INTENT(IN) ::  Jac0(N,N)
735#else
736   REAL(kind=dp), INTENT(IN) ::  Jac0(LU_NONZERO)
737#endif   
738   REAL(kind=dp), INTENT(IN) ::  gam
739   INTEGER, INTENT(IN) ::  Direction
740!~~~> Output arguments
741#ifdef FULL_ALGEBRA   
742   REAL(kind=dp), INTENT(OUT) :: Ghimj(N,N)
743#else
744   REAL(kind=dp), INTENT(OUT) :: Ghimj(LU_NONZERO)
745#endif   
746   LOGICAL, INTENT(OUT) ::  Singular
747   INTEGER, INTENT(OUT) ::  Pivot(N)
748!~~~> Inout arguments
749   REAL(kind=dp), INTENT(INOUT) :: H   ! step size is decreased when LU fails
750!~~~> Local variables
751   INTEGER  :: i, ISING, Nconsecutive
752   REAL(kind=dp) :: ghinv
753   REAL(kind=dp), PARAMETER :: ONE  = 1.0_dp, HALF = 0.5_dp
754
755   Nconsecutive = 0
756   Singular = .TRUE.
757
758   DO WHILE (Singular)
759
760!~~~>    Construct Ghimj = 1/(H*gam) - Jac0
761#ifdef FULL_ALGEBRA   
762     !slim: CALL WCOPY(N*N,Jac0,1,Ghimj,1)
763     !slim: CALL WSCAL(N*N,(-ONE),Ghimj,1)
764     Ghimj = -Jac0
765     ghinv = ONE/(Direction*H*gam)
766     DO i=1,N
767       Ghimj(i,i) = Ghimj(i,i)+ghinv
768     END DO
769#else
770     !slim: CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1)
771     !slim: CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1)
772     Ghimj(1:LU_NONZERO) = -Jac0(1:LU_NONZERO)
773     ghinv = ONE/(Direction*H*gam)
774     DO i=1,N
775       Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+ghinv
776     END DO
777#endif   
778!~~~>    Compute LU decomposition
779     CALL ros_Decomp( Ghimj, Pivot, ISING )
780     IF (ISING == 0) THEN
781!~~~>    If successful done
782        Singular = .FALSE.
783     ELSE ! ISING .ne. 0
784!~~~>    If unsuccessful half the step size; if 5 consecutive fails then return
785        ISTATUS(Nsng) = ISTATUS(Nsng) + 1
786        Nconsecutive = Nconsecutive+1
787        Singular = .TRUE.
788        PRINT*,'Warning: LU Decomposition returned ISING = ',ISING
789        IF (Nconsecutive <= 5) THEN ! Less than 5 consecutive failed decompositions
790           H = H*HALF
791        ELSE  ! More than 5 consecutive failed decompositions
792           RETURN
793        END IF  ! Nconsecutive
794      END IF    ! ISING
795
796   END DO ! WHILE Singular
797
798  END SUBROUTINE ros_PrepareMatrix
799
800
801!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
802  SUBROUTINE ros_Decomp( A, Pivot, ISING )
803!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
804!  Template for the LU decomposition
805!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
806   IMPLICIT NONE
807!~~~> Inout variables
808#ifdef FULL_ALGEBRA   
809   REAL(kind=dp), INTENT(INOUT) :: A(N,N)
810#else   
811   REAL(kind=dp), INTENT(INOUT) :: A(LU_NONZERO)
812#endif
813!~~~> Output variables
814   INTEGER, INTENT(OUT) :: Pivot(N), ISING
815
816#ifdef FULL_ALGEBRA   
817   CALL  DGETRF( N, N, A, N, Pivot, ISING )
818#else   
819   CALL KppDecomp ( A, ISING )
820   Pivot(1) = 1
821#endif
822   ISTATUS(Ndec) = ISTATUS(Ndec) + 1
823
824  END SUBROUTINE ros_Decomp
825
826
827!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
828  SUBROUTINE ros_Solve( A, Pivot, b )
829!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
830!  Template for the forward/backward substitution (using pre-computed LU decomposition)
831!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
832   IMPLICIT NONE
833!~~~> Input variables
834#ifdef FULL_ALGEBRA   
835   REAL(kind=dp), INTENT(IN) :: A(N,N)
836   INTEGER :: ISING
837#else   
838   REAL(kind=dp), INTENT(IN) :: A(LU_NONZERO)
839#endif
840   INTEGER, INTENT(IN) :: Pivot(N)
841!~~~> InOut variables
842   REAL(kind=dp), INTENT(INOUT) :: b(N)
843
844#ifdef FULL_ALGEBRA   
845   CALL  DGETRS( 'N', N , 1, A, N, Pivot, b, N, ISING )
846   IF ( Info < 0 ) THEN
847      PRINT*,"Error in DGETRS. ISING=",ISING
848   END IF 
849#else   
850   CALL KppSolve( A, b )
851#endif
852
853   ISTATUS(Nsol) = ISTATUS(Nsol) + 1
854
855  END SUBROUTINE ros_Solve
856
857
858
859!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
860  SUBROUTINE Ros2
861!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
862! --- AN L-STABLE METHOD, 2 stages, order 2
863!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
864
865   IMPLICIT NONE
866   DOUBLE PRECISION g
867
868    g = 1.0_dp + 1.0_dp/SQRT(2.0_dp)
869    rosMethod = RS2
870!~~~> Name of the method
871    ros_Name = 'ROS-2'
872!~~~> Number of stages
873    ros_S = 2
874
875!~~~> The coefficient matrices A and C are strictly lower triangular.
876!   The lower triangular (subdiagonal) elements are stored in row-wise order:
877!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
878!   The general mapping formula is:
879!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
880!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
881
882    ros_A(1) = (1.0_dp)/g
883    ros_C(1) = (-2.0_dp)/g
884!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
885!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
886    ros_NewF(1) = .TRUE.
887    ros_NewF(2) = .TRUE.
888!~~~> M_i = Coefficients for new step solution
889    ros_M(1)= (3.0_dp)/(2.0_dp*g)
890    ros_M(2)= (1.0_dp)/(2.0_dp*g)
891! E_i = Coefficients for error estimator
892    ros_E(1) = 1.0_dp/(2.0_dp*g)
893    ros_E(2) = 1.0_dp/(2.0_dp*g)
894!~~~> ros_ELO = estimator of local order - the minimum between the
895!    main and the embedded scheme orders plus one
896    ros_ELO = 2.0_dp
897!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
898    ros_Alpha(1) = 0.0_dp
899    ros_Alpha(2) = 1.0_dp
900!~~~> Gamma_i = \sum_j  gamma_{i,j}
901    ros_Gamma(1) = g
902    ros_Gamma(2) =-g
903
904 END SUBROUTINE Ros2
905
906
907!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
908  SUBROUTINE Ros3
909!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
910! --- AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations
911!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
912
913   IMPLICIT NONE
914   rosMethod = RS3
915!~~~> Name of the method
916   ros_Name = 'ROS-3'
917!~~~> Number of stages
918   ros_S = 3
919
920!~~~> The coefficient matrices A and C are strictly lower triangular.
921!   The lower triangular (subdiagonal) elements are stored in row-wise order:
922!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
923!   The general mapping formula is:
924!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
925!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
926
927   ros_A(1)= 1.0_dp
928   ros_A(2)= 1.0_dp
929   ros_A(3)= 0.0_dp
930
931   ros_C(1) = -0.10156171083877702091975600115545E+01_dp
932   ros_C(2) =  0.40759956452537699824805835358067E+01_dp
933   ros_C(3) =  0.92076794298330791242156818474003E+01_dp
934!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
935!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
936   ros_NewF(1) = .TRUE.
937   ros_NewF(2) = .TRUE.
938   ros_NewF(3) = .FALSE.
939!~~~> M_i = Coefficients for new step solution
940   ros_M(1) =  0.1E+01_dp
941   ros_M(2) =  0.61697947043828245592553615689730E+01_dp
942   ros_M(3) = -0.42772256543218573326238373806514_dp
943! E_i = Coefficients for error estimator
944   ros_E(1) =  0.5_dp
945   ros_E(2) = -0.29079558716805469821718236208017E+01_dp
946   ros_E(3) =  0.22354069897811569627360909276199_dp
947!~~~> ros_ELO = estimator of local order - the minimum between the
948!    main and the embedded scheme orders plus 1
949   ros_ELO = 3.0_dp
950!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
951   ros_Alpha(1)= 0.0_dp
952   ros_Alpha(2)= 0.43586652150845899941601945119356_dp
953   ros_Alpha(3)= 0.43586652150845899941601945119356_dp
954!~~~> Gamma_i = \sum_j  gamma_{i,j}
955   ros_Gamma(1)= 0.43586652150845899941601945119356_dp
956   ros_Gamma(2)= 0.24291996454816804366592249683314_dp
957   ros_Gamma(3)= 0.21851380027664058511513169485832E+01_dp
958
959  END SUBROUTINE Ros3
960
961!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
962
963
964!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
965  SUBROUTINE Ros4
966!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
967!     L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES
968!     L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
969!
970!      E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
971!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
972!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
973!      SPRINGER-VERLAG (1990)
974!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
975
976   IMPLICIT NONE
977
978   rosMethod = RS4
979!~~~> Name of the method
980   ros_Name = 'ROS-4'
981!~~~> Number of stages
982   ros_S = 4
983
984!~~~> The coefficient matrices A and C are strictly lower triangular.
985!   The lower triangular (subdiagonal) elements are stored in row-wise order:
986!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
987!   The general mapping formula is:
988!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
989!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
990
991   ros_A(1) = 0.2000000000000000E+01_dp
992   ros_A(2) = 0.1867943637803922E+01_dp
993   ros_A(3) = 0.2344449711399156_dp
994   ros_A(4) = ros_A(2)
995   ros_A(5) = ros_A(3)
996   ros_A(6) = 0.0_dp
997
998   ros_C(1) =-0.7137615036412310E+01_dp
999   ros_C(2) = 0.2580708087951457E+01_dp
1000   ros_C(3) = 0.6515950076447975_dp
1001   ros_C(4) =-0.2137148994382534E+01_dp
1002   ros_C(5) =-0.3214669691237626_dp
1003   ros_C(6) =-0.6949742501781779_dp
1004!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1005!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1006   ros_NewF(1)  = .TRUE.
1007   ros_NewF(2)  = .TRUE.
1008   ros_NewF(3)  = .TRUE.
1009   ros_NewF(4)  = .FALSE.
1010!~~~> M_i = Coefficients for new step solution
1011   ros_M(1) = 0.2255570073418735E+01_dp
1012   ros_M(2) = 0.2870493262186792_dp
1013   ros_M(3) = 0.4353179431840180_dp
1014   ros_M(4) = 0.1093502252409163E+01_dp
1015!~~~> E_i  = Coefficients for error estimator
1016   ros_E(1) =-0.2815431932141155_dp
1017   ros_E(2) =-0.7276199124938920E-01_dp
1018   ros_E(3) =-0.1082196201495311_dp
1019   ros_E(4) =-0.1093502252409163E+01_dp
1020!~~~> ros_ELO  = estimator of local order - the minimum between the
1021!    main and the embedded scheme orders plus 1
1022   ros_ELO  = 4.0_dp
1023!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1024   ros_Alpha(1) = 0.0_dp
1025   ros_Alpha(2) = 0.1145640000000000E+01_dp
1026   ros_Alpha(3) = 0.6552168638155900_dp
1027   ros_Alpha(4) = ros_Alpha(3)
1028!~~~> Gamma_i = \sum_j  gamma_{i,j}
1029   ros_Gamma(1) = 0.5728200000000000_dp
1030   ros_Gamma(2) =-0.1769193891319233E+01_dp
1031   ros_Gamma(3) = 0.7592633437920482_dp
1032   ros_Gamma(4) =-0.1049021087100450_dp
1033
1034  END SUBROUTINE Ros4
1035
1036!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1037  SUBROUTINE Rodas3
1038!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1039! --- A STIFFLY-STABLE METHOD, 4 stages, order 3
1040!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1041
1042   IMPLICIT NONE
1043
1044   rosMethod = RD3
1045!~~~> Name of the method
1046   ros_Name = 'RODAS-3'
1047!~~~> Number of stages
1048   ros_S = 4
1049
1050!~~~> The coefficient matrices A and C are strictly lower triangular.
1051!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1052!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1053!   The general mapping formula is:
1054!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
1055!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
1056
1057   ros_A(1) = 0.0_dp
1058   ros_A(2) = 2.0_dp
1059   ros_A(3) = 0.0_dp
1060   ros_A(4) = 2.0_dp
1061   ros_A(5) = 0.0_dp
1062   ros_A(6) = 1.0_dp
1063
1064   ros_C(1) = 4.0_dp
1065   ros_C(2) = 1.0_dp
1066   ros_C(3) =-1.0_dp
1067   ros_C(4) = 1.0_dp
1068   ros_C(5) =-1.0_dp
1069   ros_C(6) =-(8.0_dp/3.0_dp)
1070
1071!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1072!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1073   ros_NewF(1)  = .TRUE.
1074   ros_NewF(2)  = .FALSE.
1075   ros_NewF(3)  = .TRUE.
1076   ros_NewF(4)  = .TRUE.
1077!~~~> M_i = Coefficients for new step solution
1078   ros_M(1) = 2.0_dp
1079   ros_M(2) = 0.0_dp
1080   ros_M(3) = 1.0_dp
1081   ros_M(4) = 1.0_dp
1082!~~~> E_i  = Coefficients for error estimator
1083   ros_E(1) = 0.0_dp
1084   ros_E(2) = 0.0_dp
1085   ros_E(3) = 0.0_dp
1086   ros_E(4) = 1.0_dp
1087!~~~> ros_ELO  = estimator of local order - the minimum between the
1088!    main and the embedded scheme orders plus 1
1089   ros_ELO  = 3.0_dp
1090!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1091   ros_Alpha(1) = 0.0_dp
1092   ros_Alpha(2) = 0.0_dp
1093   ros_Alpha(3) = 1.0_dp
1094   ros_Alpha(4) = 1.0_dp
1095!~~~> Gamma_i = \sum_j  gamma_{i,j}
1096   ros_Gamma(1) = 0.5_dp
1097   ros_Gamma(2) = 1.5_dp
1098   ros_Gamma(3) = 0.0_dp
1099   ros_Gamma(4) = 0.0_dp
1100
1101  END SUBROUTINE Rodas3
1102
1103!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1104  SUBROUTINE Rodas4
1105!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1106!     STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES
1107!
1108!      E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
1109!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
1110!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
1111!      SPRINGER-VERLAG (1996)
1112!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1113
1114   IMPLICIT NONE
1115
1116    rosMethod = RD4
1117!~~~> Name of the method
1118    ros_Name = 'RODAS-4'
1119!~~~> Number of stages
1120    ros_S = 6
1121
1122!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1123    ros_Alpha(1) = 0.000_dp
1124    ros_Alpha(2) = 0.386_dp
1125    ros_Alpha(3) = 0.210_dp
1126    ros_Alpha(4) = 0.630_dp
1127    ros_Alpha(5) = 1.000_dp
1128    ros_Alpha(6) = 1.000_dp
1129
1130!~~~> Gamma_i = \sum_j  gamma_{i,j}
1131    ros_Gamma(1) = 0.2500000000000000_dp
1132    ros_Gamma(2) =-0.1043000000000000_dp
1133    ros_Gamma(3) = 0.1035000000000000_dp
1134    ros_Gamma(4) =-0.3620000000000023E-01_dp
1135    ros_Gamma(5) = 0.0_dp
1136    ros_Gamma(6) = 0.0_dp
1137
1138!~~~> The coefficient matrices A and C are strictly lower triangular.
1139!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1140!   A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1141!   The general mapping formula is:  A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
1142!                  C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
1143
1144    ros_A(1) = 0.1544000000000000E+01_dp
1145    ros_A(2) = 0.9466785280815826_dp
1146    ros_A(3) = 0.2557011698983284_dp
1147    ros_A(4) = 0.3314825187068521E+01_dp
1148    ros_A(5) = 0.2896124015972201E+01_dp
1149    ros_A(6) = 0.9986419139977817_dp
1150    ros_A(7) = 0.1221224509226641E+01_dp
1151    ros_A(8) = 0.6019134481288629E+01_dp
1152    ros_A(9) = 0.1253708332932087E+02_dp
1153    ros_A(10) =-0.6878860361058950_dp
1154    ros_A(11) = ros_A(7)
1155    ros_A(12) = ros_A(8)
1156    ros_A(13) = ros_A(9)
1157    ros_A(14) = ros_A(10)
1158    ros_A(15) = 1.0_dp
1159
1160    ros_C(1) =-0.5668800000000000E+01_dp
1161    ros_C(2) =-0.2430093356833875E+01_dp
1162    ros_C(3) =-0.2063599157091915_dp
1163    ros_C(4) =-0.1073529058151375_dp
1164    ros_C(5) =-0.9594562251023355E+01_dp
1165    ros_C(6) =-0.2047028614809616E+02_dp
1166    ros_C(7) = 0.7496443313967647E+01_dp
1167    ros_C(8) =-0.1024680431464352E+02_dp
1168    ros_C(9) =-0.3399990352819905E+02_dp
1169    ros_C(10) = 0.1170890893206160E+02_dp
1170    ros_C(11) = 0.8083246795921522E+01_dp
1171    ros_C(12) =-0.7981132988064893E+01_dp
1172    ros_C(13) =-0.3152159432874371E+02_dp
1173    ros_C(14) = 0.1631930543123136E+02_dp
1174    ros_C(15) =-0.6058818238834054E+01_dp
1175
1176!~~~> M_i = Coefficients for new step solution
1177    ros_M(1) = ros_A(7)
1178    ros_M(2) = ros_A(8)
1179    ros_M(3) = ros_A(9)
1180    ros_M(4) = ros_A(10)
1181    ros_M(5) = 1.0_dp
1182    ros_M(6) = 1.0_dp
1183
1184!~~~> E_i  = Coefficients for error estimator
1185    ros_E(1) = 0.0_dp
1186    ros_E(2) = 0.0_dp
1187    ros_E(3) = 0.0_dp
1188    ros_E(4) = 0.0_dp
1189    ros_E(5) = 0.0_dp
1190    ros_E(6) = 1.0_dp
1191
1192!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1193!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1194    ros_NewF(1) = .TRUE.
1195    ros_NewF(2) = .TRUE.
1196    ros_NewF(3) = .TRUE.
1197    ros_NewF(4) = .TRUE.
1198    ros_NewF(5) = .TRUE.
1199    ros_NewF(6) = .TRUE.
1200
1201!~~~> ros_ELO  = estimator of local order - the minimum between the
1202!        main and the embedded scheme orders plus 1
1203    ros_ELO = 4.0_dp
1204
1205  END SUBROUTINE Rodas4
1206!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1207  SUBROUTINE Rang3
1208!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1209! STIFFLY-STABLE W METHOD OF ORDER 3, WITH 4 STAGES
1210!
1211! J. RANG and L. ANGERMANN
1212! NEW ROSENBROCK W-METHODS OF ORDER 3
1213! FOR PARTIAL DIFFERENTIAL ALGEBRAIC
1214!        EQUATIONS OF INDEX 1                                             
1215! BIT Numerical Mathematics (2005) 45: 761-787
1216!  DOI: 10.1007/s10543-005-0035-y
1217! Table 4.1-4.2
1218!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1219
1220   IMPLICIT NONE
1221
1222    rosMethod = RG3
1223!~~~> Name of the method
1224    ros_Name = 'RANG-3'
1225!~~~> Number of stages
1226    ros_S = 4
1227
1228    ros_A(1) = 5.09052051067020d+00;
1229    ros_A(2) = 5.09052051067020d+00;
1230    ros_A(3) = 0.0d0;
1231    ros_A(4) = 4.97628111010787d+00;
1232    ros_A(5) = 2.77268164715849d-02;
1233    ros_A(6) = 2.29428036027904d-01;
1234
1235    ros_C(1) = -1.16790812312283d+01;
1236    ros_C(2) = -1.64057326467367d+01;
1237    ros_C(3) = -2.77268164715850d-01;
1238    ros_C(4) = -8.38103960500476d+00;
1239    ros_C(5) = -8.48328409199343d-01;
1240    ros_C(6) =  2.87009860433106d-01;
1241
1242    ros_M(1) =  5.22582761233094d+00;
1243    ros_M(2) = -5.56971148154165d-01;
1244    ros_M(3) =  3.57979469353645d-01;
1245    ros_M(4) =  1.72337398521064d+00;
1246
1247    ros_E(1) = -5.16845212784040d+00;
1248    ros_E(2) = -1.26351942603842d+00;
1249    ros_E(3) = -1.11022302462516d-16;
1250    ros_E(4) =  2.22044604925031d-16;
1251
1252    ros_Alpha(1) = 0.0d00;
1253    ros_Alpha(2) = 2.21878746765329d+00;
1254    ros_Alpha(3) = 2.21878746765329d+00;
1255    ros_Alpha(4) = 1.55392337535788d+00;
1256
1257    ros_Gamma(1) =  4.35866521508459d-01;
1258    ros_Gamma(2) = -1.78292094614483d+00;
1259    ros_Gamma(3) = -2.46541900496934d+00;
1260    ros_Gamma(4) = -8.05529997906370d-01;
1261
1262
1263!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1264!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1265    ros_NewF(1) = .TRUE.
1266    ros_NewF(2) = .TRUE.
1267    ros_NewF(3) = .TRUE.
1268    ros_NewF(4) = .TRUE.
1269
1270!~~~> ros_ELO  = estimator of local order - the minimum between the
1271!        main and the embedded scheme orders plus 1
1272    ros_ELO = 3.0_dp
1273
1274  END SUBROUTINE Rang3
1275!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1276
1277!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1278!   End of the set of internal Rosenbrock subroutines
1279!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1280END SUBROUTINE Rosenbrock
1281!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1282
1283
1284!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1285SUBROUTINE FunTemplate( T, Y, Ydot )
1286!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1287!  Template for the ODE function call.
1288!  Updates the rate coefficients (and possibly the fixed species) at each call
1289!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1290 USE chem_gasphase_mod_Parameters, ONLY: NVAR, LU_NONZERO
1291 USE chem_gasphase_mod_Global, ONLY: FIX, RCONST, TIME
1292 USE chem_gasphase_mod_Function, ONLY: Fun
1293 USE chem_gasphase_mod_Rates, ONLY: Update_SUN, Update_RCONST
1294!~~~> Input variables
1295   REAL(kind=dp) :: T, Y(NVAR)
1296!~~~> Output variables
1297   REAL(kind=dp) :: Ydot(NVAR)
1298!~~~> Local variables
1299   REAL(kind=dp) :: Told
1300
1301   Told = TIME
1302   TIME = T
1303   CALL Update_SUN()
1304   CALL Update_RCONST()
1305   CALL Fun( Y, FIX, RCONST, Ydot )
1306   TIME = Told
1307
1308END SUBROUTINE FunTemplate
1309
1310
1311!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1312SUBROUTINE JacTemplate( T, Y, Jcb )
1313!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1314!  Template for the ODE Jacobian call.
1315!  Updates the rate coefficients (and possibly the fixed species) at each call
1316!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1317 USE chem_gasphase_mod_Parameters, ONLY: NVAR, LU_NONZERO
1318 USE chem_gasphase_mod_Global, ONLY: FIX, RCONST, TIME
1319 USE chem_gasphase_mod_Jacobian, ONLY: Jac_SP, LU_IROW, LU_ICOL
1320 USE chem_gasphase_mod_LinearAlgebra
1321 USE chem_gasphase_mod_Rates, ONLY: Update_SUN, Update_RCONST
1322!~~~> Input variables
1323    REAL(kind=dp) :: T, Y(NVAR)
1324!~~~> Output variables
1325#ifdef FULL_ALGEBRA   
1326    REAL(kind=dp) :: JV(LU_NONZERO), Jcb(NVAR,NVAR)
1327#else
1328    REAL(kind=dp) :: Jcb(LU_NONZERO)
1329#endif   
1330!~~~> Local variables
1331    REAL(kind=dp) :: Told
1332#ifdef FULL_ALGEBRA   
1333    INTEGER :: i, j
1334#endif   
1335
1336    Told = TIME
1337    TIME = T
1338    CALL Update_SUN()
1339    CALL Update_RCONST()
1340#ifdef FULL_ALGEBRA   
1341    CALL Jac_SP(Y, FIX, RCONST, JV)
1342    DO j=1,NVAR
1343      DO i=1,NVAR
1344         Jcb(i,j) = 0.0_dp
1345      END DO
1346    END DO
1347    DO i=1,LU_NONZERO
1348       Jcb(LU_IROW(i),LU_ICOL(i)) = JV(i)
1349    END DO
1350#else
1351    CALL Jac_SP( Y, FIX, RCONST, Jcb )
1352#endif   
1353    TIME = Told
1354
1355END SUBROUTINE JacTemplate
1356
1357END MODULE chem_gasphase_mod_Integrator
1358
1359
1360
1361
1362! End of INTEGRATE function
1363! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1364
1365
Note: See TracBrowser for help on using the repository browser.