source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/rosenbrock_tlm.f90 @ 3698

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

Merge of branch palm4u into trunk

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