Home

Context Navigation

runge_kutta_tlm.f90

Last change on this file was 2696, checked in by kanani, 6 years ago
Merge of branch palm4u into trunk
File size: 79.9 KB

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

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

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int.modified_WCOPY/runge_kutta_tlm.f90

Download in other formats: