Home

Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sdirk_adj.f90 @ 4829

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

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

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

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/sdirk_adj.f90 @ 4829

Download in other formats: