Home

Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

runge_kutta.f90 @ 4797

Last change on this file since 4797 was 2696, checked in by kanani, 7 years ago
Merge of branch palm4u into trunk
File size: 68.0 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	INTEGER :: ISING
1037	#else
1038	KPP_REAL :: E1(LU_NONZERO)
1039	COMPLEX(kind=dp) :: E2(LU_NONZERO)
1040	#endif
1041	KPP_REAL :: R1(N),R2(N),R3(N)
1042	KPP_REAL :: H, x1, x2, x3
1043	COMPLEX(kind=dp) :: BC(N)
1044	INTEGER :: i
1045	!
1046	! Z <- h^{-1) T^{-1) A^{-1) x Z
1047	DO i=1,N
1048	x1 = R1(i)/H; x2 = R2(i)/H; x3 = R3(i)/H
1049	R1(i) = rkTinvAinv(1,1)x1 + rkTinvAinv(1,2)x2 + rkTinvAinv(1,3)*x3
1050	R2(i) = rkTinvAinv(2,1)x1 + rkTinvAinv(2,2)x2 + rkTinvAinv(2,3)*x3
1051	R3(i) = rkTinvAinv(3,1)x1 + rkTinvAinv(3,2)x2 + rkTinvAinv(3,3)*x3
1052	END DO
1053
1054	#ifdef FULL_ALGEBRA
1055	CALL DGETRS ('N',N,1,E1,N,IP1,R1,N,ISING)
1056	#else
1057	CALL KppSolve (E1,R1)
1058	#endif
1059	!
1060	DO i=1,N
1061	BC(i) = DCMPLX(R2(i),R3(i))
1062	END DO
1063	#ifdef FULL_ALGEBRA
1064	CALL ZGETRS ('N',N,1,E2,N,IP2,BC,N,ISING)
1065	#else
1066	CALL KppSolveCmplx (E2,BC)
1067	#endif
1068	DO i=1,N
1069	R2(i) = DBLE( BC(i) )
1070	R3(i) = AIMAG( BC(i) )
1071	END DO
1072
1073	! Z <- T x Z
1074	DO i=1,N
1075	x1 = R1(i); x2 = R2(i); x3 = R3(i)
1076	R1(i) = rkT(1,1)x1 + rkT(1,2)x2 + rkT(1,3)*x3
1077	R2(i) = rkT(2,1)x1 + rkT(2,2)x2 + rkT(2,3)*x3
1078	R3(i) = rkT(3,1)x1 + rkT(3,2)x2 + rkT(3,3)*x3
1079	END DO
1080
1081	ISTATUS(Nsol) = ISTATUS(Nsol) + 1
1082
1083	END SUBROUTINE RK_Solve
1084
1085
1086	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1087	SUBROUTINE RK_ErrorEstimate(N,H,T,Y,F0, &
1088	E1,IP1,Z1,Z2,Z3,SCAL,Err, &
1089	FirstStep,Reject)
1090	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1091	IMPLICIT NONE
1092
1093	INTEGER :: N
1094	#ifdef FULL_ALGEBRA
1095	KPP_REAL :: E1(NVAR,NVAR)
1096	INTEGER :: ISING
1097	#else
1098	KPP_REAL :: E1(LU_NONZERO)
1099	#endif
1100	KPP_REAL :: SCAL(N),Z1(N),Z2(N),Z3(N),F1(N),F2(N), &
1101	F0(N),Y(N),TMP(N),T,H
1102	INTEGER :: IP1(N), i
1103	LOGICAL FirstStep,Reject
1104	KPP_REAL :: HrkE1,HrkE2,HrkE3,Err
1105
1106	HrkE1 = rkE(1)/H
1107	HrkE2 = rkE(2)/H
1108	HrkE3 = rkE(3)/H
1109
1110	DO i=1,N
1111	F2(i) = HrkE1Z1(i)+HrkE2Z2(i)+HrkE3*Z3(i)
1112	TMP(i) = rkE(0)*F0(i) + F2(i)
1113	END DO
1114
1115
1116	#ifdef FULL_ALGEBRA
1117	CALL DGETRS ('N',N,1,E1,N,IP1,TMP,N,ISING)
1118	IF ((rkMethod==R1A).OR.(rkMethod==GAU).OR.(rkMethod==L3A)) THEN
1119	CALL DGETRS ('N',N,1,E1,N,IP1,TMP,N,ISING)
1120	END IF
1121	IF (rkMethod==GAU) THEN
1122	CALL DGETRS ('N',N,1,E1,N,IP1,TMP,N,ISING)
1123	END IF
1124	#else
1125	CALL KppSolve (E1, TMP)
1126	IF ((rkMethod==R1A).OR.(rkMethod==GAU).OR.(rkMethod==L3A)) THEN
1127	CALL KppSolve (E1,TMP)
1128	END IF
1129	IF (rkMethod==GAU) THEN
1130	CALL KppSolve (E1,TMP)
1131	END IF
1132	#endif
1133
1134	Err = RK_ErrorNorm(N,SCAL,TMP)
1135	!
1136	IF (Err < ONE) RETURN
1137	firej:IF (FirstStep.OR.Reject) THEN
1138	DO i=1,N
1139	TMP(i)=Y(i)+TMP(i)
1140	END DO
1141	CALL FUN_CHEM(T,TMP,F1)
1142	ISTATUS(Nfun) = ISTATUS(Nfun) + 1
1143	DO i=1,N
1144	TMP(i)=F1(i)+F2(i)
1145	END DO
1146
1147	#ifdef FULL_ALGEBRA
1148	CALL DGETRS ('N',N,1,E1,N,IP1,TMP,N,ISING)
1149	#else
1150	CALL KppSolve (E1, TMP)
1151	#endif
1152	Err = RK_ErrorNorm(N,SCAL,TMP)
1153	END IF firej
1154
1155	END SUBROUTINE RK_ErrorEstimate
1156
1157
1158	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1159	KPP_REAL FUNCTION RK_ErrorNorm(N,SCAL,DY)
1160	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1161	IMPLICIT NONE
1162
1163	INTEGER :: N
1164	KPP_REAL :: SCAL(N),DY(N)
1165	INTEGER :: i
1166
1167	RK_ErrorNorm = ZERO
1168	DO i=1,N
1169	RK_ErrorNorm = RK_ErrorNorm + (DY(i)SCAL(i))*2
1170	END DO
1171	RK_ErrorNorm = MAX( SQRT(RK_ErrorNorm/N), 1.0d-10 )
1172
1173	END FUNCTION RK_ErrorNorm
1174
1175	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1176	SUBROUTINE Radau2A_Coefficients
1177	! The coefficients of the 3-stage Radau-2A method
1178	! (given to ~30 accurate digits)
1179	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1180	IMPLICIT NONE
1181	! The coefficients of the Radau2A method
1182	KPP_REAL :: b0
1183
1184	! b0 = 1.0d0
1185	IF (SdirkError) THEN
1186	b0 = 0.2d-1
1187	ELSE
1188	b0 = 0.5d-1
1189	END IF
1190
1191	! The coefficients of the Radau2A method
1192	rkMethod = R2A
1193
1194	rkA(1,1) = 1.968154772236604258683861429918299d-1
1195	rkA(1,2) = -6.55354258501983881085227825696087d-2
1196	rkA(1,3) = 2.377097434822015242040823210718965d-2
1197	rkA(2,1) = 3.944243147390872769974116714584975d-1
1198	rkA(2,2) = 2.920734116652284630205027458970589d-1
1199	rkA(2,3) = -4.154875212599793019818600988496743d-2
1200	rkA(3,1) = 3.764030627004672750500754423692808d-1
1201	rkA(3,2) = 5.124858261884216138388134465196080d-1
1202	rkA(3,3) = 1.111111111111111111111111111111111d-1
1203
1204	rkB(1) = 3.764030627004672750500754423692808d-1
1205	rkB(2) = 5.124858261884216138388134465196080d-1
1206	rkB(3) = 1.111111111111111111111111111111111d-1
1207
1208	rkC(1) = 1.550510257216821901802715925294109d-1
1209	rkC(2) = 6.449489742783178098197284074705891d-1
1210	rkC(3) = 1.0d0
1211
1212	! New solution: H* Sum B_jf(Z_j) = Sum D_jZ_j
1213	rkD(1) = 0.0d0
1214	rkD(2) = 0.0d0
1215	rkD(3) = 1.0d0
1216
1217	! Classical error estimator:
1218	! H* Sum (B_j-Bhat_j)f(Z_j) = HE(0)f(0) + Sum E_jZ_j
1219	rkE(0) = 1.0d0*b0
1220	rkE(1) = -10.04880939982741556246032950764708d0*b0
1221	rkE(2) = 1.382142733160748895793662840980412d0*b0
1222	rkE(3) = -.3333333333333333333333333333333333d0*b0
1223
1224	! Sdirk error estimator
1225	rkBgam(0) = b0
1226	rkBgam(1) = .3764030627004672750500754423692807d0-1.558078204724922382431975370686279d0*b0
1227	rkBgam(2) = .8914115380582557157653087040196118d0*b0+.5124858261884216138388134465196077d0
1228	rkBgam(3) = -.1637777184845662566367174924883037d0-.3333333333333333333333333333333333d0*b0
1229	rkBgam(4) = .2748888295956773677478286035994148d0
1230
1231	! H* Sum Bgam_jf(Z_j) = HBgam(0)f(0) + Sum Theta_jZ_j
1232	rkTheta(1) = -1.520677486405081647234271944611547d0-10.04880939982741556246032950764708d0*b0
1233	rkTheta(2) = 2.070455145596436382729929151810376d0+1.382142733160748895793662840980413d0*b0
1234	rkTheta(3) = -.3333333333333333333333333333333333d0*b0-.3744441479783868387391430179970741d0
1235
1236	! Local order of error estimator
1237	IF (b0==0.0d0) THEN
1238	rkELO = 6.0d0
1239	ELSE
1240	rkELO = 4.0d0
1241	END IF
1242
1243	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1244	!~~~> Diagonalize the RK matrix:
1245	! rkTinv * inv(rkA) * rkT =
1246	! \| rkGamma 0 0 \|
1247	! \| 0 rkAlpha -rkBeta \|
1248	! \| 0 rkBeta rkAlpha \|
1249	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1250
1251	rkGamma = 3.637834252744495732208418513577775d0
1252	rkAlpha = 2.681082873627752133895790743211112d0
1253	rkBeta = 3.050430199247410569426377624787569d0
1254
1255	rkT(1,1) = 9.443876248897524148749007950641664d-2
1256	rkT(1,2) = -1.412552950209542084279903838077973d-1
1257	rkT(1,3) = -3.00291941051474244918611170890539d-2
1258	rkT(2,1) = 2.502131229653333113765090675125018d-1
1259	rkT(2,2) = 2.041293522937999319959908102983381d-1
1260	rkT(2,3) = 3.829421127572619377954382335998733d-1
1261	rkT(3,1) = 1.0d0
1262	rkT(3,2) = 1.0d0
1263	rkT(3,3) = 0.0d0
1264
1265	rkTinv(1,1) = 4.178718591551904727346462658512057d0
1266	rkTinv(1,2) = 3.27682820761062387082533272429617d-1
1267	rkTinv(1,3) = 5.233764454994495480399309159089876d-1
1268	rkTinv(2,1) = -4.178718591551904727346462658512057d0
1269	rkTinv(2,2) = -3.27682820761062387082533272429617d-1
1270	rkTinv(2,3) = 4.766235545005504519600690840910124d-1
1271	rkTinv(3,1) = -5.02872634945786875951247343139544d-1
1272	rkTinv(3,2) = 2.571926949855605429186785353601676d0
1273	rkTinv(3,3) = -5.960392048282249249688219110993024d-1
1274
1275	rkTinvAinv(1,1) = 1.520148562492775501049204957366528d+1
1276	rkTinvAinv(1,2) = 1.192055789400527921212348994770778d0
1277	rkTinvAinv(1,3) = 1.903956760517560343018332287285119d0
1278	rkTinvAinv(2,1) = -9.669512977505946748632625374449567d0
1279	rkTinvAinv(2,2) = -8.724028436822336183071773193986487d0
1280	rkTinvAinv(2,3) = 3.096043239482439656981667712714881d0
1281	rkTinvAinv(3,1) = -1.409513259499574544876303981551774d+1
1282	rkTinvAinv(3,2) = 5.895975725255405108079130152868952d0
1283	rkTinvAinv(3,3) = -1.441236197545344702389881889085515d-1
1284
1285	rkAinvT(1,1) = .3435525649691961614912493915818282d0
1286	rkAinvT(1,2) = -.4703191128473198422370558694426832d0
1287	rkAinvT(1,3) = .3503786597113668965366406634269080d0
1288	rkAinvT(2,1) = .9102338692094599309122768354288852d0
1289	rkAinvT(2,2) = 1.715425895757991796035292755937326d0
1290	rkAinvT(2,3) = .4040171993145015239277111187301784d0
1291	rkAinvT(3,1) = 3.637834252744495732208418513577775d0
1292	rkAinvT(3,2) = 2.681082873627752133895790743211112d0
1293	rkAinvT(3,3) = -3.050430199247410569426377624787569d0
1294
1295	END SUBROUTINE Radau2A_Coefficients
1296
1297
1298
1299
1300	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1301	SUBROUTINE Lobatto3C_Coefficients
1302	! The coefficients of the 3-stage Lobatto-3C method
1303	! (given to ~30 accurate digits)
1304	! The parameter b0 can be chosen to tune the error estimator
1305	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1306	IMPLICIT NONE
1307	KPP_REAL :: b0
1308
1309	rkMethod = L3C
1310
1311	! b0 = 1.0d0
1312	IF (SdirkError) THEN
1313	b0 = 0.2d0
1314	ELSE
1315	b0 = 0.5d0
1316	END IF
1317	! The coefficients of the Lobatto3C method
1318
1319	rkA(1,1) = .1666666666666666666666666666666667d0
1320	rkA(1,2) = -.3333333333333333333333333333333333d0
1321	rkA(1,3) = .1666666666666666666666666666666667d0
1322	rkA(2,1) = .1666666666666666666666666666666667d0
1323	rkA(2,2) = .4166666666666666666666666666666667d0
1324	rkA(2,3) = -.8333333333333333333333333333333333d-1
1325	rkA(3,1) = .1666666666666666666666666666666667d0
1326	rkA(3,2) = .6666666666666666666666666666666667d0
1327	rkA(3,3) = .1666666666666666666666666666666667d0
1328
1329	rkB(1) = .1666666666666666666666666666666667d0
1330	rkB(2) = .6666666666666666666666666666666667d0
1331	rkB(3) = .1666666666666666666666666666666667d0
1332
1333	rkC(1) = 0.0d0
1334	rkC(2) = 0.5d0
1335	rkC(3) = 1.0d0
1336
1337	! Classical error estimator, embedded solution:
1338	rkBhat(0) = b0
1339	rkBhat(1) = .16666666666666666666666666666666667d0-b0
1340	rkBhat(2) = .66666666666666666666666666666666667d0
1341	rkBhat(3) = .16666666666666666666666666666666667d0
1342
1343	! New solution: h Sum_j b_j f(Z_j) = sum d_j Z_j
1344	rkD(1) = 0.0d0
1345	rkD(2) = 0.0d0
1346	rkD(3) = 1.0d0
1347
1348	! Classical error estimator:
1349	! H* Sum (B_j-Bhat_j)f(Z_j) = HE(0)f(0) + Sum E_jZ_j
1350	rkE(0) = .3808338772072650364017425226487022*b0
1351	rkE(1) = -1.142501631621795109205227567946107*b0
1352	rkE(2) = -1.523335508829060145606970090594809*b0
1353	rkE(3) = .3808338772072650364017425226487022*b0
1354
1355	! Sdirk error estimator
1356	rkBgam(0) = b0
1357	rkBgam(1) = .1666666666666666666666666666666667d0-1.d0*b0
1358	rkBgam(2) = .6666666666666666666666666666666667d0
1359	rkBgam(3) = -.2141672105405983697350758559820354d0
1360	rkBgam(4) = .3808338772072650364017425226487021d0
1361
1362	! H* Sum Bgam_jf(Z_j) = HBgam(0)f(0) + Sum Theta_jZ_j
1363	rkTheta(1) = -3.d0*b0-.3808338772072650364017425226487021d0
1364	rkTheta(2) = -4.d0*b0+1.523335508829060145606970090594808d0
1365	rkTheta(3) = -.142501631621795109205227567946106d0+b0
1366
1367	! Local order of error estimator
1368	IF (b0==0.0d0) THEN
1369	rkELO = 5.0d0
1370	ELSE
1371	rkELO = 4.0d0
1372	END IF
1373
1374	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1375	!~~~> Diagonalize the RK matrix:
1376	! rkTinv * inv(rkA) * rkT =
1377	! \| rkGamma 0 0 \|
1378	! \| 0 rkAlpha -rkBeta \|
1379	! \| 0 rkBeta rkAlpha \|
1380	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1381
1382	rkGamma = 2.625816818958466716011888933765284d0
1383	rkAlpha = 1.687091590520766641994055533117359d0
1384	rkBeta = 2.508731754924880510838743672432351d0
1385
1386	rkT(1,1) = 1.d0
1387	rkT(1,2) = 1.d0
1388	rkT(1,3) = 0.d0
1389	rkT(2,1) = .4554100411010284672111720348287483d0
1390	rkT(2,2) = -.6027050205505142336055860174143743d0
1391	rkT(2,3) = -.4309321229203225731070721341350346d0
1392	rkT(3,1) = 2.195823345445647152832799205549709d0
1393	rkT(3,2) = -1.097911672722823576416399602774855d0
1394	rkT(3,3) = .7850032632435902184104551358922130d0
1395
1396	rkTinv(1,1) = .4205559181381766909344950150991349d0
1397	rkTinv(1,2) = .3488903392193734304046467270632057d0
1398	rkTinv(1,3) = .1915253879645878102698098373933487d0
1399	rkTinv(2,1) = .5794440818618233090655049849008650d0
1400	rkTinv(2,2) = -.3488903392193734304046467270632057d0
1401	rkTinv(2,3) = -.1915253879645878102698098373933487d0
1402	rkTinv(3,1) = -.3659705575742745254721332009249516d0
1403	rkTinv(3,2) = -1.463882230297098101888532803699806d0
1404	rkTinv(3,3) = .4702733607340189781407813565524989d0
1405
1406	rkTinvAinv(1,1) = 1.104302803159744452668648155627548d0
1407	rkTinvAinv(1,2) = .916122120694355522658740710823143d0
1408	rkTinvAinv(1,3) = .5029105849749601702795812241441172d0
1409	rkTinvAinv(2,1) = 1.895697196840255547331351844372453d0
1410	rkTinvAinv(2,2) = 3.083877879305644477341259289176857d0
1411	rkTinvAinv(2,3) = -1.502910584974960170279581224144117d0
1412	rkTinvAinv(3,1) = .8362439183082935036129145574774502d0
1413	rkTinvAinv(3,2) = -3.344975673233174014451658229909802d0
1414	rkTinvAinv(3,3) = .312908409479233358005944466882642d0
1415
1416	rkAinvT(1,1) = 2.625816818958466716011888933765282d0
1417	rkAinvT(1,2) = 1.687091590520766641994055533117358d0
1418	rkAinvT(1,3) = -2.508731754924880510838743672432351d0
1419	rkAinvT(2,1) = 1.195823345445647152832799205549710d0
1420	rkAinvT(2,2) = -2.097911672722823576416399602774855d0
1421	rkAinvT(2,3) = .7850032632435902184104551358922130d0
1422	rkAinvT(3,1) = 5.765829871932827589653709477334136d0
1423	rkAinvT(3,2) = .1170850640335862051731452613329320d0
1424	rkAinvT(3,3) = 4.078738281412060947659653944216779d0
1425
1426	END SUBROUTINE Lobatto3C_Coefficients
1427
1428	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1429	SUBROUTINE Gauss_Coefficients
1430	! The coefficients of the 3-stage Gauss method
1431	! (given to ~30 accurate digits)
1432	! The parameter b3 can be chosen by the user
1433	! to tune the error estimator
1434	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1435	IMPLICIT NONE
1436	KPP_REAL :: b0
1437	! The coefficients of the Gauss method
1438
1439
1440	rkMethod = GAU
1441
1442	! b0 = 4.0d0
1443	b0 = 0.1d0
1444
1445	! The coefficients of the Gauss method
1446
1447	rkA(1,1) = .1388888888888888888888888888888889d0
1448	rkA(1,2) = -.359766675249389034563954710966045d-1
1449	rkA(1,3) = .97894440153083260495800422294756d-2
1450	rkA(2,1) = .3002631949808645924380249472131556d0
1451	rkA(2,2) = .2222222222222222222222222222222222d0
1452	rkA(2,3) = -.224854172030868146602471694353778d-1
1453	rkA(3,1) = .2679883337624694517281977355483022d0
1454	rkA(3,2) = .4804211119693833479008399155410489d0
1455	rkA(3,3) = .1388888888888888888888888888888889d0
1456
1457	rkB(1) = .2777777777777777777777777777777778d0
1458	rkB(2) = .4444444444444444444444444444444444d0
1459	rkB(3) = .2777777777777777777777777777777778d0
1460
1461	rkC(1) = .1127016653792583114820734600217600d0
1462	rkC(2) = .5000000000000000000000000000000000d0
1463	rkC(3) = .8872983346207416885179265399782400d0
1464
1465	! Classical error estimator, embedded solution:
1466	rkBhat(0) = b0
1467	rkBhat(1) =-1.4788305577012361475298775666303999d0*b0 &
1468	+.27777777777777777777777777777777778d0
1469	rkBhat(2) = .44444444444444444444444444444444444d0 &
1470	+.66666666666666666666666666666666667d0*b0
1471	rkBhat(3) = -.18783610896543051913678910003626672d0*b0 &
1472	+.27777777777777777777777777777777778d0
1473
1474	! New solution: h Sum_j b_j f(Z_j) = sum d_j Z_j
1475	rkD(1) = .1666666666666666666666666666666667d1
1476	rkD(2) = -.1333333333333333333333333333333333d1
1477	rkD(3) = .1666666666666666666666666666666667d1
1478
1479	! Classical error estimator:
1480	! H* Sum (B_j-Bhat_j)f(Z_j) = HE(0)f(0) + Sum E_jZ_j
1481	rkE(0) = .2153144231161121782447335303806954d0*b0
1482	rkE(1) = -2.825278112319014084275808340593191d0*b0
1483	rkE(2) = .2870858974881495709929780405075939d0*b0
1484	rkE(3) = -.4558086256248162565397206448274867d-1*b0
1485
1486	! Sdirk error estimator
1487	rkBgam(0) = 0.d0
1488	rkBgam(1) = .2373339543355109188382583162660537d0
1489	rkBgam(2) = .5879873931885192299409334646982414d0
1490	rkBgam(3) = -.4063577064014232702392531134499046d-1
1491	rkBgam(4) = .2153144231161121782447335303806955d0
1492
1493	! H* Sum Bgam_jf(Z_j) = HBgam(0)f(0) + Sum Theta_jZ_j
1494	rkTheta(1) = -2.594040933093095272574031876464493d0
1495	rkTheta(2) = 1.824611539036311947589425112250199d0
1496	rkTheta(3) = .1856563166634371860478043996459493d0
1497
1498	! ELO = local order of classical error estimator
1499	rkELO = 4.0d0
1500
1501	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1502	!~~~> Diagonalize the RK matrix:
1503	! rkTinv * inv(rkA) * rkT =
1504	! \| rkGamma 0 0 \|
1505	! \| 0 rkAlpha -rkBeta \|
1506	! \| 0 rkBeta rkAlpha \|
1507	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1508
1509	rkGamma = 4.644370709252171185822941421408064d0
1510	rkAlpha = 3.677814645373914407088529289295970d0
1511	rkBeta = 3.508761919567443321903661209182446d0
1512
1513	rkT(1,1) = .7215185205520017032081769924397664d-1
1514	rkT(1,2) = -.8224123057363067064866206597516454d-1
1515	rkT(1,3) = -.6012073861930850173085948921439054d-1
1516	rkT(2,1) = .1188325787412778070708888193730294d0
1517	rkT(2,2) = .5306509074206139504614411373957448d-1
1518	rkT(2,3) = .3162050511322915732224862926182701d0
1519	rkT(3,1) = 1.d0
1520	rkT(3,2) = 1.d0
1521	rkT(3,3) = 0.d0
1522
1523	rkTinv(1,1) = 5.991698084937800775649580743981285d0
1524	rkTinv(1,2) = 1.139214295155735444567002236934009d0
1525	rkTinv(1,3) = .4323121137838583855696375901180497d0
1526	rkTinv(2,1) = -5.991698084937800775649580743981285d0
1527	rkTinv(2,2) = -1.139214295155735444567002236934009d0
1528	rkTinv(2,3) = .5676878862161416144303624098819503d0
1529	rkTinv(3,1) = -1.246213273586231410815571640493082d0
1530	rkTinv(3,2) = 2.925559646192313662599230367054972d0
1531	rkTinv(3,3) = -.2577352012734324923468722836888244d0
1532
1533	rkTinvAinv(1,1) = 27.82766708436744962047620566703329d0
1534	rkTinvAinv(1,2) = 5.290933503982655311815946575100597d0
1535	rkTinvAinv(1,3) = 2.007817718512643701322151051660114d0
1536	rkTinvAinv(2,1) = -17.66368928942422710690385180065675d0
1537	rkTinvAinv(2,2) = -14.45491129892587782538830044147713d0
1538	rkTinvAinv(2,3) = 2.992182281487356298677848948339886d0
1539	rkTinvAinv(3,1) = -25.60678350282974256072419392007303d0
1540	rkTinvAinv(3,2) = 6.762434375611708328910623303779923d0
1541	rkTinvAinv(3,3) = 1.043979339483109825041215970036771d0
1542
1543	rkAinvT(1,1) = .3350999483034677402618981153470483d0
1544	rkAinvT(1,2) = -.5134173605009692329246186488441294d0
1545	rkAinvT(1,3) = .6745196507033116204327635673208923d-1
1546	rkAinvT(2,1) = .5519025480108928886873752035738885d0
1547	rkAinvT(2,2) = 1.304651810077110066076640761092008d0
1548	rkAinvT(2,3) = .9767507983414134987545585703726984d0
1549	rkAinvT(3,1) = 4.644370709252171185822941421408064d0
1550	rkAinvT(3,2) = 3.677814645373914407088529289295970d0
1551	rkAinvT(3,3) = -3.508761919567443321903661209182446d0
1552
1553	END SUBROUTINE Gauss_Coefficients
1554
1555
1556
1557	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1558	SUBROUTINE Radau1A_Coefficients
1559	! The coefficients of the 3-stage Gauss method
1560	! (given to ~30 accurate digits)
1561	! The parameter b3 can be chosen by the user
1562	! to tune the error estimator
1563	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1564	IMPLICIT NONE
1565	! KPP_REAL :: b0 = 0.3d0
1566	KPP_REAL :: b0 = 0.1d0
1567
1568	! The coefficients of the Radau1A method
1569
1570	rkMethod = R1A
1571
1572	rkA(1,1) = .1111111111111111111111111111111111d0
1573	rkA(1,2) = -.1916383190435098943442935597058829d0
1574	rkA(1,3) = .8052720793239878323318244859477174d-1
1575	rkA(2,1) = .1111111111111111111111111111111111d0
1576	rkA(2,2) = .2920734116652284630205027458970589d0
1577	rkA(2,3) = -.481334970546573839513422644787591d-1
1578	rkA(3,1) = .1111111111111111111111111111111111d0
1579	rkA(3,2) = .5370223859435462728402311533676479d0
1580	rkA(3,3) = .1968154772236604258683861429918299d0
1581
1582	rkB(1) = .1111111111111111111111111111111111d0
1583	rkB(2) = .5124858261884216138388134465196080d0
1584	rkB(3) = .3764030627004672750500754423692808d0
1585
1586	rkC(1) = 0.d0
1587	rkC(2) = .3550510257216821901802715925294109d0
1588	rkC(3) = .8449489742783178098197284074705891d0
1589
1590	! Classical error estimator, embedded solution:
1591	rkBhat(0) = b0
1592	rkBhat(1) = .11111111111111111111111111111111111d0-b0
1593	rkBhat(2) = .51248582618842161383881344651960810d0
1594	rkBhat(3) = .37640306270046727505007544236928079d0
1595
1596	! New solution: H* Sum B_jf(Z_j) = Sum D_jZ_j
1597	rkD(1) = .3333333333333333333333333333333333d0
1598	rkD(2) = -.8914115380582557157653087040196127d0
1599	rkD(3) = .1558078204724922382431975370686279d1
1600
1601	! Classical error estimator:
1602	! H* Sum (b_j-bhat_j) f(Z_j) = HE(0)F(0) + Sum E_j Z_j
1603	rkE(0) = .2748888295956773677478286035994148d0*b0
1604	rkE(1) = -1.374444147978386838739143017997074d0*b0
1605	rkE(2) = -1.335337922441686804550326197041126d0*b0
1606	rkE(3) = .235782604058977333559011782643466d0*b0
1607
1608	! Sdirk error estimator
1609	rkBgam(0) = 0.0d0
1610	rkBgam(1) = .1948150124588532186183490991130616d-1
1611	rkBgam(2) = .7575249005733381398986810981093584d0
1612	rkBgam(3) = -.518952314149008295083446116200793d-1
1613	rkBgam(4) = .2748888295956773677478286035994148d0
1614
1615	! H* Sum Bgam_jf(Z_j) = HBgam(0)f(0) + Sum Theta_jZ_j
1616	rkTheta(1) = -1.224370034375505083904362087063351d0
1617	rkTheta(2) = .9340045331532641409047527962010133d0
1618	rkTheta(3) = .4656990124352088397561234800640929d0
1619
1620	! ELO = local order of classical error estimator
1621	rkELO = 4.0d0
1622
1623	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1624	!~~~> Diagonalize the RK matrix:
1625	! rkTinv * inv(rkA) * rkT =
1626	! \| rkGamma 0 0 \|
1627	! \| 0 rkAlpha -rkBeta \|
1628	! \| 0 rkBeta rkAlpha \|
1629	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1630
1631	rkGamma = 3.637834252744495732208418513577775d0
1632	rkAlpha = 2.681082873627752133895790743211112d0
1633	rkBeta = 3.050430199247410569426377624787569d0
1634
1635	rkT(1,1) = .424293819848497965354371036408369d0
1636	rkT(1,2) = -.3235571519651980681202894497035503d0
1637	rkT(1,3) = -.522137786846287839586599927945048d0
1638	rkT(2,1) = .57594609499806128896291585429339d-1
1639	rkT(2,2) = .3148663231849760131614374283783d-2
1640	rkT(2,3) = .452429247674359778577728510381731d0
1641	rkT(3,1) = 1.d0
1642	rkT(3,2) = 1.d0
1643	rkT(3,3) = 0.d0
1644
1645	rkTinv(1,1) = 1.233523612685027760114769983066164d0
1646	rkTinv(1,2) = 1.423580134265707095505388133369554d0
1647	rkTinv(1,3) = .3946330125758354736049045150429624d0
1648	rkTinv(2,1) = -1.233523612685027760114769983066164d0
1649	rkTinv(2,2) = -1.423580134265707095505388133369554d0
1650	rkTinv(2,3) = .6053669874241645263950954849570376d0
1651	rkTinv(3,1) = -.1484438963257383124456490049673414d0
1652	rkTinv(3,2) = 2.038974794939896109682070471785315d0
1653	rkTinv(3,3) = -.544501292892686735299355831692542d-1
1654
1655	rkTinvAinv(1,1) = 4.487354449794728738538663081025420d0
1656	rkTinvAinv(1,2) = 5.178748573958397475446442544234494d0
1657	rkTinvAinv(1,3) = 1.435609490412123627047824222335563d0
1658	rkTinvAinv(2,1) = -2.854361287939276673073807031221493d0
1659	rkTinvAinv(2,2) = -1.003648660720543859000994063139137d+1
1660	rkTinvAinv(2,3) = 1.789135380979465422050817815017383d0
1661	rkTinvAinv(3,1) = -4.160768067752685525282947313530352d0
1662	rkTinvAinv(3,2) = 1.124128569859216916690209918405860d0
1663	rkTinvAinv(3,3) = 1.700644430961823796581896350418417d0
1664
1665	rkAinvT(1,1) = 1.543510591072668287198054583233180d0
1666	rkAinvT(1,2) = -2.460228411937788329157493833295004d0
1667	rkAinvT(1,3) = -.412906170450356277003910443520499d0
1668	rkAinvT(2,1) = .209519643211838264029272585946993d0
1669	rkAinvT(2,2) = 1.388545667194387164417459732995766d0
1670	rkAinvT(2,3) = 1.20339553005832004974976023130002d0
1671	rkAinvT(3,1) = 3.637834252744495732208418513577775d0
1672	rkAinvT(3,2) = 2.681082873627752133895790743211112d0
1673	rkAinvT(3,3) = -3.050430199247410569426377624787569d0
1674
1675	END SUBROUTINE Radau1A_Coefficients
1676
1677
1678	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1679	SUBROUTINE Lobatto3A_Coefficients
1680	! The coefficients of the 4-stage Lobatto-3A method
1681	! (given to ~30 accurate digits)
1682	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1683	IMPLICIT NONE
1684
1685	! The coefficients of the Lobatto-3A method
1686
1687	rkMethod = L3A
1688
1689	rkA(0,0) = 0.0d0
1690	rkA(0,1) = 0.0d0
1691	rkA(0,2) = 0.0d0
1692	rkA(0,3) = 0.0d0
1693	rkA(1,0) = .11030056647916491413674311390609397d0
1694	rkA(1,1) = .1896994335208350858632568860939060d0
1695	rkA(1,2) = -.339073642291438837776604807792215d-1
1696	rkA(1,3) = .1030056647916491413674311390609397d-1
1697	rkA(2,0) = .73032766854168419196590219427239365d-1
1698	rkA(2,1) = .4505740308958105504443271474458881d0
1699	rkA(2,2) = .2269672331458315808034097805727606d0
1700	rkA(2,3) = -.2696723314583158080340978057276063d-1
1701	rkA(3,0) = .83333333333333333333333333333333333d-1
1702	rkA(3,1) = .4166666666666666666666666666666667d0
1703	rkA(3,2) = .4166666666666666666666666666666667d0
1704	rkA(3,3) = .8333333333333333333333333333333333d-1
1705
1706	rkB(0) = .83333333333333333333333333333333333d-1
1707	rkB(1) = .4166666666666666666666666666666667d0
1708	rkB(2) = .4166666666666666666666666666666667d0
1709	rkB(3) = .8333333333333333333333333333333333d-1
1710
1711	rkC(0) = 0.0d0
1712	rkC(1) = .2763932022500210303590826331268724d0
1713	rkC(2) = .7236067977499789696409173668731276d0
1714	rkC(3) = 1.0d0
1715
1716	! New solution: HSum B_jf(Z_j) = Sum D_j*Z_j
1717	rkD(0) = 0.0d0
1718	rkD(1) = 0.0d0
1719	rkD(2) = 0.0d0
1720	rkD(3) = 1.0d0
1721
1722	! Classical error estimator, embedded solution:
1723	rkBhat(0) = .90909090909090909090909090909090909d-1
1724	rkBhat(1) = .39972675774621371442114262372173276d0
1725	rkBhat(2) = .43360657558711961891219070961160058d0
1726	rkBhat(3) = .15151515151515151515151515151515152d-1
1727
1728	! Classical error estimator:
1729	! H* Sum (B_j-Bhat_j)f(Z_j) = HE(0)f(0) + Sum E_jZ_j
1730
1731	rkE(0) = .1957403846510110711315759367097231d-1
1732	rkE(1) = -.1986820345632580910316020806676438d0
1733	rkE(2) = .1660586371214229125096727578826900d0
1734	rkE(3) = -.9787019232550553556578796835486154d-1
1735
1736	! Sdirk error estimator:
1737	rkF(0) = 0.0d0
1738	rkF(1) = -.66535815876916686607437314126436349d0
1739	rkF(2) = 1.7419302743497277572980407931678409d0
1740	rkF(3) = -1.2918865386966730694684011822841728d0
1741
1742	! ELO = local order of classical error estimator
1743	rkELO = 4.0d0
1744
1745	! Sdirk error estimator:
1746	rkBgam(0) = .2950472755430528877214995073815946d-1
1747	rkBgam(1) = .5370310883226113978352873633882769d0
1748	rkBgam(2) = .2963022450107219354980459699450564d0
1749	rkBgam(3) = -.7815248400375080035021681445218837d-1
1750	rkBgam(4) = .2153144231161121782447335303806956d0
1751
1752	! H* Sum Bgam_jf(Z_j) = HBgam(0)f(0) + Sum Theta_jZ_j
1753	rkTheta(0) = 0.0d0
1754	rkTheta(1) = -.6653581587691668660743731412643631d0
1755	rkTheta(2) = 1.741930274349727757298040793167842d0
1756	rkTheta(3) = -.291886538696673069468401182284174d0
1757
1758
1759	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1760	!~~~> Diagonalize the RK matrix:
1761	! rkTinv * inv(rkA) * rkT =
1762	! \| rkGamma 0 0 \|
1763	! \| 0 rkAlpha -rkBeta \|
1764	! \| 0 rkBeta rkAlpha \|
1765	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1766
1767	rkGamma = 4.644370709252171185822941421408063d0
1768	rkAlpha = 3.677814645373914407088529289295968d0
1769	rkBeta = 3.508761919567443321903661209182446d0
1770
1771	rkT(1,1) = .5303036326129938105898786144870856d-1
1772	rkT(1,2) = -.7776129960563076320631956091016914d-1
1773	rkT(1,3) = .6043307469475508514468017399717112d-2
1774	rkT(2,1) = .2637242522173698467283726114649606d0
1775	rkT(2,2) = .2193839918662961493126393244533346d0
1776	rkT(2,3) = .3198765142300936188514264752235344d0
1777	rkT(3,1) = 1.d0
1778	rkT(3,2) = 1.d0
1779	rkT(3,3) = 0.d0
1780
1781	rkTinv(1,1) = 7.695032983257654470769069079238553d0
1782	rkTinv(1,2) = -.1453793830957233720334601186354032d0
1783	rkTinv(1,3) = .6302696746849084900422461036874826d0
1784	rkTinv(2,1) = -7.695032983257654470769069079238553d0
1785	rkTinv(2,2) = .1453793830957233720334601186354032d0
1786	rkTinv(2,3) = .3697303253150915099577538963125174d0
1787	rkTinv(3,1) = -1.066660885401270392058552736086173d0
1788	rkTinv(3,2) = 3.146358406832537460764521760668932d0
1789	rkTinv(3,3) = -.7732056038202974770406168510664738d0
1790
1791	rkTinvAinv(1,1) = 35.73858579417120341641749040405149d0
1792	rkTinvAinv(1,2) = -.675195748578927863668368190236025d0
1793	rkTinvAinv(1,3) = 2.927206016036483646751158874041632d0
1794	rkTinvAinv(2,1) = -24.55824590667225493437162206039511d0
1795	rkTinvAinv(2,2) = -10.50514413892002061837750015342036
1796	rkTinvAinv(2,3) = 4.072793983963516353248841125958369d0
1797	rkTinvAinv(3,1) = -30.92301972744621647251975054630589d0
1798	rkTinvAinv(3,2) = 12.08182467154052413351908559269928d0
1799	rkTinvAinv(3,3) = -1.546411207640594954081233702132946d0
1800
1801	rkAinvT(1,1) = .2462926658317812882584158369803835d0
1802	rkAinvT(1,2) = -.2647871194157644619747121197289574d0
1803	rkAinvT(1,3) = .2950720515900466654896406799284586d0
1804	rkAinvT(2,1) = 1.224833192317784474576995878738004d0
1805	rkAinvT(2,2) = 1.929224190340981580557006261869763d0
1806	rkAinvT(2,3) = .4066803323234419988910915619080306d0
1807	rkAinvT(3,1) = 4.644370709252171185822941421408064d0
1808	rkAinvT(3,2) = 3.677814645373914407088529289295968d0
1809	rkAinvT(3,3) = -3.508761919567443321903661209182446d0
1810
1811	END SUBROUTINE Lobatto3A_Coefficients
1812
1813	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1814	END SUBROUTINE RungeKutta ! and all its internal procedures
1815	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1816
1817	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1818	SUBROUTINE FUN_CHEM(T, V, FCT)
1819	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1820
1821	USE KPP_ROOT_Parameters
1822	USE KPP_ROOT_Global
1823	USE KPP_ROOT_Function, ONLY: Fun
1824	USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO
1825
1826	IMPLICIT NONE
1827
1828	KPP_REAL :: V(NVAR), FCT(NVAR)
1829	KPP_REAL :: T, Told
1830
1831	Told = TIME
1832	TIME = T
1833	CALL Update_SUN()
1834	CALL Update_RCONST()
1835	CALL Update_PHOTO()
1836	TIME = Told
1837
1838	CALL Fun(V, FIX, RCONST, FCT)
1839
1840	END SUBROUTINE FUN_CHEM
1841
1842
1843	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1844	SUBROUTINE JAC_CHEM (T, V, JF)
1845	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1846
1847	USE KPP_ROOT_Parameters
1848	USE KPP_ROOT_Global
1849	USE KPP_ROOT_JacobianSP
1850	USE KPP_ROOT_Jacobian, ONLY: Jac_SP
1851	USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO
1852
1853	IMPLICIT NONE
1854
1855	KPP_REAL :: V(NVAR), T , Told
1856	#ifdef FULL_ALGEBRA
1857	KPP_REAL :: JV(LU_NONZERO), JF(NVAR,NVAR)
1858	INTEGER :: i, j
1859	#else
1860	KPP_REAL :: JF(LU_NONZERO)
1861	#endif
1862
1863	Told = TIME
1864	TIME = T
1865	CALL Update_SUN()
1866	CALL Update_RCONST()
1867	CALL Update_PHOTO()
1868	TIME = Told
1869
1870	#ifdef FULL_ALGEBRA
1871	CALL Jac_SP(V, FIX, RCONST, JV)
1872	DO j=1,NVAR
1873	DO i=1,NVAR
1874	JF(i,j) = 0.0d0
1875	END DO
1876	END DO
1877	DO i=1,LU_NONZERO
1878	JF(LU_IROW(i),LU_ICOL(i)) = JV(i)
1879	END DO
1880	#else
1881	CALL Jac_SP(V, FIX, RCONST, JF)
1882	#endif
1883
1884	END SUBROUTINE JAC_CHEM
1885
1886
1887	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1888
1889	END MODULE KPP_ROOT_Integrator
1890
1891

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

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/runge_kutta.f90 @ 4797

Download in other formats: