Home

Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sdirk_soa.f90 @ 4901

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

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

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

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int.modified_WCOPY/sdirk_soa.f90 @ 4901

Download in other formats: