Home

Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

rosenbrock_tlm.f90 @ 3997

Last change on this file since 3997 was 2696, checked in by kanani, 6 years ago
Merge of branch palm4u into trunk
File size: 47.2 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
19	MODULE 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
41	CONTAINS ! Functions in the module KPP_ROOT_Integrator
42
43	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
44	SUBROUTINE 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
108	END SUBROUTINE INTEGRATE_TLM
109
110
111	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
112	SUBROUTINE 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
413	CONTAINS ! 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
521	TimeLoop: 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:N,itlm), Fcn0_tlm(1:N,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
561	UntilAccepted: 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
571	Stage: 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:Nj),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:Nj,itlm),1,Ynew_tlm(1:N,itlm),1)
592	END DO
593	END DO
594	Tau = T + ros_Alpha(istage)DirectionH
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:N,itlm), Fcn_tlm(1:N,itlm) )
601	END DO
602	END IF ! if istage == 1 elseif ros_NewF(istage)
603	CALL WCOPY(N,Fcn,1,K(ioffset+1:ioffset+N),1)
604	DO itlm=1,NTLM
605	CALL WCOPY(N,Fcn_tlm(1:N,itlm),1,K_tlm(ioffset+1:ioffset+N,itlm),1)
606	END DO
607	DO j = 1, istage-1
608	HC = ros_C((istage-1)(istage-2)/2+j)/(DirectionH)
609	CALL WAXPY(N,HC,K(N(j-1)+1:Nj),1,K(ioffset+1:ioffset+N),1)
610	DO itlm=1,NTLM
611	CALL WAXPY(N,HC,K_tlm(N(j-1)+1:Nj,itlm),1,K_tlm(ioffset+1:ioffset+N,itlm),1)
612	END DO
613	END DO
614	IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN
615	HG = DirectionHros_Gamma(istage)
616	CALL WAXPY(N,HG,dFdT(1:N),1,K(ioffset+1:ioffset+N),1)
617	DO itlm=1,NTLM
618	CALL Jac_SP_Vec ( dJdT, Ynew_tlm(1:N,itlm), Tmp )
619	CALL WAXPY(N,HG,Tmp,1,K_tlm(ioffset+1:ioffset+N,itlm),1)
620	END DO
621	END IF
622	CALL ros_Solve(Ghimj, Pivot, K(ioffset+1:ioffset+N))
623	DO itlm=1,NTLM
624	CALL Hess_Vec ( Hes0, K(ioffset+1:ioffset+N), Y_tlm(1:N,itlm), Tmp )
625	CALL WAXPY(N,ONE,Tmp,1,K_tlm(ioffset+1:ioffset+N,itlm),1)
626	CALL ros_Solve(Ghimj, Pivot, K_tlm(ioffset+1:ioffset+N,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:Nj),1,Ynew,1)
636	END DO
637	DO itlm=1,NTLM
638	CALL WCOPY(N,Y_tlm(1:N,itlm),1,Ynew_tlm(1:N,itlm),1)
639	DO j=1,ros_S
640	CALL WAXPY(N,ros_M(j),K_tlm(N(j-1)+1:Nj,itlm),1,Ynew_tlm(1:N,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:Nj),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:Nj,itlm),1,Yerr_tlm(1:N,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/(DirectionHgam)
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	KPP_REAL, INTENT(IN) :: A(LU_NONZERO)
907	INTEGER, INTENT(IN) :: Pivot(N)
908	!~~~> InOut variables
909	KPP_REAL, INTENT(INOUT) :: b(N)
910
911	CALL KppSolve( A, b )
912	!~~~> Note: for a full matrix use Lapack:
913	! NRHS = 1
914	! CALL DGETRS( 'N', N , NRHS, A, N, Pivot, b, N, INFO )
915
916	ISTATUS(Nsol) = ISTATUS(Nsol) + 1
917
918	END SUBROUTINE ros_Solve
919
920
921
922	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
923	SUBROUTINE Ros2
924	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
925	! --- AN L-STABLE METHOD, 2 stages, order 2
926	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
927
928	IMPLICIT NONE
929	DOUBLE PRECISION g
930
931	g = 1.0d0 + 1.0d0/SQRT(2.0d0)
932
933	rosMethod = RS2
934	!~~~> Name of the method
935	ros_Name = 'ROS-2'
936	!~~~> Number of stages
937	ros_S = 2
938
939	!~~~> The coefficient matrices A and C are strictly lower triangular.
940	! The lower triangular (subdiagonal) elements are stored in row-wise order:
941	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
942	! The general mapping formula is:
943	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
944	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
945
946	ros_A(1) = (1.d0)/g
947	ros_C(1) = (-2.d0)/g
948	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
949	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
950	ros_NewF(1) = .TRUE.
951	ros_NewF(2) = .TRUE.
952	!~~~> M_i = Coefficients for new step solution
953	ros_M(1)= (3.d0)/(2.d0*g)
954	ros_M(2)= (1.d0)/(2.d0*g)
955	! E_i = Coefficients for error estimator
956	ros_E(1) = 1.d0/(2.d0*g)
957	ros_E(2) = 1.d0/(2.d0*g)
958	!~~~> ros_ELO = estimator of local order - the minimum between the
959	! main and the embedded scheme orders plus one
960	ros_ELO = 2.0d0
961	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
962	ros_Alpha(1) = 0.0d0
963	ros_Alpha(2) = 1.0d0
964	!~~~> Gamma_i = \sum_j gamma_{i,j}
965	ros_Gamma(1) = g
966	ros_Gamma(2) =-g
967
968	END SUBROUTINE Ros2
969
970
971	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
972	SUBROUTINE Ros3
973	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
974	! --- AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations
975	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
976
977	IMPLICIT NONE
978
979	rosMethod = RS3
980	!~~~> Name of the method
981	ros_Name = 'ROS-3'
982	!~~~> Number of stages
983	ros_S = 3
984
985	!~~~> The coefficient matrices A and C are strictly lower triangular.
986	! The lower triangular (subdiagonal) elements are stored in row-wise order:
987	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
988	! The general mapping formula is:
989	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
990	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
991
992	ros_A(1)= 1.d0
993	ros_A(2)= 1.d0
994	ros_A(3)= 0.d0
995
996	ros_C(1) = -0.10156171083877702091975600115545d+01
997	ros_C(2) = 0.40759956452537699824805835358067d+01
998	ros_C(3) = 0.92076794298330791242156818474003d+01
999	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1000	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1001	ros_NewF(1) = .TRUE.
1002	ros_NewF(2) = .TRUE.
1003	ros_NewF(3) = .FALSE.
1004	!~~~> M_i = Coefficients for new step solution
1005	ros_M(1) = 0.1d+01
1006	ros_M(2) = 0.61697947043828245592553615689730d+01
1007	ros_M(3) = -0.42772256543218573326238373806514d+00
1008	! E_i = Coefficients for error estimator
1009	ros_E(1) = 0.5d+00
1010	ros_E(2) = -0.29079558716805469821718236208017d+01
1011	ros_E(3) = 0.22354069897811569627360909276199d+00
1012	!~~~> ros_ELO = estimator of local order - the minimum between the
1013	! main and the embedded scheme orders plus 1
1014	ros_ELO = 3.0d0
1015	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1016	ros_Alpha(1)= 0.0d+00
1017	ros_Alpha(2)= 0.43586652150845899941601945119356d+00
1018	ros_Alpha(3)= 0.43586652150845899941601945119356d+00
1019	!~~~> Gamma_i = \sum_j gamma_{i,j}
1020	ros_Gamma(1)= 0.43586652150845899941601945119356d+00
1021	ros_Gamma(2)= 0.24291996454816804366592249683314d+00
1022	ros_Gamma(3)= 0.21851380027664058511513169485832d+01
1023
1024	END SUBROUTINE Ros3
1025
1026	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1027
1028
1029	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1030	SUBROUTINE Ros4
1031	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1032	! L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES
1033	! L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
1034	!
1035	! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
1036	! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
1037	! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
1038	! SPRINGER-VERLAG (1990)
1039	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1040
1041	IMPLICIT NONE
1042
1043	rosMethod = RS4
1044	!~~~> Name of the method
1045	ros_Name = 'ROS-4'
1046	!~~~> Number of stages
1047	ros_S = 4
1048
1049	!~~~> The coefficient matrices A and C are strictly lower triangular.
1050	! The lower triangular (subdiagonal) elements are stored in row-wise order:
1051	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1052	! The general mapping formula is:
1053	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
1054	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
1055
1056	ros_A(1) = 0.2000000000000000d+01
1057	ros_A(2) = 0.1867943637803922d+01
1058	ros_A(3) = 0.2344449711399156d+00
1059	ros_A(4) = ros_A(2)
1060	ros_A(5) = ros_A(3)
1061	ros_A(6) = 0.0D0
1062
1063	ros_C(1) =-0.7137615036412310d+01
1064	ros_C(2) = 0.2580708087951457d+01
1065	ros_C(3) = 0.6515950076447975d+00
1066	ros_C(4) =-0.2137148994382534d+01
1067	ros_C(5) =-0.3214669691237626d+00
1068	ros_C(6) =-0.6949742501781779d+00
1069	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1070	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1071	ros_NewF(1) = .TRUE.
1072	ros_NewF(2) = .TRUE.
1073	ros_NewF(3) = .TRUE.
1074	ros_NewF(4) = .FALSE.
1075	!~~~> M_i = Coefficients for new step solution
1076	ros_M(1) = 0.2255570073418735d+01
1077	ros_M(2) = 0.2870493262186792d+00
1078	ros_M(3) = 0.4353179431840180d+00
1079	ros_M(4) = 0.1093502252409163d+01
1080	!~~~> E_i = Coefficients for error estimator
1081	ros_E(1) =-0.2815431932141155d+00
1082	ros_E(2) =-0.7276199124938920d-01
1083	ros_E(3) =-0.1082196201495311d+00
1084	ros_E(4) =-0.1093502252409163d+01
1085	!~~~> ros_ELO = estimator of local order - the minimum between the
1086	! main and the embedded scheme orders plus 1
1087	ros_ELO = 4.0d0
1088	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1089	ros_Alpha(1) = 0.D0
1090	ros_Alpha(2) = 0.1145640000000000d+01
1091	ros_Alpha(3) = 0.6552168638155900d+00
1092	ros_Alpha(4) = ros_Alpha(3)
1093	!~~~> Gamma_i = \sum_j gamma_{i,j}
1094	ros_Gamma(1) = 0.5728200000000000d+00
1095	ros_Gamma(2) =-0.1769193891319233d+01
1096	ros_Gamma(3) = 0.7592633437920482d+00
1097	ros_Gamma(4) =-0.1049021087100450d+00
1098
1099	END SUBROUTINE Ros4
1100
1101	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1102	SUBROUTINE Rodas3
1103	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1104	! --- A STIFFLY-STABLE METHOD, 4 stages, order 3
1105	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1106
1107	IMPLICIT NONE
1108
1109	rosMethod = RD3
1110	!~~~> Name of the method
1111	ros_Name = 'RODAS-3'
1112	!~~~> Number of stages
1113	ros_S = 4
1114
1115	!~~~> The coefficient matrices A and C are strictly lower triangular.
1116	! The lower triangular (subdiagonal) elements are stored in row-wise order:
1117	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1118	! The general mapping formula is:
1119	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
1120	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
1121
1122	ros_A(1) = 0.0d+00
1123	ros_A(2) = 2.0d+00
1124	ros_A(3) = 0.0d+00
1125	ros_A(4) = 2.0d+00
1126	ros_A(5) = 0.0d+00
1127	ros_A(6) = 1.0d+00
1128
1129	ros_C(1) = 4.0d+00
1130	ros_C(2) = 1.0d+00
1131	ros_C(3) =-1.0d+00
1132	ros_C(4) = 1.0d+00
1133	ros_C(5) =-1.0d+00
1134	ros_C(6) =-(8.0d+00/3.0d+00)
1135
1136	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1137	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1138	ros_NewF(1) = .TRUE.
1139	ros_NewF(2) = .FALSE.
1140	ros_NewF(3) = .TRUE.
1141	ros_NewF(4) = .TRUE.
1142	!~~~> M_i = Coefficients for new step solution
1143	ros_M(1) = 2.0d+00
1144	ros_M(2) = 0.0d+00
1145	ros_M(3) = 1.0d+00
1146	ros_M(4) = 1.0d+00
1147	!~~~> E_i = Coefficients for error estimator
1148	ros_E(1) = 0.0d+00
1149	ros_E(2) = 0.0d+00
1150	ros_E(3) = 0.0d+00
1151	ros_E(4) = 1.0d+00
1152	!~~~> ros_ELO = estimator of local order - the minimum between the
1153	! main and the embedded scheme orders plus 1
1154	ros_ELO = 3.0d+00
1155	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1156	ros_Alpha(1) = 0.0d+00
1157	ros_Alpha(2) = 0.0d+00
1158	ros_Alpha(3) = 1.0d+00
1159	ros_Alpha(4) = 1.0d+00
1160	!~~~> Gamma_i = \sum_j gamma_{i,j}
1161	ros_Gamma(1) = 0.5d+00
1162	ros_Gamma(2) = 1.5d+00
1163	ros_Gamma(3) = 0.0d+00
1164	ros_Gamma(4) = 0.0d+00
1165
1166	END SUBROUTINE Rodas3
1167
1168	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1169	SUBROUTINE Rodas4
1170	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1171	! STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES
1172	!
1173	! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
1174	! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
1175	! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
1176	! SPRINGER-VERLAG (1996)
1177	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1178
1179	IMPLICIT NONE
1180
1181	rosMethod = RD4
1182	!~~~> Name of the method
1183	ros_Name = 'RODAS-4'
1184	!~~~> Number of stages
1185	ros_S = 6
1186
1187	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1188	ros_Alpha(1) = 0.000d0
1189	ros_Alpha(2) = 0.386d0
1190	ros_Alpha(3) = 0.210d0
1191	ros_Alpha(4) = 0.630d0
1192	ros_Alpha(5) = 1.000d0
1193	ros_Alpha(6) = 1.000d0
1194
1195	!~~~> Gamma_i = \sum_j gamma_{i,j}
1196	ros_Gamma(1) = 0.2500000000000000d+00
1197	ros_Gamma(2) =-0.1043000000000000d+00
1198	ros_Gamma(3) = 0.1035000000000000d+00
1199	ros_Gamma(4) =-0.3620000000000023d-01
1200	ros_Gamma(5) = 0.0d0
1201	ros_Gamma(6) = 0.0d0
1202
1203	!~~~> The coefficient matrices A and C are strictly lower triangular.
1204	! The lower triangular (subdiagonal) elements are stored in row-wise order:
1205	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1206	! The general mapping formula is: A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
1207	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
1208
1209	ros_A(1) = 0.1544000000000000d+01
1210	ros_A(2) = 0.9466785280815826d+00
1211	ros_A(3) = 0.2557011698983284d+00
1212	ros_A(4) = 0.3314825187068521d+01
1213	ros_A(5) = 0.2896124015972201d+01
1214	ros_A(6) = 0.9986419139977817d+00
1215	ros_A(7) = 0.1221224509226641d+01
1216	ros_A(8) = 0.6019134481288629d+01
1217	ros_A(9) = 0.1253708332932087d+02
1218	ros_A(10) =-0.6878860361058950d+00
1219	ros_A(11) = ros_A(7)
1220	ros_A(12) = ros_A(8)
1221	ros_A(13) = ros_A(9)
1222	ros_A(14) = ros_A(10)
1223	ros_A(15) = 1.0d+00
1224
1225	ros_C(1) =-0.5668800000000000d+01
1226	ros_C(2) =-0.2430093356833875d+01
1227	ros_C(3) =-0.2063599157091915d+00
1228	ros_C(4) =-0.1073529058151375d+00
1229	ros_C(5) =-0.9594562251023355d+01
1230	ros_C(6) =-0.2047028614809616d+02
1231	ros_C(7) = 0.7496443313967647d+01
1232	ros_C(8) =-0.1024680431464352d+02
1233	ros_C(9) =-0.3399990352819905d+02
1234	ros_C(10) = 0.1170890893206160d+02
1235	ros_C(11) = 0.8083246795921522d+01
1236	ros_C(12) =-0.7981132988064893d+01
1237	ros_C(13) =-0.3152159432874371d+02
1238	ros_C(14) = 0.1631930543123136d+02
1239	ros_C(15) =-0.6058818238834054d+01
1240
1241	!~~~> M_i = Coefficients for new step solution
1242	ros_M(1) = ros_A(7)
1243	ros_M(2) = ros_A(8)
1244	ros_M(3) = ros_A(9)
1245	ros_M(4) = ros_A(10)
1246	ros_M(5) = 1.0d+00
1247	ros_M(6) = 1.0d+00
1248
1249	!~~~> E_i = Coefficients for error estimator
1250	ros_E(1) = 0.0d+00
1251	ros_E(2) = 0.0d+00
1252	ros_E(3) = 0.0d+00
1253	ros_E(4) = 0.0d+00
1254	ros_E(5) = 0.0d+00
1255	ros_E(6) = 1.0d+00
1256
1257	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1258	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1259	ros_NewF(1) = .TRUE.
1260	ros_NewF(2) = .TRUE.
1261	ros_NewF(3) = .TRUE.
1262	ros_NewF(4) = .TRUE.
1263	ros_NewF(5) = .TRUE.
1264	ros_NewF(6) = .TRUE.
1265
1266	!~~~> ros_ELO = estimator of local order - the minimum between the
1267	! main and the embedded scheme orders plus 1
1268	ros_ELO = 4.0d0
1269
1270	END SUBROUTINE Rodas4
1271
1272	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1273	END SUBROUTINE RosenbrockTLM
1274	! and all its internal procedures
1275	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1276
1277
1278	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1279	SUBROUTINE FunTemplate( T, Y, Ydot )
1280	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1281	! Template for the ODE function call.
1282	! Updates the rate coefficients (and possibly the fixed species) at each call
1283	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1284
1285	IMPLICIT NONE
1286	!~~~> Input variables
1287	KPP_REAL :: T, Y(NVAR)
1288	!~~~> Output variables
1289	KPP_REAL :: Ydot(NVAR)
1290	!~~~> Local variables
1291	KPP_REAL :: Told
1292
1293	Told = TIME
1294	TIME = T
1295	CALL Update_SUN()
1296	CALL Update_RCONST()
1297	CALL Fun( Y, FIX, RCONST, Ydot )
1298	TIME = Told
1299
1300	END SUBROUTINE FunTemplate
1301
1302
1303	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1304	SUBROUTINE JacTemplate( T, Y, Jcb )
1305	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1306	! Template for the ODE Jacobian call.
1307	! Updates the rate coefficients (and possibly the fixed species) at each call
1308	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1309	IMPLICIT NONE
1310
1311	!~~~> Input variables
1312	KPP_REAL :: T, Y(NVAR)
1313	!~~~> Output variables
1314	KPP_REAL :: Jcb(LU_NONZERO)
1315	!~~~> Local variables
1316	KPP_REAL :: Told
1317
1318	Told = TIME
1319	TIME = T
1320	CALL Update_SUN()
1321	CALL Update_RCONST()
1322	CALL Jac_SP( Y, FIX, RCONST, Jcb )
1323	TIME = Told
1324
1325	END SUBROUTINE JacTemplate
1326
1327
1328
1329	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1330	SUBROUTINE HessTemplate( T, Y, Hes )
1331	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1332	! Template for the ODE Hessian call.
1333	! Updates the rate coefficients (and possibly the fixed species) at each call
1334	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1335	IMPLICIT NONE
1336
1337	!~~~> Input variables
1338	KPP_REAL :: T, Y(NVAR)
1339	!~~~> Output variables
1340	KPP_REAL :: Hes(NHESS)
1341	!~~~> Local variables
1342	KPP_REAL :: Told
1343
1344	Told = TIME
1345	TIME = T
1346	CALL Update_SUN()
1347	CALL Update_RCONST()
1348	CALL Hessian( Y, FIX, RCONST, Hes )
1349	TIME = Told
1350
1351	END SUBROUTINE HessTemplate
1352
1353	END MODULE KPP_ROOT_Integrator
1354
1355
1356
1357

Note: See TracBrowser for help on using the repository browser.

Context Navigation

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

Download in other formats: