Home

Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sdirk_tlm.f90 @ 2968

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

Line
1	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
2	! SDIRK-TLM - Tangent Linear Model of Singly-Diagonally-Implicit RK !
3	! * Sdirk 2a, 2b: L-stable, 2 stages, order 2 !
4	! * Sdirk 3a: L-stable, 3 stages, order 2, adj-invariant !
5	! * Sdirk 4a, 4b: L-stable, 5 stages, order 4 !
6	! By default the code employs the KPP sparse linear algebra routines !
7	! Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK) !
8	! !
9	! (C) Adrian Sandu, July 2005 !
10	! Virginia Polytechnic Institute and State University !
11	! Contact: sandu@cs.vt.edu !
12	! Revised by Philipp Miehe and Adrian Sandu, May 2006 !
13	! This implementation is part of KPP - the Kinetic PreProcessor !
14	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
15
16	MODULE KPP_ROOT_Integrator
17
18	USE KPP_ROOT_Precision
19	USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME
20	USE KPP_ROOT_Parameters, ONLY: NVAR, NSPEC, NFIX, LU_NONZERO
21	USE KPP_ROOT_JacobianSP, ONLY: LU_DIAG
22	USE KPP_ROOT_Jacobian, ONLY: Jac_SP_Vec
23	USE KPP_ROOT_LinearAlgebra, ONLY: KppDecomp, &
24	KppSolve, Set2zero, WLAMCH, WCOPY, WAXPY, WSCAL, WADD
25
26	IMPLICIT NONE
27	PUBLIC
28	SAVE
29
30	!~~~> Statistics on the work performed by the SDIRK method
31	INTEGER, PARAMETER :: Nfun=1, Njac=2, Nstp=3, Nacc=4, &
32	Nrej=5, Ndec=6, Nsol=7, Nsng=8, &
33	Ntexit=1, Nhexit=2, Nhnew=3
34
35	CONTAINS
36
37	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
38	SUBROUTINE INTEGRATE_TLM( NTLM, Y, Y_tlm, TIN, TOUT, ATOL_tlm,RTOL_tlm, &
39	ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U, IERR_U )
40	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
41
42	USE KPP_ROOT_Parameters, ONLY: NVAR,ind_O3
43	USE KPP_ROOT_Global, ONLY: ATOL,RTOL,VAR
44
45	IMPLICIT NONE
46
47	!~~~> Y - Concentrations
48	KPP_REAL :: Y(NVAR)
49	!~~~> NTLM - No. of sensitivity coefficients
50	INTEGER NTLM
51	!~~~> Y_tlm - Sensitivities of concentrations
52	! Note: Y_tlm (1:NVAR,j) contains sensitivities of
53	! Y(1:NVAR) w.r.t. the j-th parameter, j=1...NTLM
54	KPP_REAL :: Y_tlm(NVAR,NTLM)
55	KPP_REAL :: TIN ! TIN - Start Time
56	KPP_REAL :: TOUT ! TOUT - End Time
57	KPP_REAL, INTENT(IN), OPTIONAL :: RTOL_tlm(NVAR,NTLM),ATOL_tlm(NVAR,NTLM)
58	INTEGER, INTENT(IN), OPTIONAL :: ICNTRL_U(20)
59	KPP_REAL, INTENT(IN), OPTIONAL :: RCNTRL_U(20)
60	INTEGER, INTENT(OUT), OPTIONAL :: ISTATUS_U(20)
61	KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20)
62	INTEGER, INTENT(OUT), OPTIONAL :: IERR_U
63
64	INTEGER :: IERR
65
66	KPP_REAL :: RCNTRL(20), RSTATUS(20), T1, T2
67	INTEGER :: ICNTRL(20), ISTATUS(20)
68	INTEGER, SAVE :: Ntotal = 0
69
70	ICNTRL(1:20) = 0
71	RCNTRL(1:20) = 0.0_dp
72
73	!~~~> fine-tune the integrator:
74	ICNTRL(2) = 0 ! 0=vector tolerances, 1=scalar tolerances
75	ICNTRL(5) = 8 ! Max no. of Newton iterations
76	ICNTRL(6) = 0 ! Starting values for Newton are interpolated (0) or zero (1)
77	ICNTRL(7) = 0 ! How to solve TLM: 0=modified Newton, 1=direct
78	ICNTRL(9) = 0 ! TLM Newton Iterations influence
79	ICNTRL(12) = 0 ! TLM Truncation Error influence
80
81	!~~~> if optional parameters are given, and if they are >0,
82	! then use them to overwrite default settings
83	IF (PRESENT(ICNTRL_U)) THEN
84	WHERE(ICNTRL_U(:) > 0) ICNTRL(:) = ICNTRL_U(:)
85	END IF
86	IF (PRESENT(RCNTRL_U)) THEN
87	WHERE(RCNTRL_U(:) > 0) RCNTRL(:) = RCNTRL_U(:)
88	END IF
89
90
91	T1 = TIN; T2 = TOUT
92	CALL SdirkTLM( NVAR, NTLM, T1, T2, Y, Y_tlm, RTOL, ATOL, &
93	RTOL_tlm, ATOL_tlm, RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR )
94
95	!~~~> Debug option: print number of steps
96	! Ntotal = Ntotal + ISTATUS(Nstp)
97	! PRINT*,'NSTEPS=',ISTATUS(Nstp),' (',Ntotal,')',' O3=', VAR(ind_O3)
98
99	! if optional parameters are given for output
100	! use them to store information in them
101	IF (PRESENT(ISTATUS_U)) ISTATUS_U(:) = ISTATUS(:)
102	IF (PRESENT(RSTATUS_U)) RSTATUS_U(:) = RSTATUS(:)
103	IF (PRESENT(IERR_U)) IERR_U = IERR
104
105	IF (IERR < 0) THEN
106	PRINT *,'SDIRK-TLM: Unsuccessful exit at T=', TIN,' (IERR=',IERR,')'
107	ENDIF
108
109	END SUBROUTINE INTEGRATE_TLM
110
111
112
113	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
114	SUBROUTINE SdirkTLM(N, NTLM, Tinitial, Tfinal, Y, Y_tlm, RelTol, AbsTol, &
115	RelTol_tlm, AbsTol_tlm, &
116	RCNTRL, ICNTRL, RSTATUS, ISTATUS, Ierr)
117	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
118	!
119	! Solves the system y'=F(t,y) using a Singly-Diagonally-Implicit
120	! Runge-Kutta (SDIRK) method.
121	!
122	! This implementation is based on the book and the code Sdirk4:
123	!
124	! E. Hairer and G. Wanner
125	! "Solving ODEs II. Stiff and differential-algebraic problems".
126	! Springer series in computational mathematics, Springer-Verlag, 1996.
127	! This code is based on the SDIRK4 routine in the above book.
128	!
129	! Methods:
130	! * Sdirk 2a, 2b: L-stable, 2 stages, order 2
131	! * Sdirk 3a: L-stable, 3 stages, order 2, adjoint-invariant
132	! * Sdirk 4a, 4b: L-stable, 5 stages, order 4
133	!
134	! (C) Adrian Sandu, July 2005
135	! Virginia Polytechnic Institute and State University
136	! Contact: sandu@cs.vt.edu
137	! Revised by Philipp Miehe and Adrian Sandu, May 2006
138	! This implementation is part of KPP - the Kinetic PreProcessor
139	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
140	!
141	!~~~> INPUT ARGUMENTS:
142	!
143	!- Y(NVAR) = vector of initial conditions (at T=Tinitial)
144	!- [Tinitial,Tfinal] = time range of integration
145	! (if Tinitial>Tfinal the integration is performed backwards in time)
146	!- RelTol, AbsTol = user precribed accuracy
147	!- SUBROUTINE ode_Fun( T, Y, Ydot ) = ODE function,
148	! returns Ydot = Y' = F(T,Y)
149	!- SUBROUTINE ode_Fun( T, Y, Ydot ) = Jacobian of the ODE function,
150	! returns Jcb = dF/dY
151	!- ICNTRL(1:20) = integer inputs parameters
152	!- RCNTRL(1:20) = real inputs parameters
153	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
154	!
155	!~~~> OUTPUT ARGUMENTS:
156	!
157	!- Y(NVAR) -> vector of final states (at T->Tfinal)
158	!- ISTATUS(1:20) -> integer output parameters
159	!- RSTATUS(1:20) -> real output parameters
160	!- Ierr -> job status upon return
161	! success (positive value) or
162	! failure (negative value)
163	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
164	!
165	!~~~> INPUT PARAMETERS:
166	!
167	! Note: For input parameters equal to zero the default values of the
168	! corresponding variables are used.
169	!
170	! Note: For input parameters equal to zero the default values of the
171	! corresponding variables are used.
172	!~~~>
173	! ICNTRL(1) = not used
174	!
175	! ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors
176	! = 1: AbsTol, RelTol are scalars
177	!
178	! ICNTRL(3) = Method
179	!
180	! ICNTRL(4) -> maximum number of integration steps
181	! For ICNTRL(4)=0 the default value of 100000 is used
182	!
183	! ICNTRL(5) -> maximum number of Newton iterations
184	! For ICNTRL(5)=0 the default value of 8 is used
185	!
186	! ICNTRL(6) -> starting values of Newton iterations:
187	! ICNTRL(6)=0 : starting values are interpolated (the default)
188	! ICNTRL(6)=1 : starting values are zero
189	!
190	! ICNTRL(7) -> method to solve TLM equations:
191	! ICNTRL(7)=0 : modified Newton re-using LU (the default)
192	! ICNTRL(7)=1 : direct solution (additional one LU factorization per stage)
193	!
194	! ICNTRL(9) -> switch for TLM Newton iteration error estimation strategy
195	! ICNTRL(9) = 0: base number of iterations as forward solution
196	! ICNTRL(9) = 1: use RTOL_tlm and ATOL_tlm to calculate
197	! error estimation for TLM at Newton stages
198	!
199	! ICNTRL(12) -> switch for TLM truncation error control
200	! ICNTRL(12) = 0: TLM error is not used
201	! ICNTRL(12) = 1: TLM error is computed and used
202	!
203	!
204	!~~~> Real parameters
205	!
206	! RCNTRL(1) -> Hmin, lower bound for the integration step size
207	! It is strongly recommended to keep Hmin = ZERO
208	! RCNTRL(2) -> Hmax, upper bound for the integration step size
209	! RCNTRL(3) -> Hstart, starting value for the integration step size
210	!
211	! RCNTRL(4) -> FacMin, lower bound on step decrease factor (default=0.2)
212	! RCNTRL(5) -> FacMax, upper bound on step increase factor (default=6)
213	! RCNTRL(6) -> FacRej, step decrease factor after multiple rejections
214	! (default=0.1)
215	! RCNTRL(7) -> FacSafe, by which the new step is slightly smaller
216	! than the predicted value (default=0.9)
217	! RCNTRL(8) -> ThetaMin. If Newton convergence rate smaller
218	! than ThetaMin the Jacobian is not recomputed;
219	! (default=0.001)
220	! RCNTRL(9) -> NewtonTol, stopping criterion for Newton's method
221	! (default=0.03)
222	! RCNTRL(10) -> Qmin
223	! RCNTRL(11) -> Qmax. If Qmin < Hnew/Hold < Qmax, then the
224	! step size is kept constant and the LU factorization
225	! reused (default Qmin=1, Qmax=1.2)
226	!
227	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
228	!
229	!~~~> OUTPUT PARAMETERS:
230	!
231	! Note: each call to Rosenbrock adds the current no. of fcn calls
232	! to previous value of ISTATUS(1), and similar for the other params.
233	! Set ISTATUS(1:10) = 0 before call to avoid this accumulation.
234	!
235	! ISTATUS(1) = No. of function calls
236	! ISTATUS(2) = No. of jacobian calls
237	! ISTATUS(3) = No. of steps
238	! ISTATUS(4) = No. of accepted steps
239	! ISTATUS(5) = No. of rejected steps (except at the beginning)
240	! ISTATUS(6) = No. of LU decompositions
241	! ISTATUS(7) = No. of forward/backward substitutions
242	! ISTATUS(8) = No. of singular matrix decompositions
243	!
244	! RSTATUS(1) -> Texit, the time corresponding to the
245	! computed Y upon return
246	! RSTATUS(2) -> Hexit,last accepted step before return
247	! RSTATUS(3) -> Hnew, last predicted step before return
248	! For multiple restarts, use Hnew as Hstart in the following run
249	!
250	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
251	IMPLICIT NONE
252
253	! Arguments
254	INTEGER, INTENT(IN) :: N, NTLM, ICNTRL(20)
255	KPP_REAL, INTENT(IN) :: Tinitial, Tfinal, &
256	RelTol(N), AbsTol(N), RCNTRL(20), &
257	RelTol_tlm(N,NTLM), AbsTol_tlm(N,NTLM)
258	KPP_REAL, INTENT(INOUT) :: Y(NVAR), Y_tlm(N,NTLM)
259	INTEGER, INTENT(OUT) :: Ierr
260	INTEGER, INTENT(INOUT) :: ISTATUS(20)
261	KPP_REAL, INTENT(OUT) :: RSTATUS(20)
262
263	!~~~> SDIRK method coefficients, up to 5 stages
264	INTEGER, PARAMETER :: Smax = 5
265	INTEGER, PARAMETER :: S2A=1, S2B=2, S3A=3, S4A=4, S4B=5
266	KPP_REAL :: rkGamma, rkA(Smax,Smax), rkB(Smax), rkC(Smax), &
267	rkD(Smax), rkE(Smax), rkBhat(Smax), rkELO, &
268	rkAlpha(Smax,Smax), rkTheta(Smax,Smax)
269	INTEGER :: sdMethod, rkS ! The number of stages
270
271	! Local variables
272	KPP_REAL :: Hmin, Hmax, Hstart, Roundoff, &
273	FacMin, Facmax, FacSafe, FacRej, &
274	ThetaMin, NewtonTol, Qmin, Qmax
275	INTEGER :: ITOL, NewtonMaxit, Max_no_steps, i
276	LOGICAL :: StartNewton, DirectTLM, TLMNewtonEst, TLMtruncErr
277	KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
278
279
280	Ierr = 0
281	ISTATUS(1:20) = 0
282	RSTATUS(1:20) = ZERO
283
284	!~~~> For Scalar tolerances (ICNTRL(2).NE.0) the code uses AbsTol(1) and RelTol(1)
285	! For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:NVAR) and RelTol(1:NVAR)
286	IF (ICNTRL(2) == 0) THEN
287	ITOL = 1
288	ELSE
289	ITOL = 0
290	END IF
291
292	!~~~> ICNTRL(3) - method selection
293	SELECT CASE (ICNTRL(3))
294	CASE (0,1)
295	CALL Sdirk2a
296	CASE (2)
297	CALL Sdirk2b
298	CASE (3)
299	CALL Sdirk3a
300	CASE (4)
301	CALL Sdirk4a
302	CASE (5)
303	CALL Sdirk4b
304	CASE DEFAULT
305	CALL Sdirk2a
306	END SELECT
307
308	!~~~> The maximum number of time steps admitted
309	IF (ICNTRL(4) == 0) THEN
310	Max_no_steps = 200000
311	ELSEIF (ICNTRL(4) > 0) THEN
312	Max_no_steps=ICNTRL(4)
313	ELSE
314	PRINT * ,'User-selected ICNTRL(4)=',ICNTRL(4)
315	CALL SDIRK_ErrorMsg(-1,Tinitial,ZERO,Ierr)
316	END IF
317	!~~~> The maximum number of Newton iterations admitted
318	IF(ICNTRL(5) == 0)THEN
319	NewtonMaxit=8
320	ELSE
321	NewtonMaxit=ICNTRL(5)
322	IF(NewtonMaxit <= 0)THEN
323	PRINT * ,'User-selected ICNTRL(5)=',ICNTRL(5)
324	CALL SDIRK_ErrorMsg(-2,Tinitial,ZERO,Ierr)
325	END IF
326	END IF
327	!~~~> StartNewton: Use extrapolation for starting values of Newton iterations
328	IF (ICNTRL(6) == 0) THEN
329	StartNewton = .TRUE.
330	ELSE
331	StartNewton = .FALSE.
332	END IF
333	!~~~> Solve TLM equations directly or by Newton iterations
334	DirectTLM = (ICNTRL(7) == 1)
335	!~~~> Newton iteration error control selection
336	IF (ICNTRL(9) == 0) THEN
337	TLMNewtonEst = .FALSE.
338	ELSE
339	TLMNewtonEst = .TRUE.
340	END IF
341	!~~~> TLM truncation error control selection
342	IF (ICNTRL(12) == 0) THEN
343	TLMtruncErr = .FALSE.
344	ELSE
345	TLMtruncErr = .TRUE.
346	END IF
347
348	!~~~> Unit roundoff (1+Roundoff>1)
349	Roundoff = WLAMCH('E')
350
351	!~~~> Lower bound on the step size: (positive value)
352	IF (RCNTRL(1) == ZERO) THEN
353	Hmin = ZERO
354	ELSEIF (RCNTRL(1) > ZERO) THEN
355	Hmin = RCNTRL(1)
356	ELSE
357	PRINT * , 'User-selected RCNTRL(1)=', RCNTRL(1)
358	CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr)
359	END IF
360
361	!~~~> Upper bound on the step size: (positive value)
362	IF (RCNTRL(2) == ZERO) THEN
363	Hmax = ABS(Tfinal-Tinitial)
364	ELSEIF (RCNTRL(2) > ZERO) THEN
365	Hmax = MIN(ABS(RCNTRL(2)),ABS(Tfinal-Tinitial))
366	ELSE
367	PRINT * , 'User-selected RCNTRL(2)=', RCNTRL(2)
368	CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr)
369	END IF
370
371	!~~~> Starting step size: (positive value)
372	IF (RCNTRL(3) == ZERO) THEN
373	Hstart = MAX(Hmin,Roundoff)
374	ELSEIF (RCNTRL(3) > ZERO) THEN
375	Hstart = MIN(ABS(RCNTRL(3)),ABS(Tfinal-Tinitial))
376	ELSE
377	PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3)
378	CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,Ierr)
379	END IF
380
381	!~~~> Step size can be changed s.t. FacMin < Hnew/Hexit < FacMax
382	IF (RCNTRL(4) == ZERO) THEN
383	FacMin = 0.2_dp
384	ELSEIF (RCNTRL(4) > ZERO) THEN
385	FacMin = RCNTRL(4)
386	ELSE
387	PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4)
388	CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr)
389	END IF
390	IF (RCNTRL(5) == ZERO) THEN
391	FacMax = 10.0_dp
392	ELSEIF (RCNTRL(5) > ZERO) THEN
393	FacMax = RCNTRL(5)
394	ELSE
395	PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5)
396	CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr)
397	END IF
398	!~~~> FacRej: Factor to decrease step after 2 succesive rejections
399	IF (RCNTRL(6) == ZERO) THEN
400	FacRej = 0.1_dp
401	ELSEIF (RCNTRL(6) > ZERO) THEN
402	FacRej = RCNTRL(6)
403	ELSE
404	PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6)
405	CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr)
406	END IF
407	!~~~> FacSafe: Safety Factor in the computation of new step size
408	IF (RCNTRL(7) == ZERO) THEN
409	FacSafe = 0.9_dp
410	ELSEIF (RCNTRL(7) > ZERO) THEN
411	FacSafe = RCNTRL(7)
412	ELSE
413	PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7)
414	CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,Ierr)
415	END IF
416
417	!~~~> ThetaMin: decides whether the Jacobian should be recomputed
418	IF(RCNTRL(8) == 0.D0)THEN
419	ThetaMin = 1.0d-3
420	ELSE
421	ThetaMin = RCNTRL(8)
422	END IF
423
424	!~~~> Stopping criterion for Newton's method
425	IF(RCNTRL(9) == ZERO)THEN
426	NewtonTol = 3.0d-2
427	ELSE
428	NewtonTol = RCNTRL(9)
429	END IF
430
431	!~~~> Qmin, Qmax: IF Qmin < Hnew/Hold < Qmax, STEP SIZE = CONST.
432	IF(RCNTRL(10) == ZERO)THEN
433	Qmin=ONE
434	ELSE
435	Qmin=RCNTRL(10)
436	END IF
437	IF(RCNTRL(11) == ZERO)THEN
438	Qmax=1.2D0
439	ELSE
440	Qmax=RCNTRL(11)
441	END IF
442
443	!~~~> Check if tolerances are reasonable
444	IF (ITOL == 0) THEN
445	IF (AbsTol(1) <= ZERO .OR. RelTol(1) <= 10.D0*Roundoff) THEN
446	PRINT * , ' Scalar AbsTol = ',AbsTol(1)
447	PRINT * , ' Scalar RelTol = ',RelTol(1)
448	CALL SDIRK_ErrorMsg(-5,Tinitial,ZERO,Ierr)
449	END IF
450	ELSE
451	DO i=1,N
452	IF (AbsTol(i) <= 0.D0.OR.RelTol(i) <= 10.D0*Roundoff) THEN
453	PRINT * , ' AbsTol(',i,') = ',AbsTol(i)
454	PRINT * , ' RelTol(',i,') = ',RelTol(i)
455	CALL SDIRK_ErrorMsg(-5,Tinitial,ZERO,Ierr)
456	END IF
457	END DO
458	END IF
459
460	IF (Ierr < 0) RETURN
461
462	CALL SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr )
463
464
465	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
466	CONTAINS ! PROCEDURES INTERNAL TO SDIRK
467	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
468
469
470	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
471	SUBROUTINE SDIRK_IntegratorTLM( N,NTLM,Tinitial,Tfinal,Y,Y_tlm,Ierr )
472	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
473
474	USE KPP_ROOT_Parameters
475	IMPLICIT NONE
476
477	!~~~> Arguments:
478	INTEGER, INTENT(IN) :: N, NTLM
479	KPP_REAL, INTENT(INOUT) :: Y(N), Y_tlm(N,NTLM)
480	KPP_REAL, INTENT(IN) :: Tinitial, Tfinal
481	INTEGER, INTENT(OUT) :: Ierr
482
483	!~~~> Local variables:
484	KPP_REAL :: Z(NVAR,rkS), G(NVAR), TMP(NVAR), &
485	NewtonRate, SCAL(NVAR), DZ(NVAR), &
486	T, H, Theta, Hratio, NewtonPredictedErr, &
487	Qnewton, Err, Fac, Hnew, Tdirection, &
488	NewtonIncrement, NewtonIncrementOld, &
489	SCAL_tlm(NVAR), Yerr(N), Yerr_tlm(N,NTLM), ThetaTLM
490	KPP_REAL :: Z_tlm(NVAR,rkS,NTLM)
491	INTEGER :: itlm, j, IER, istage, NewtonIter, saveNiter, NewtonIterTLM
492	INTEGER :: IP(NVAR), IP_tlm(NVAR)
493	LOGICAL :: Reject, FirstStep, SkipJac, SkipLU, NewtonDone
494
495	#ifdef FULL_ALGEBRA
496	KPP_REAL, DIMENSION(NVAR,NVAR) :: FJAC, E, Jac, E_tlm
497	#else
498	KPP_REAL, DIMENSION(LU_NONZERO) :: FJAC, E, Jac, E_tlm
499	#endif
500	KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
501
502
503	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
504	!~~~> Initializations
505	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
506
507	T = Tinitial
508	Tdirection = SIGN(ONE,Tfinal-Tinitial)
509	H = MAX(ABS(Hmin),ABS(Hstart))
510	IF (ABS(H) <= 10.D0*Roundoff) H=1.0D-6
511	H=MIN(ABS(H),Hmax)
512	H=SIGN(H,Tdirection)
513	SkipLU = .FALSE.
514	SkipJac = .FALSE.
515	Reject = .FALSE.
516	FirstStep=.TRUE.
517
518	CALL SDIRK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL)
519
520	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
521	!~~~> Time loop begins
522	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
523	Tloop: DO WHILE ( (Tfinal-T)*Tdirection - Roundoff > ZERO )
524
525
526	!~~~> Compute E = 1/(h*gamma)-Jac and its LU decomposition
527	IF ( .NOT.SkipLU ) THEN ! This time around skip the Jac update and LU
528	CALL SDIRK_PrepareMatrix ( H, T, Y, FJAC, &
529	SkipJac, SkipLU, E, IP, Reject, IER )
530	IF (IER /= 0) THEN
531	CALL SDIRK_ErrorMsg(-8,T,H,Ierr); RETURN
532	END IF
533	END IF
534
535	IF (ISTATUS(Nstp) > Max_no_steps) THEN
536	CALL SDIRK_ErrorMsg(-6,T,H,Ierr); RETURN
537	END IF
538	IF ( (T+0.1d0*H == T) .OR. (ABS(H) <= Roundoff) ) THEN
539	CALL SDIRK_ErrorMsg(-7,T,H,Ierr); RETURN
540	END IF
541
542	stages:DO istage = 1, rkS
543
544	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
545	!~~~> Simplified Newton iterations
546	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
547
548	!~~~> Starting values for Newton iterations
549	CALL Set2zero(N,Z(1,istage))
550
551	!~~~> Prepare the loop-independent part of the right-hand side
552	CALL Set2zero(N,G)
553	IF (istage > 1) THEN
554	DO j = 1, istage-1
555	! Gj(:) = sum_j Theta(i,j)Zj(:) = H sum_j A(i,j)*Fun(Zj(:))
556	CALL WAXPY(N,rkTheta(istage,j),Z(1,j),1,G,1)
557	! Zi(:) = sum_j Alpha(i,j)*Zj(:)
558	IF (StartNewton) THEN
559	CALL WAXPY(N,rkAlpha(istage,j),Z(1,j),1,Z(1,istage),1)
560	END IF
561	END DO
562	END IF
563
564	!~~~> Initializations for Newton iteration
565	NewtonDone = .FALSE.
566	Fac = 0.5d0 ! Step reduction factor if too many iterations
567
568	NewtonLoop:DO NewtonIter = 1, NewtonMaxit
569
570	!~~~> Prepare the loop-dependent part of the right-hand side
571	CALL WADD(N,Y,Z(1,istage),TMP) ! TMP <- Y + Zi
572	CALL FUN_CHEM(T+rkC(istage)*H,TMP,DZ) ! DZ <- Fun(Y+Zi)
573	ISTATUS(Nfun) = ISTATUS(Nfun) + 1
574	! DZ(1:N) = G(1:N) - Z(1:N,istage) + (HrkGamma)DZ(1:N)
575	CALL WSCAL(N, H*rkGamma, DZ, 1)
576	CALL WAXPY (N, -ONE, Z(1,istage), 1, DZ, 1)
577	CALL WAXPY (N, ONE, G,1, DZ,1)
578
579	!~~~> Solve the linear system
580	CALL SDIRK_Solve ( H, N, E, IP, IER, DZ )
581
582	!~~~> Check convergence of Newton iterations
583	CALL SDIRK_ErrorNorm(N, DZ, SCAL, NewtonIncrement)
584	IF ( NewtonIter == 1 ) THEN
585	Theta = ABS(ThetaMin)
586	NewtonRate = 2.0d0
587	ELSE
588	Theta = NewtonIncrement/NewtonIncrementOld
589	IF (Theta < 0.99d0) THEN
590	NewtonRate = Theta/(ONE-Theta)
591	! Predict error at the end of Newton process
592	NewtonPredictedErr = NewtonIncrement &
593	Theta*(NewtonMaxit-NewtonIter)/(ONE-Theta)
594	IF (NewtonPredictedErr >= NewtonTol) THEN
595	! Non-convergence of Newton: predicted error too large
596	Qnewton = MIN(10.0d0,NewtonPredictedErr/NewtonTol)
597	Fac = 0.8d0Qnewton*(-ONE/(1+NewtonMaxit-NewtonIter))
598	EXIT NewtonLoop
599	END IF
600	ELSE ! Non-convergence of Newton: Theta too large
601	EXIT NewtonLoop
602	END IF
603	END IF
604	NewtonIncrementOld = NewtonIncrement
605	! Update solution: Z(:) <-- Z(:)+DZ(:)
606	CALL WAXPY(N,ONE,DZ,1,Z(1,istage),1)
607
608	! Check error in Newton iterations
609	NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol)
610	IF (NewtonDone) THEN
611	! Tune error in TLM variables by defining the minimal number of Newton iterations.
612	saveNiter = NewtonIter+1
613	EXIT NewtonLoop
614	END IF
615
616	END DO NewtonLoop
617
618	IF (.NOT.NewtonDone) THEN
619	!CALL RK_ErrorMsg(-12,T,H,Ierr);
620	H = Fac*H; Reject=.TRUE.
621	SkipJac = .TRUE.; SkipLU = .FALSE.
622	CYCLE Tloop
623	END IF
624
625	!~~~> End of simplified Newton iterations for forward variables
626
627	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
628	!~~~> Solve for TLM variables
629	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
630
631	!~~~> Direct solution for TLM variables
632	DirTLM:IF (DirectTLM) THEN
633
634	TMP(1:N) = Y(1:N) + Z(1:N,istage)
635	SkipJac = .FALSE.
636	CALL SDIRK_PrepareMatrix ( H, T+rkC(istage)*H, TMP, Jac, &
637	SkipJac, SkipLU, E_tlm, IP_tlm, Reject, IER )
638	IF (IER /= 0) CYCLE TLoop
639
640	TlmL: DO itlm = 1, NTLM
641	G(1:N) = Y_tlm(1:N,itlm)
642	IF (istage > 1) THEN
643	! Gj(:) = sum_j Theta(i,j)*Zj_tlm(:)
644	! = H * sum_j A(i,j)Jac(Zj(:))(Yj_tlm+Zj_tlm)
645	DO j = 1, istage-1
646	CALL WAXPY(N,rkTheta(istage,j),Z_tlm(1,j,itlm),1,G,1)
647	END DO
648	END IF
649	CALL SDIRK_Solve ( H, N, E_tlm, IP_tlm, IER, G )
650	Z_tlm(1:N,istage,itlm) = G(1:N) - Y_tlm(1:N,itlm)
651	END DO TlmL
652
653	ELSE DirTLM
654
655	!~~~> Jacobian of the current stage solution
656	TMP(1:N) = Y(1:N) + Z(1:N,istage)
657	CALL JAC_CHEM(T+rkC(istage)*H,TMP,Jac)
658	ISTATUS(Njac) = ISTATUS(Njac) + 1
659
660	!~~~> Simplified Newton iterations for TLM variables
661	TlmLoop:DO itlm = 1,NTLM
662	NewtonRate = MAX(NewtonRate,Roundoff)**0.8d0
663
664	!~~~> Starting values for Newton iterations
665	CALL Set2zero(N,Z_tlm(1,istage,itlm))
666
667	!~~~> Prepare the loop-independent part of the right-hand side
668	#ifdef FULL_ALGEBRA
669	DZ = MATMUL(Jac,Y_tlm(1,itlm)) ! DZ <- Jac(Y+Z)*Y_tlm
670	#else
671	CALL Jac_SP_Vec ( Jac, Y_tlm(1,itlm), DZ )
672	#endif
673	G(1:N) = (HrkGamma)DZ(1:N)
674	IF (istage > 1) THEN
675	! Gj(:) = sum_j Theta(i,j)*Zj_tlm(:)
676	! = H * sum_j A(i,j)Jac(Zj(:))(Yj_tlm+Zj_tlm)
677	DO j = 1, istage-1
678	CALL WAXPY(N,rkTheta(istage,j),Z_tlm(1,j,itlm),1,G,1)
679	END DO
680	END IF
681
682
683	!~~~> Initializations for Newton iteration
684	IF (TLMNewtonEst) THEN
685	NewtonDone = .FALSE.
686	Fac = 0.5d0 ! Step reduction factor if too many iterations
687
688	CALL SDIRK_ErrorScale(N,ITOL,AbsTol_tlm(1,itlm),RelTol_tlm(1,itlm), &
689	Y_tlm(1,itlm),SCAL_tlm)
690	END IF
691
692	NewtonLoopTLM:DO NewtonIterTLM = 1, NewtonMaxit
693
694	!~~~> Prepare the loop-dependent part of the right-hand side
695	#ifdef FULL_ALGEBRA
696	DZ = MATMUL(Jac,Z_tlm(1,istage,itlm)) ! DZ <- Jac(Y+Z)*Z_tlm
697	#else
698	CALL Jac_SP_Vec ( Jac, Z_tlm(1,istage,itlm), DZ )
699	#endif
700	DZ(1:N) = (HrkGamma)DZ(1:N)+G(1:N)-Z_tlm(1:N,istage,itlm)
701
702	CALL SDIRK_Solve ( H, N, E, IP, IER, DZ )
703
704	IF (TLMNewtonEst) THEN
705	!~~~> Check convergence of Newton iterations
706	CALL SDIRK_ErrorNorm(N, DZ, SCAL_tlm, NewtonIncrement)
707	IF ( NewtonIterTLM <= 1 ) THEN
708	ThetaTLM = ABS(ThetaMin)
709	NewtonRate = 2.0d0
710	ELSE
711	ThetaTLM = NewtonIncrement/NewtonIncrementOld
712	IF (ThetaTLM < 0.99d0) THEN
713	NewtonRate = ThetaTLM/(ONE-ThetaTLM)
714	! Predict error at the end of Newton process
715	NewtonPredictedErr = NewtonIncrement &
716	ThetaTLM*(NewtonMaxit-NewtonIterTLM)/(ONE-ThetaTLM)
717	IF (NewtonPredictedErr >= NewtonTol) THEN
718	! Non-convergence of Newton: predicted error too large
719	Qnewton = MIN(10.0d0,NewtonPredictedErr/NewtonTol)
720	Fac = 0.8d0Qnewton*(-ONE/(1+NewtonMaxit-NewtonIterTLM))
721	EXIT NewtonLoopTLM
722	END IF
723	ELSE ! Non-convergence of Newton: Theta too large
724	EXIT NewtonLoopTLM
725	END IF
726	END IF
727	NewtonIncrementOld = NewtonIncrement
728	END IF !(TLMNewtonEst)
729
730	! Update solution: Z_tlm(:) <-- Z_tlm(:)+DZ(:)
731	CALL WAXPY(N,ONE,DZ,1,Z_tlm(1,istage,itlm),1)
732
733	! Check error in Newton iterations
734	IF (TLMNewtonEst) THEN
735	NewtonDone = (NewtonRate*NewtonIncrement <= NewtonTol)
736	IF (NewtonDone) EXIT NewtonLoopTLM
737	ELSE
738	! Minimum number of iterations same as FWD iterations
739	IF (NewtonIterTLM>=saveNiter) EXIT NewtonLoopTLM
740	END IF
741
742	END DO NewtonLoopTLM
743
744	IF ((TLMNewtonEst) .AND. (.NOT.NewtonDone)) THEN
745	!CALL RK_ErrorMsg(-12,T,H,Ierr);
746	H = Fac*H; Reject=.TRUE.
747	SkipJac = .TRUE.; SkipLU = .FALSE.
748	CYCLE Tloop
749	END IF
750
751	END DO TlmLoop
752	!~~~> End of simplified Newton iterations for TLM
753
754	END IF DirTLM
755
756	END DO stages
757
758	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
759	!~~~> Error estimation
760	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
761	ISTATUS(Nstp) = ISTATUS(Nstp) + 1
762	CALL Set2zero(N,Yerr)
763	DO i = 1,rkS
764	IF (rkE(i)/=ZERO) CALL WAXPY(N,rkE(i),Z(1,i),1,Yerr,1)
765	END DO
766
767	CALL SDIRK_Solve ( H, N, E, IP, IER, Yerr )
768	CALL SDIRK_ErrorNorm(N, Yerr, SCAL, Err)
769
770	IF (TLMtruncErr) THEN
771	CALL Set2zero(NVAR*NTLM,Yerr_tlm)
772	DO itlm=1,NTLM
773	DO j=1,rkS
774	IF (rkE(j) /= ZERO) CALL WAXPY(N,rkE(j),Z_tlm(1,j,itlm),1,Yerr_tlm(1,itlm),1)
775	END DO
776	CALL SDIRK_Solve (H, N, E, IP, IER, Yerr_tlm(1,itlm))
777	END DO
778	CALL SDIRK_ErrorNorm_tlm(N,NTLM, Yerr_tlm, Err)
779	END IF
780
781	!~~~> Computation of new step size Hnew
782	Fac = FacSafe(Err)*(-ONE/rkELO)
783	Fac = MAX(FacMin,MIN(FacMax,Fac))
784	Hnew = H*Fac
785
786	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
787	!~~~> Accept/Reject step
788	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
789	accept: IF ( Err < ONE ) THEN !~~~> Step is accepted
790
791	FirstStep=.FALSE.
792	ISTATUS(Nacc) = ISTATUS(Nacc) + 1
793
794	!~~~> Update time and solution
795	T = T + H
796	! Y(:) <-- Y(:) + Sum_j rkD(j)*Z_j(:)
797	DO i = 1,rkS
798	IF (rkD(i)/=ZERO) THEN
799	CALL WAXPY(N,rkD(i),Z(1,i),1,Y,1)
800	DO itlm = 1, NTLM
801	CALL WAXPY(N,rkD(i),Z_tlm(1,i,itlm),1,Y_tlm(1,itlm),1)
802	END DO
803	END IF
804	END DO
805
806	!~~~> Update scaling coefficients
807	CALL SDIRK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL)
808
809	!~~~> Next time step
810	Hnew = Tdirection*MIN(ABS(Hnew),Hmax)
811	! Last T and H
812	RSTATUS(Ntexit) = T
813	RSTATUS(Nhexit) = H
814	RSTATUS(Nhnew) = Hnew
815	! No step increase after a rejection
816	IF (Reject) Hnew = Tdirection*MIN(ABS(Hnew),ABS(H))
817	Reject = .FALSE.
818	IF ((T+Hnew/Qmin-Tfinal)*Tdirection > ZERO) THEN
819	H = Tfinal-T
820	ELSE
821	Hratio=Hnew/H
822	! If step not changed too much keep Jacobian and reuse LU
823	SkipLU = ( (Theta <= ThetaMin) .AND. (Hratio >= Qmin) &
824	.AND. (Hratio <= Qmax) )
825	! For TLM: do not skip LU (decrease TLM error)
826	SkipLU = .FALSE.
827	IF (.NOT.SkipLU) H = Hnew
828	END IF
829	! If convergence is fast enough, do not update Jacobian
830	! SkipJac = (Theta <= ThetaMin)
831	SkipJac = .FALSE.
832
833	ELSE accept !~~~> Step is rejected
834
835	IF (FirstStep .OR. Reject) THEN
836	H = FacRej*H
837	ELSE
838	H = Hnew
839	END IF
840	Reject = .TRUE.
841	SkipJac = .TRUE.
842	SkipLU = .FALSE.
843	IF (ISTATUS(Nacc) >= 1) ISTATUS(Nrej) = ISTATUS(Nrej) + 1
844
845	END IF accept
846
847	END DO Tloop
848
849	! Successful return
850	Ierr = 1
851
852	END SUBROUTINE SDIRK_IntegratorTLM
853
854
855	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
856	SUBROUTINE SDIRK_ErrorScale(N, ITOL, AbsTol, RelTol, Y, SCAL)
857	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
858	IMPLICIT NONE
859	INTEGER :: N, i, ITOL
860	KPP_REAL :: AbsTol(NVAR), RelTol(NVAR), &
861	Y(NVAR), SCAL(NVAR)
862	IF (ITOL == 0) THEN
863	DO i=1,N
864	SCAL(i) = ONE / ( AbsTol(1)+RelTol(1)*ABS(Y(i)) )
865	END DO
866	ELSE
867	DO i=1,N
868	SCAL(i) = ONE / ( AbsTol(i)+RelTol(i)*ABS(Y(i)) )
869	END DO
870	END IF
871	END SUBROUTINE SDIRK_ErrorScale
872
873
874	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
875	SUBROUTINE SDIRK_ErrorNorm(N, Y, SCAL, Err)
876	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
877	!
878	INTEGER :: i, N
879	KPP_REAL :: Y(N), SCAL(N), Err
880	Err = ZERO
881	DO i=1,N
882	Err = Err+(Y(i)SCAL(i))*2
883	END DO
884	Err = MAX( SQRT(Err/DBLE(N)), 1.0d-10 )
885	!
886	END SUBROUTINE SDIRK_ErrorNorm
887
888	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
889	SUBROUTINE SDIRK_ErrorNorm_tlm(N,NTLM, Y_tlm, FWD_Err)
890	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
891	!
892	INTEGER :: itlm, NTLM, N
893	KPP_REAL :: Y_tlm(N,NTLM), SCAL_tlm(N), FWD_Err, Err
894
895	DO itlm=1,NTLM
896	CALL SDIRK_ErrorScale(N,ITOL,AbsTol_tlm(1,itlm),RelTol_tlm(1,itlm), &
897	Y_tlm(1,itlm),SCAL_tlm)
898	CALL SDIRK_ErrorNorm(N, Y_tlm(1,itlm), SCAL_tlm, Err)
899	FWD_Err = MAX(FWD_Err, Err)
900	END DO
901	!
902	END SUBROUTINE SDIRK_ErrorNorm_tlm
903
904
905	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
906	SUBROUTINE SDIRK_ErrorMsg(Code,T,H,Ierr)
907	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
908	! Handles all error messages
909	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
910
911	KPP_REAL, INTENT(IN) :: T, H
912	INTEGER, INTENT(IN) :: Code
913	INTEGER, INTENT(OUT) :: Ierr
914
915	Ierr = Code
916	PRINT * , &
917	'Forced exit from SDIRK due to the following error:'
918
919	SELECT CASE (Code)
920	CASE (-1)
921	PRINT * , '--> Improper value for maximal no of steps'
922	CASE (-2)
923	PRINT * , '--> Improper value for maximal no of Newton iterations'
924	CASE (-3)
925	PRINT * , '--> Hmin/Hmax/Hstart must be positive'
926	CASE (-4)
927	PRINT * , '--> FacMin/FacMax/FacRej must be positive'
928	CASE (-5)
929	PRINT * , '--> Improper tolerance values'
930	CASE (-6)
931	PRINT * , '--> No of steps exceeds maximum bound'
932	CASE (-7)
933	PRINT * , '--> Step size too small: T + 10*H = T', &
934	' or H < Roundoff'
935	CASE (-8)
936	PRINT * , '--> Matrix is repeatedly singular'
937	CASE DEFAULT
938	PRINT *, 'Unknown Error code: ', Code
939	END SELECT
940
941	PRINT *, "T=", T, "and H=", H
942
943	END SUBROUTINE SDIRK_ErrorMsg
944
945
946	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
947	SUBROUTINE SDIRK_PrepareMatrix ( H, T, Y, FJAC, &
948	SkipJac, SkipLU, E, IP, Reject, ISING )
949	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
950	!~~~> Compute the matrix E = I - 1/(HGamma)Jac, and its decomposition
951	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
952
953	IMPLICIT NONE
954
955	KPP_REAL, INTENT(INOUT) :: H
956	KPP_REAL, INTENT(IN) :: T, Y(NVAR)
957	LOGICAL, INTENT(INOUT) :: SkipJac,SkipLU,Reject
958	INTEGER, INTENT(OUT) :: ISING, IP(NVAR)
959	#ifdef FULL_ALGEBRA
960	KPP_REAL, INTENT(INOUT) :: FJAC(NVAR,NVAR)
961	KPP_REAL, INTENT(OUT) :: E(NVAR,NVAR)
962	#else
963	KPP_REAL, INTENT(INOUT) :: FJAC(LU_NONZERO)
964	KPP_REAL, INTENT(OUT) :: E(LU_NONZERO)
965	#endif
966	KPP_REAL :: HGammaInv
967	INTEGER :: i, j, ConsecutiveSng
968
969	ConsecutiveSng = 0
970	ISING = 1
971
972	Hloop: DO WHILE (ISING /= 0)
973
974	HGammaInv = ONE/(H*rkGamma)
975
976	!~~~> Compute the Jacobian
977	! IF (SkipJac) THEN
978	! SkipJac = .FALSE.
979	! ELSE
980	IF (.NOT. SkipJac) THEN
981	CALL JAC_CHEM( T, Y, FJAC )
982	ISTATUS(Njac) = ISTATUS(Njac) + 1
983	END IF
984
985	#ifdef FULL_ALGEBRA
986	DO j=1,NVAR
987	DO i=1,NVAR
988	E(i,j) = -FJAC(i,j)
989	END DO
990	E(j,j) = E(j,j)+HGammaInv
991	END DO
992	CALL DGETRF( NVAR, NVAR, E, NVAR, IP, ISING )
993	#else
994	DO i = 1,LU_NONZERO
995	E(i) = -FJAC(i)
996	END DO
997	DO i = 1,NVAR
998	j = LU_DIAG(i); E(j) = E(j) + HGammaInv
999	END DO
1000	CALL KppDecomp ( E, ISING)
1001	IP(1) = 1
1002	#endif
1003	ISTATUS(Ndec) = ISTATUS(Ndec) + 1
1004
1005	IF (ISING /= 0) THEN
1006	WRITE (6,*) ' MATRIX IS SINGULAR, ISING=',ISING,' T=',T,' H=',H
1007	ISTATUS(Nsng) = ISTATUS(Nsng) + 1; ConsecutiveSng = ConsecutiveSng + 1
1008	IF (ConsecutiveSng >= 6) RETURN ! Failure
1009	H = 0.5d0*H
1010	SkipJac = .TRUE.
1011	SkipLU = .FALSE.
1012	Reject = .TRUE.
1013	END IF
1014
1015	END DO Hloop
1016
1017	END SUBROUTINE SDIRK_PrepareMatrix
1018
1019	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1020	SUBROUTINE SDIRK_Solve ( H, N, E, IP, ISING, RHS )
1021	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1022	!~~~> Solves the system (HGamma-Jac)x = R
1023	! using the LU decomposition of E = I - 1/(HGamma)Jac
1024	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1025	IMPLICIT NONE
1026	INTEGER, INTENT(IN) :: N, IP(N), ISING
1027	KPP_REAL, INTENT(IN) :: H
1028	#ifdef FULL_ALGEBRA
1029	KPP_REAL, INTENT(IN) :: E(NVAR,NVAR)
1030	#else
1031	KPP_REAL, INTENT(IN) :: E(LU_NONZERO)
1032	#endif
1033	KPP_REAL, INTENT(INOUT) :: RHS(N)
1034	KPP_REAL :: HGammaInv
1035
1036	HGammaInv = ONE/(H*rkGamma)
1037	CALL WSCAL(N,HGammaInv,RHS,1)
1038	#ifdef FULL_ALGEBRA
1039	CALL DGETRS( 'N', N, 1, E, N, IP, RHS, N, ISING )
1040	#else
1041	CALL KppSolve(E, RHS)
1042	#endif
1043	ISTATUS(Nsol) = ISTATUS(Nsol) + 1
1044
1045	END SUBROUTINE SDIRK_Solve
1046
1047	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1048	SUBROUTINE Sdirk4a
1049	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1050	sdMethod = S4A
1051	! Number of stages
1052	rkS = 5
1053
1054	! Method coefficients
1055	rkGamma = .2666666666666666666666666666666667d0
1056
1057	rkA(1,1) = .2666666666666666666666666666666667d0
1058	rkA(2,1) = .5000000000000000000000000000000000d0
1059	rkA(2,2) = .2666666666666666666666666666666667d0
1060	rkA(3,1) = .3541539528432732316227461858529820d0
1061	rkA(3,2) = -.5415395284327323162274618585298197d-1
1062	rkA(3,3) = .2666666666666666666666666666666667d0
1063	rkA(4,1) = .8515494131138652076337791881433756d-1
1064	rkA(4,2) = -.6484332287891555171683963466229754d-1
1065	rkA(4,3) = .7915325296404206392428857585141242d-1
1066	rkA(4,4) = .2666666666666666666666666666666667d0
1067	rkA(5,1) = 2.100115700566932777970612055999074d0
1068	rkA(5,2) = -.7677800284445976813343102185062276d0
1069	rkA(5,3) = 2.399816361080026398094746205273880d0
1070	rkA(5,4) = -2.998818699869028161397714709433394d0
1071	rkA(5,5) = .2666666666666666666666666666666667d0
1072
1073	rkB(1) = 2.100115700566932777970612055999074d0
1074	rkB(2) = -.7677800284445976813343102185062276d0
1075	rkB(3) = 2.399816361080026398094746205273880d0
1076	rkB(4) = -2.998818699869028161397714709433394d0
1077	rkB(5) = .2666666666666666666666666666666667d0
1078
1079	rkBhat(1)= 2.885264204387193942183851612883390d0
1080	rkBhat(2)= -.1458793482962771337341223443218041d0
1081	rkBhat(3)= 2.390008682465139866479830743628554d0
1082	rkBhat(4)= -4.129393538556056674929560012190140d0
1083	rkBhat(5)= 0.d0
1084
1085	rkC(1) = .2666666666666666666666666666666667d0
1086	rkC(2) = .7666666666666666666666666666666667d0
1087	rkC(3) = .5666666666666666666666666666666667d0
1088	rkC(4) = .3661315380631796996374935266701191d0
1089	rkC(5) = 1.d0
1090
1091	! Ynew = Yold + hSum_i {rkB_ik_i} = Yold + Sum_i {rkD_i*Z_i}
1092	rkD(1) = 0.d0
1093	rkD(2) = 0.d0
1094	rkD(3) = 0.d0
1095	rkD(4) = 0.d0
1096	rkD(5) = 1.d0
1097
1098	! Err = h * Sum_i {(rkB_i-rkBhat_i)k_i} = Sum_i {rkE_iZ_i}
1099	rkE(1) = -.6804000050475287124787034884002302d0
1100	rkE(2) = 1.558961944525217193393931795738823d0
1101	rkE(3) = -13.55893003128907927748632408763868d0
1102	rkE(4) = 15.48522576958521253098585004571302d0
1103	rkE(5) = 1.d0
1104
1105	! Local order of Err estimate
1106	rkElo = 4
1107
1108	! hSum_j {rkA_ijk_j} = Sum_j {rkTheta_ij*Z_j}
1109	rkTheta(2,1) = 1.875000000000000000000000000000000d0
1110	rkTheta(3,1) = 1.708847304091539528432732316227462d0
1111	rkTheta(3,2) = -.2030773231622746185852981969486824d0
1112	rkTheta(4,1) = .2680325578937783958847157206823118d0
1113	rkTheta(4,2) = -.1828840955527181631794050728644549d0
1114	rkTheta(4,3) = .2968246986151577397160821594427966d0
1115	rkTheta(5,1) = .9096171815241460655379433581446771d0
1116	rkTheta(5,2) = -3.108254967778352416114774430509465d0
1117	rkTheta(5,3) = 12.33727431701306195581826123274001d0
1118	rkTheta(5,4) = -11.24557012450885560524143016037523d0
1119
1120	! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j}
1121	rkAlpha(2,1) = 2.875000000000000000000000000000000d0
1122	rkAlpha(3,1) = .8500000000000000000000000000000000d0
1123	rkAlpha(3,2) = .4434782608695652173913043478260870d0
1124	rkAlpha(4,1) = .7352046091658870564637910527807370d0
1125	rkAlpha(4,2) = -.9525565003057343527941920657462074d-1
1126	rkAlpha(4,3) = .4290111305453813852259481840631738d0
1127	rkAlpha(5,1) = -16.10898993405067684831655675112808d0
1128	rkAlpha(5,2) = 6.559571569643355712998131800797873d0
1129	rkAlpha(5,3) = -15.90772144271326504260996815012482d0
1130	rkAlpha(5,4) = 25.34908987169226073668861694892683d0
1131
1132	!~~~> Coefficients for continuous solution
1133	! rkD(1,1)= 24.74416644927758d0
1134	! rkD(1,2)= -4.325375951824688d0
1135	! rkD(1,3)= 41.39683763286316d0
1136	! rkD(1,4)= -61.04144619901784d0
1137	! rkD(1,5)= -3.391332232917013d0
1138	! rkD(2,1)= -51.98245719616925d0
1139	! rkD(2,2)= 10.52501981094525d0
1140	! rkD(2,3)= -154.2067922191855d0
1141	! rkD(2,4)= 214.3082125319825d0
1142	! rkD(2,5)= 14.71166018088679d0
1143	! rkD(3,1)= 33.14347947522142d0
1144	! rkD(3,2)= -19.72986789558523d0
1145	! rkD(3,3)= 230.4878502285804d0
1146	! rkD(3,4)= -287.6629744338197d0
1147	! rkD(3,5)= -18.99932366302254d0
1148	! rkD(4,1)= -5.905188728329743d0
1149	! rkD(4,2)= 13.53022403646467d0
1150	! rkD(4,3)= -117.6778956422581d0
1151	! rkD(4,4)= 134.3962081008550d0
1152	! rkD(4,5)= 8.678995715052762d0
1153	!
1154	! DO i=1,4 ! CONTi <-- Sum_j rkD(i,j)*Zj
1155	! CALL Set2zero(N,CONT(1,i))
1156	! DO j = 1,rkS
1157	! CALL WAXPY(N,rkD(i,j),Z(1,j),1,CONT(1,i),1)
1158	! END DO
1159	! END DO
1160
1161	rkELO = 4.0d0
1162
1163	END SUBROUTINE Sdirk4a
1164
1165	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1166	SUBROUTINE Sdirk4b
1167	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1168	sdMethod = S4B
1169	! Number of stages
1170	rkS = 5
1171
1172	! Method coefficients
1173	rkGamma = .25d0
1174
1175	rkA(1,1) = 0.25d0
1176	rkA(2,1) = 0.5d00
1177	rkA(2,2) = 0.25d0
1178	rkA(3,1) = 0.34d0
1179	rkA(3,2) =-0.40d-1
1180	rkA(3,3) = 0.25d0
1181	rkA(4,1) = 0.2727941176470588235294117647058824d0
1182	rkA(4,2) =-0.5036764705882352941176470588235294d-1
1183	rkA(4,3) = 0.2757352941176470588235294117647059d-1
1184	rkA(4,4) = 0.25d0
1185	rkA(5,1) = 1.041666666666666666666666666666667d0
1186	rkA(5,2) =-1.020833333333333333333333333333333d0
1187	rkA(5,3) = 7.812500000000000000000000000000000d0
1188	rkA(5,4) =-7.083333333333333333333333333333333d0
1189	rkA(5,5) = 0.25d0
1190
1191	rkB(1) = 1.041666666666666666666666666666667d0
1192	rkB(2) = -1.020833333333333333333333333333333d0
1193	rkB(3) = 7.812500000000000000000000000000000d0
1194	rkB(4) = -7.083333333333333333333333333333333d0
1195	rkB(5) = 0.250000000000000000000000000000000d0
1196
1197	rkBhat(1)= 1.069791666666666666666666666666667d0
1198	rkBhat(2)= -0.894270833333333333333333333333333d0
1199	rkBhat(3)= 7.695312500000000000000000000000000d0
1200	rkBhat(4)= -7.083333333333333333333333333333333d0
1201	rkBhat(5)= 0.212500000000000000000000000000000d0
1202
1203	rkC(1) = 0.25d0
1204	rkC(2) = 0.75d0
1205	rkC(3) = 0.55d0
1206	rkC(4) = 0.50d0
1207	rkC(5) = 1.00d0
1208
1209	! Ynew = Yold + hSum_i {rkB_ik_i} = Yold + Sum_i {rkD_i*Z_i}
1210	rkD(1) = 0.0d0
1211	rkD(2) = 0.0d0
1212	rkD(3) = 0.0d0
1213	rkD(4) = 0.0d0
1214	rkD(5) = 1.0d0
1215
1216	! Err = h * Sum_i {(rkB_i-rkBhat_i)k_i} = Sum_i {rkE_iZ_i}
1217	rkE(1) = 0.5750d0
1218	rkE(2) = 0.2125d0
1219	rkE(3) = -4.6875d0
1220	rkE(4) = 4.2500d0
1221	rkE(5) = 0.1500d0
1222
1223	! Local order of Err estimate
1224	rkElo = 4
1225
1226	! hSum_j {rkA_ijk_j} = Sum_j {rkTheta_ij*Z_j}
1227	rkTheta(2,1) = 2.d0
1228	rkTheta(3,1) = 1.680000000000000000000000000000000d0
1229	rkTheta(3,2) = -.1600000000000000000000000000000000d0
1230	rkTheta(4,1) = 1.308823529411764705882352941176471d0
1231	rkTheta(4,2) = -.1838235294117647058823529411764706d0
1232	rkTheta(4,3) = 0.1102941176470588235294117647058824d0
1233	rkTheta(5,1) = -3.083333333333333333333333333333333d0
1234	rkTheta(5,2) = -4.291666666666666666666666666666667d0
1235	rkTheta(5,3) = 34.37500000000000000000000000000000d0
1236	rkTheta(5,4) = -28.33333333333333333333333333333333d0
1237
1238	! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j}
1239	rkAlpha(2,1) = 3.
1240	rkAlpha(3,1) = .8800000000000000000000000000000000d0
1241	rkAlpha(3,2) = .4400000000000000000000000000000000d0
1242	rkAlpha(4,1) = .1666666666666666666666666666666667d0
1243	rkAlpha(4,2) = -.8333333333333333333333333333333333d-1
1244	rkAlpha(4,3) = .9469696969696969696969696969696970d0
1245	rkAlpha(5,1) = -6.d0
1246	rkAlpha(5,2) = 9.d0
1247	rkAlpha(5,3) = -56.81818181818181818181818181818182d0
1248	rkAlpha(5,4) = 54.d0
1249
1250	END SUBROUTINE Sdirk4b
1251
1252	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1253	SUBROUTINE Sdirk2a
1254	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1255	sdMethod = S2A
1256	! Number of stages
1257	rkS = 2
1258
1259	! Method coefficients
1260	rkGamma = .2928932188134524755991556378951510d0
1261
1262	rkA(1,1) = .2928932188134524755991556378951510d0
1263	rkA(2,1) = .7071067811865475244008443621048490d0
1264	rkA(2,2) = .2928932188134524755991556378951510d0
1265
1266	rkB(1) = .7071067811865475244008443621048490d0
1267	rkB(2) = .2928932188134524755991556378951510d0
1268
1269	rkBhat(1)= .6666666666666666666666666666666667d0
1270	rkBhat(2)= .3333333333333333333333333333333333d0
1271
1272	rkC(1) = 0.292893218813452475599155637895151d0
1273	rkC(2) = 1.0d0
1274
1275	! Ynew = Yold + hSum_i {rkB_ik_i} = Yold + Sum_i {rkD_i*Z_i}
1276	rkD(1) = 0.0d0
1277	rkD(2) = 1.0d0
1278
1279	! Err = h * Sum_i {(rkB_i-rkBhat_i)k_i} = Sum_i {rkE_iZ_i}
1280	rkE(1) = 0.4714045207910316829338962414032326d0
1281	rkE(2) = -0.1380711874576983496005629080698993d0
1282
1283	! Local order of Err estimate
1284	rkElo = 2
1285
1286	! hSum_j {rkA_ijk_j} = Sum_j {rkTheta_ij*Z_j}
1287	rkTheta(2,1) = 2.414213562373095048801688724209698d0
1288
1289	! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j}
1290	rkAlpha(2,1) = 3.414213562373095048801688724209698d0
1291
1292	END SUBROUTINE Sdirk2a
1293
1294	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1295	SUBROUTINE Sdirk2b
1296	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1297	sdMethod = S2B
1298	! Number of stages
1299	rkS = 2
1300
1301	! Method coefficients
1302	rkGamma = 1.707106781186547524400844362104849d0
1303
1304	rkA(1,1) = 1.707106781186547524400844362104849d0
1305	rkA(2,1) = -.707106781186547524400844362104849d0
1306	rkA(2,2) = 1.707106781186547524400844362104849d0
1307
1308	rkB(1) = -.707106781186547524400844362104849d0
1309	rkB(2) = 1.707106781186547524400844362104849d0
1310
1311	rkBhat(1)= .6666666666666666666666666666666667d0
1312	rkBhat(2)= .3333333333333333333333333333333333d0
1313
1314	rkC(1) = 1.707106781186547524400844362104849d0
1315	rkC(2) = 1.0d0
1316
1317	! Ynew = Yold + hSum_i {rkB_ik_i} = Yold + Sum_i {rkD_i*Z_i}
1318	rkD(1) = 0.0d0
1319	rkD(2) = 1.0d0
1320
1321	! Err = h * Sum_i {(rkB_i-rkBhat_i)k_i} = Sum_i {rkE_iZ_i}
1322	rkE(1) = -.4714045207910316829338962414032326d0
1323	rkE(2) = .8047378541243650162672295747365659d0
1324
1325	! Local order of Err estimate
1326	rkElo = 2
1327
1328	! hSum_j {rkA_ijk_j} = Sum_j {rkTheta_ij*Z_j}
1329	rkTheta(2,1) = -.414213562373095048801688724209698d0
1330
1331	! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j}
1332	rkAlpha(2,1) = .5857864376269049511983112757903019d0
1333
1334	END SUBROUTINE Sdirk2b
1335
1336
1337	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1338	SUBROUTINE Sdirk3a
1339	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1340	sdMethod = S3A
1341	! Number of stages
1342	rkS = 3
1343
1344	! Method coefficients
1345	rkGamma = .2113248654051871177454256097490213d0
1346
1347	rkA(1,1) = .2113248654051871177454256097490213d0
1348	rkA(2,1) = .2113248654051871177454256097490213d0
1349	rkA(2,2) = .2113248654051871177454256097490213d0
1350	rkA(3,1) = .2113248654051871177454256097490213d0
1351	rkA(3,2) = .5773502691896257645091487805019573d0
1352	rkA(3,3) = .2113248654051871177454256097490213d0
1353
1354	rkB(1) = .2113248654051871177454256097490213d0
1355	rkB(2) = .5773502691896257645091487805019573d0
1356	rkB(3) = .2113248654051871177454256097490213d0
1357
1358	rkBhat(1)= .2113248654051871177454256097490213d0
1359	rkBhat(2)= .6477918909913548037576239837516312d0
1360	rkBhat(3)= .1408832436034580784969504064993475d0
1361
1362	rkC(1) = .2113248654051871177454256097490213d0
1363	rkC(2) = .4226497308103742354908512194980427d0
1364	rkC(3) = 1.d0
1365
1366	! Ynew = Yold + hSum_i {rkB_ik_i} = Yold + Sum_i {rkD_i*Z_i}
1367	rkD(1) = 0.d0
1368	rkD(2) = 0.d0
1369	rkD(3) = 1.d0
1370
1371	! Err = h * Sum_i {(rkB_i-rkBhat_i)k_i} = Sum_i {rkE_iZ_i}
1372	rkE(1) = 0.9106836025229590978424821138352906d0
1373	rkE(2) = -1.244016935856292431175815447168624d0
1374	rkE(3) = 0.3333333333333333333333333333333333d0
1375
1376	! Local order of Err estimate
1377	rkElo = 2
1378
1379	! hSum_j {rkA_ijk_j} = Sum_j {rkTheta_ij*Z_j}
1380	rkTheta(2,1) = 1.0d0
1381	rkTheta(3,1) = -1.732050807568877293527446341505872d0
1382	rkTheta(3,2) = 2.732050807568877293527446341505872d0
1383
1384	! Starting value for Newton iterations: Z_i^0 = Sum_j {rkAlpha_ij*Z_j}
1385	rkAlpha(2,1) = 2.0d0
1386	rkAlpha(3,1) = -12.92820323027550917410978536602349d0
1387	rkAlpha(3,2) = 8.83012701892219323381861585376468d0
1388
1389	END SUBROUTINE Sdirk3a
1390
1391	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1392	END SUBROUTINE SdirkTLM ! AND ALL ITS INTERNAL PROCEDURES
1393	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1394
1395
1396	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1397	SUBROUTINE FUN_CHEM( T, Y, P )
1398	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1399
1400	USE KPP_ROOT_Parameters, ONLY: NVAR
1401	USE KPP_ROOT_Global, ONLY: TIME, FIX, RCONST
1402	USE KPP_ROOT_Function, ONLY: Fun
1403	USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO
1404
1405	KPP_REAL :: T, Told
1406	KPP_REAL :: Y(NVAR), P(NVAR)
1407
1408	Told = TIME
1409	TIME = T
1410	CALL Update_SUN()
1411	CALL Update_RCONST()
1412
1413	CALL Fun( Y, FIX, RCONST, P )
1414
1415	TIME = Told
1416
1417	END SUBROUTINE FUN_CHEM
1418
1419
1420	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1421	SUBROUTINE JAC_CHEM( T, Y, JV )
1422	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1423
1424
1425	USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO
1426	USE KPP_ROOT_Global, ONLY: TIME, FIX, RCONST
1427	USE KPP_ROOT_Jacobian, ONLY: Jac_SP,LU_IROW,LU_ICOL
1428	USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO
1429
1430	KPP_REAL :: T, Told
1431	KPP_REAL :: Y(NVAR)
1432	#ifdef FULL_ALGEBRA
1433	KPP_REAL :: JS(LU_NONZERO), JV(NVAR,NVAR)
1434	INTEGER :: i, j
1435	#else
1436	KPP_REAL :: JV(LU_NONZERO)
1437	#endif
1438
1439	Told = TIME
1440	TIME = T
1441	CALL Update_SUN()
1442	CALL Update_RCONST()
1443
1444	#ifdef FULL_ALGEBRA
1445	CALL Jac_SP(Y, FIX, RCONST, JS)
1446	DO j=1,NVAR
1447	DO i=1,NVAR
1448	JV(i,j) = 0.0D0
1449	END DO
1450	END DO
1451	DO i=1,LU_NONZERO
1452	JV(LU_IROW(i),LU_ICOL(i)) = JS(i)
1453	END DO
1454	#else
1455	CALL Jac_SP(Y, FIX, RCONST, JV)
1456	#endif
1457	TIME = Told
1458
1459	END SUBROUTINE JAC_CHEM
1460
1461	END MODULE KPP_ROOT_Integrator

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

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/sdirk_tlm.f90 @ 2968

Download in other formats: