Home

Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

runge_kutta.f90 @ 3534

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

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

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.f90 @ 3534

Download in other formats: