Home

Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sdirk.f90 @ 4592

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

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

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

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int.modified_WCOPY/sdirk.f90 @ 4592

Download in other formats: