Home

Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

rosenbrock_tlm.f90 @ 4116

Last change on this file since 4116 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
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,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
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),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)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,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)/(DirectionH)
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 = DirectionHros_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/(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	#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	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1279	END SUBROUTINE RosenbrockTLM
1280	! and all its internal procedures
1281	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1282
1283
1284	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1285	SUBROUTINE 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
1306	END SUBROUTINE FunTemplate
1307
1308
1309	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1310	SUBROUTINE 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
1331	END SUBROUTINE JacTemplate
1332
1333
1334
1335	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1336	SUBROUTINE 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
1357	END SUBROUTINE HessTemplate
1358
1359	END MODULE KPP_ROOT_Integrator
1360
1361
1362
1363

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

Context Navigation

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

Download in other formats: