Home

Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

kpp_sdirk4.f90 @ 4086

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

Line
1	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
2	! SDIRK - Singly-Diagonally-Implicit Runge-Kutta method !
3	! (L-stable, 5 stages, order 4) !
4	! By default the code employs the KPP sparse linear algebra routines !
5	! Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK) !
6	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
7	! A. Sandu - version of July 10, 2005
8
9	MODULE KPP_ROOT_Integrator
10
11	USE KPP_ROOT_Precision
12	USE KPP_ROOT_Global, ONLY: FIX, RCONST, TIME
13	USE KPP_ROOT_Parameters, ONLY: NVAR, NSPEC, NFIX, LU_NONZERO
14	USE KPP_ROOT_JacobianSP, ONLY: LU_DIAG
15	USE KPP_ROOT_LinearAlgebra, ONLY: KppDecomp, KppSolve, &
16	Set2zero, WLAMCH, WAXPY, WCOPY
17
18	IMPLICIT NONE
19	PUBLIC
20	SAVE
21
22	!~~~> Statistics on the work performed by the SDIRK method
23	INTEGER :: Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng
24	INTEGER, PARAMETER :: ifun=1, ijac=2, istp=3, iacc=4, &
25	irej=5, idec=6, isol=7, isng=8, itexit=1, ihexit=2
26	! SDIRK method coefficients
27	KPP_REAL :: rkAlpha(5,4), rkBeta(5,4), rkD(4,5), &
28	rkGamma, rkA(5,5), rkB(5), rkC(5)
29
30	! description of the error numbers IERR
31	CHARACTER(LEN=50), PARAMETER, DIMENSION(-8:1) :: IERR_NAMES = (/ &
32	'Matrix is repeatedly singular ', & ! -8
33	'Step size too small: T + 10*H = T or H < Roundoff ', & ! -7
34	'No of steps exceeds maximum bound ', & ! -6
35	'Improper tolerance values ', & ! -5
36	'FacMin/FacMax/FacRej must be positive ', & ! -4
37	'Hmin/Hmax/Hstart must be positive ', & ! -3
38	'Improper value for maximal no of Newton iterations', & ! -2
39	'Improper value for maximal no of steps ', & ! -1
40	' ', & ! 0 (not used)
41	'Success ' /) ! 1
42
43	CONTAINS
44
45	SUBROUTINE INTEGRATE( TIN, TOUT, &
46	ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U, IERR_U )
47
48	USE KPP_ROOT_Parameters
49	USE KPP_ROOT_Global
50	IMPLICIT NONE
51
52	KPP_REAL, INTENT(IN) :: TIN ! Start Time
53	KPP_REAL, INTENT(IN) :: TOUT ! End Time
54	! Optional input parameters and statistics
55	INTEGER, INTENT(IN), OPTIONAL :: ICNTRL_U(20)
56	KPP_REAL, INTENT(IN), OPTIONAL :: RCNTRL_U(20)
57	INTEGER, INTENT(OUT), OPTIONAL :: ISTATUS_U(20)
58	KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20)
59	INTEGER, INTENT(OUT), OPTIONAL :: IERR_U
60
61	!INTEGER, SAVE :: Ntotal = 0
62	KPP_REAL :: RCNTRL(20), RSTATUS(20)
63	INTEGER :: ICNTRL(20), ISTATUS(20), IERR
64
65	ICNTRL(:) = 0
66	RCNTRL(:) = 0.0_dp
67	ISTATUS(:) = 0
68
69	! If optional parameters are given, and if they are >0,
70	! then they overwrite default settings.
71	IF (PRESENT(ICNTRL_U)) THEN
72	WHERE(ICNTRL_U(:) > 0) ICNTRL(:) = ICNTRL_U(:)
73	END IF
74	IF (PRESENT(RCNTRL_U)) THEN
75	WHERE(RCNTRL_U(:) > 0) RCNTRL(:) = RCNTRL_U(:)
76	END IF
77
78	CALL SDIRK( NVAR,TIN,TOUT,VAR,RTOL,ATOL, &
79	RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR )
80
81	! mz_rs_20050716: IERR and ISTATUS(istp) are returned to the user who then
82	! decides what to do about it, i.e. either stop the run or ignore it.
83	!!$ IF (IERR < 0) THEN
84	!!$ PRINT *,'SDIRK: Unsuccessful exit at T=',TIN,' (IERR=',IERR,')'
85	!!$ ENDIF
86	!!$ Ntotal = Ntotal + Nstp
87	!!$ PRINT*,'NSTEPS=',Nstp, '(',Ntotal,')'
88
89	STEPMIN = RSTATUS(ihexit) ! Save last step
90
91	! if optional parameters are given for output they to return information
92	IF (PRESENT(ISTATUS_U)) ISTATUS_U(:) = ISTATUS(:)
93	IF (PRESENT(RSTATUS_U)) THEN
94	RSTATUS_U(:) = RSTATUS(:)
95	RSTATUS_U(1) = TOUT ! final time
96	END IF
97	IF (PRESENT(IERR_U)) IERR_U = IERR
98
99	END SUBROUTINE INTEGRATE
100
101	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
102	SUBROUTINE SDIRK(N, Tinitial, Tfinal, Y, RelTol, AbsTol, &
103	RCNTRL, ICNTRL, RSTATUS, ISTATUS, IDID)
104	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
105	!
106	! Solves the system y'=F(t,y) using a Singly-Diagonally-Implicit
107	! Runge-Kutta (SDIRK) method.
108	!
109	! For details on SDIRK methods and their implementation consult:
110	! E. Hairer and G. Wanner
111	! "Solving ODEs II. Stiff and differential-algebraic problems".
112	! Springer series in computational mathematics, Springer-Verlag, 1996.
113	! This code is based on the SDIRK4 routine in the above book.
114	!
115	! (C) Adrian Sandu, July 2005
116	! Virginia Polytechnic Institute and State University
117	! Contact: sandu@cs.vt.edu
118	! This implementation is part of KPP - the Kinetic PreProcessor
119	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
120	!
121	!~~~> INPUT ARGUMENTS:
122	!
123	!- Y(NVAR) = vector of initial conditions (at T=Tinitial)
124	!- [Tinitial,Tfinal] = time range of integration
125	! (if Tinitial>Tfinal the integration is performed backwards in time)
126	!- RelTol, AbsTol = user precribed accuracy
127	!- SUBROUTINE ode_Fun( T, Y, Ydot ) = ODE function,
128	! returns Ydot = Y' = F(T,Y)
129	!- SUBROUTINE ode_Fun( T, Y, Ydot ) = Jacobian of the ODE function,
130	! returns Jcb = dF/dY
131	!- ICNTRL(1:20) = integer inputs parameters
132	!- RCNTRL(1:20) = real inputs parameters
133	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
134	!
135	!~~~> OUTPUT ARGUMENTS:
136	!
137	!- Y(NVAR) -> vector of final states (at T->Tfinal)
138	!- ISTATUS(1:20) -> integer output parameters
139	!- RSTATUS(1:20) -> real output parameters
140	!- IDID -> job status upon return
141	! success (positive value) or
142	! failure (negative value)
143	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
144	!
145	!~~~> INPUT PARAMETERS:
146	!
147	! Note: For input parameters equal to zero the default values of the
148	! corresponding variables are used.
149	!
150	! Note: For input parameters equal to zero the default values of the
151	! corresponding variables are used.
152	!~~~>
153	! ICNTRL(1) = not used
154	!
155	! ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors
156	! = 1: AbsTol, RelTol are scalars
157	!
158	! ICNTRL(3) = not used
159	!
160	! ICNTRL(4) -> maximum number of integration steps
161	! For ICNTRL(4)=0 the default value of 100000 is used
162	!
163	! ICNTRL(5) -> maximum number of Newton iterations
164	! For ICNTRL(5)=0 the default value of 8 is used
165	!
166	! ICNTRL(6) -> starting values of Newton iterations:
167	! ICNTRL(6)=0 : starting values are interpolated (the default)
168	! ICNTRL(6)=1 : starting values are zero
169	!
170	!~~~> Real parameters
171	!
172	! RCNTRL(1) -> Hmin, lower bound for the integration step size
173	! It is strongly recommended to keep Hmin = ZERO
174	! RCNTRL(2) -> Hmax, upper bound for the integration step size
175	! RCNTRL(3) -> Hstart, starting value for the integration step size
176	!
177	! RCNTRL(4) -> FacMin, lower bound on step decrease factor (default=0.2)
178	! RCNTRL(5) -> FacMax, upper bound on step increase factor (default=6)
179	! RCNTRL(6) -> FacRej, step decrease factor after multiple rejections
180	! (default=0.1)
181	! RCNTRL(7) -> FacSafe, by which the new step is slightly smaller
182	! than the predicted value (default=0.9)
183	! RCNTRL(8) -> ThetaMin. If Newton convergence rate smaller
184	! than ThetaMin the Jacobian is not recomputed;
185	! (default=0.001)
186	! RCNTRL(9) -> NewtonTol, stopping criterion for Newton's method
187	! (default=0.03)
188	! RCNTRL(10) -> Qmin
189	! RCNTRL(11) -> Qmax. If Qmin < Hnew/Hold < Qmax, then the
190	! step size is kept constant and the LU factorization
191	! reused (default Qmin=1, Qmax=1.2)
192	!
193	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
194	!
195	!~~~> OUTPUT PARAMETERS:
196	!
197	! Note: each call to Rosenbrock adds the current no. of fcn calls
198	! to previous value of ISTATUS(1), and similar for the other params.
199	! Set ISTATUS(1:10) = 0 before call to avoid this accumulation.
200	!
201	! ISTATUS(1) = No. of function calls
202	! ISTATUS(2) = No. of jacobian calls
203	! ISTATUS(3) = No. of steps
204	! ISTATUS(4) = No. of accepted steps
205	! ISTATUS(5) = No. of rejected steps (except at the beginning)
206	! ISTATUS(6) = No. of LU decompositions
207	! ISTATUS(7) = No. of forward/backward substitutions
208	! ISTATUS(8) = No. of singular matrix decompositions
209	!
210	! RSTATUS(1) -> Texit, the time corresponding to the
211	! computed Y upon return
212	! RSTATUS(2) -> Hexit, last predicted step before exit
213	! For multiple restarts, use Hexit as Hstart in the following run
214	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
215	IMPLICIT NONE
216
217	! Arguments
218	INTEGER, INTENT(IN) :: N, ICNTRL(20)
219	KPP_REAL, INTENT(IN) :: Tinitial, Tfinal, &
220	RelTol(NVAR), AbsTol(NVAR), RCNTRL(20)
221	KPP_REAL, INTENT(INOUT) :: Y(NVAR)
222	INTEGER, INTENT(OUT) :: IDID
223	INTEGER, INTENT(INOUT) :: ISTATUS(20)
224	KPP_REAL, INTENT(OUT) :: RSTATUS(20)
225
226	! Local variables
227	KPP_REAL :: Hmin, Hmax, Hstart, Roundoff, &
228	FacMin, Facmax, FacSafe, FacRej, &
229	ThetaMin, NewtonTol, Qmin, Qmax, &
230	Texit, Hexit
231	INTEGER :: ITOL, NewtonMaxit, Max_no_steps, i
232	KPP_REAL, PARAMETER :: ZERO = 0.0d0
233
234	!~~~> Initialize statistics
235	Nfun = ISTATUS(ifun)
236	Njac = ISTATUS(ijac)
237	Nstp = ISTATUS(istp)
238	Nacc = ISTATUS(iacc)
239	Nrej = ISTATUS(irej)
240	Ndec = ISTATUS(idec)
241	Nsol = ISTATUS(isol)
242	Nsng = ISTATUS(isng)
243	! Nfun=0; Njac=0; Nstp=0; Nacc=0
244	! Nrej=0; Ndec=0; Nsol=0; Nsng=0
245
246	IDID = 0
247
248	!~~~> For Scalar tolerances (ICNTRL(2).NE.0) the code uses AbsTol(1) and RelTol(1)
249	! For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:NVAR) and RelTol(1:NVAR)
250	IF (ICNTRL(2) == 0) THEN
251	ITOL = 1
252	ELSE
253	ITOL = 0
254	END IF
255
256	!~~~> The maximum number of time steps admitted
257	IF (ICNTRL(3) == 0) THEN
258	Max_no_steps = 100000
259	ELSEIF (Max_no_steps > 0) THEN
260	Max_no_steps=ICNTRL(3)
261	ELSE
262	PRINT * ,'User-selected ICNTRL(3)=',ICNTRL(3)
263	CALL SDIRK_ErrorMsg(-1,Tinitial,ZERO,IDID)
264	END IF
265
266
267	!~~~> The maximum number of Newton iterations admitted
268	IF(ICNTRL(4) == 0)THEN
269	NewtonMaxit=8
270	ELSE
271	NewtonMaxit=ICNTRL(4)
272	IF(NewtonMaxit <= 0)THEN
273	PRINT * ,'User-selected ICNTRL(4)=',ICNTRL(4)
274	CALL SDIRK_ErrorMsg(-2,Tinitial,ZERO,IDID)
275	END IF
276	END IF
277
278	!~~~> Unit roundoff (1+Roundoff>1)
279	Roundoff = WLAMCH('E')
280
281	!~~~> Lower bound on the step size: (positive value)
282	IF (RCNTRL(1) == ZERO) THEN
283	Hmin = ZERO
284	ELSEIF (RCNTRL(1) > ZERO) THEN
285	Hmin = RCNTRL(1)
286	ELSE
287	PRINT * , 'User-selected RCNTRL(1)=', RCNTRL(1)
288	CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,IDID)
289	END IF
290
291	!~~~> Upper bound on the step size: (positive value)
292	IF (RCNTRL(2) == ZERO) THEN
293	Hmax = ABS(Tfinal-Tinitial)
294	ELSEIF (RCNTRL(2) > ZERO) THEN
295	Hmax = MIN(ABS(RCNTRL(2)),ABS(Tfinal-Tinitial))
296	ELSE
297	PRINT * , 'User-selected RCNTRL(2)=', RCNTRL(2)
298	CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,IDID)
299	END IF
300
301	!~~~> Starting step size: (positive value)
302	IF (RCNTRL(3) == ZERO) THEN
303	Hstart = MAX(Hmin,Roundoff)
304	ELSEIF (RCNTRL(3) > ZERO) THEN
305	Hstart = MIN(ABS(RCNTRL(3)),ABS(Tfinal-Tinitial))
306	ELSE
307	PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3)
308	CALL SDIRK_ErrorMsg(-3,Tinitial,ZERO,IDID)
309	END IF
310
311	!~~~> Step size can be changed s.t. FacMin < Hnew/Hexit < FacMax
312	IF (RCNTRL(4) == ZERO) THEN
313	FacMin = 0.2_dp
314	ELSEIF (RCNTRL(4) > ZERO) THEN
315	FacMin = RCNTRL(4)
316	ELSE
317	PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4)
318	CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,IDID)
319	END IF
320	IF (RCNTRL(5) == ZERO) THEN
321	FacMax = 10.0_dp
322	ELSEIF (RCNTRL(5) > ZERO) THEN
323	FacMax = RCNTRL(5)
324	ELSE
325	PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5)
326	CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,IDID)
327	END IF
328	!~~~> FacRej: Factor to decrease step after 2 succesive rejections
329	IF (RCNTRL(6) == ZERO) THEN
330	FacRej = 0.1_dp
331	ELSEIF (RCNTRL(6) > ZERO) THEN
332	FacRej = RCNTRL(6)
333	ELSE
334	PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6)
335	CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,IDID)
336	END IF
337	!~~~> FacSafe: Safety Factor in the computation of new step size
338	IF (RCNTRL(7) == ZERO) THEN
339	FacSafe = 0.9_dp
340	ELSEIF (RCNTRL(7) > ZERO) THEN
341	FacSafe = RCNTRL(7)
342	ELSE
343	PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7)
344	CALL SDIRK_ErrorMsg(-4,Tinitial,ZERO,IDID)
345	END IF
346
347	!~~~> DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED;
348	IF(RCNTRL(8) == 0.D0)THEN
349	ThetaMin = 1.0d-3
350	ELSE
351	ThetaMin = RCNTRL(8)
352	END IF
353
354	!~~~> STOPPING CRITERION FOR NEWTON'S METHOD
355	IF(RCNTRL(9) == 0.0d0)THEN
356	NewtonTol = 3.0d-2
357	ELSE
358	NewtonTol =RCNTRL(9)
359	END IF
360
361	!~~~> Qmin AND Qmax: IF Qmin < Hnew/Hold < Qmax, STEP SIZE = CONST.
362	IF(RCNTRL(10) == 0.D0)THEN
363	Qmin=1.D0
364	ELSE
365	Qmin=RCNTRL(10)
366	END IF
367	IF(RCNTRL(11) == 0.D0)THEN
368	Qmax=1.2D0
369	ELSE
370	Qmax=RCNTRL(11)
371	END IF
372
373	!~~~> Check if tolerances are reasonable
374	IF (ITOL == 0) THEN
375	IF (AbsTol(1) <= 0.D0.OR.RelTol(1) <= 10.D0*Roundoff) THEN
376	PRINT * , ' Scalar AbsTol = ',AbsTol(1)
377	PRINT * , ' Scalar RelTol = ',RelTol(1)
378	CALL SDIRK_ErrorMsg(-5,Tinitial,ZERO,IDID)
379	END IF
380	ELSE
381	DO i=1,N
382	IF (AbsTol(i) <= 0.D0.OR.RelTol(i) <= 10.D0*Roundoff) THEN
383	PRINT * , ' AbsTol(',i,') = ',AbsTol(i)
384	PRINT * , ' RelTol(',i,') = ',RelTol(i)
385	CALL SDIRK_ErrorMsg(-5,Tinitial,ZERO,IDID)
386	END IF
387	END DO
388	END IF
389
390	IF (IDID < 0) RETURN
391
392
393	!~~~> CALL TO CORE INTEGRATOR
394	CALL SDIRK_Integrator( N,Tinitial,Tfinal,Y,Hmin,Hmax,Hstart, &
395	RelTol,AbsTol,ITOL, Max_no_steps, NewtonMaxit, &
396	Roundoff, FacMin, FacMax, FacRej, FacSafe, ThetaMin, &
397	NewtonTol, Qmin, Qmax, Hexit, Texit, IDID )
398
399	!~~~> Collect run statistics
400	ISTATUS(ifun) = Nfun
401	ISTATUS(ijac) = Njac
402	ISTATUS(istp) = Nstp
403	ISTATUS(iacc) = Nacc
404	ISTATUS(irej) = Nrej
405	ISTATUS(idec) = Ndec
406	ISTATUS(isol) = Nsol
407	ISTATUS(isng) = Nsng
408	!~~~> Last T and H
409	RSTATUS(:) = 0.0_dp
410	RSTATUS(itexit) = Texit
411	RSTATUS(ihexit) = Hexit
412
413	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
414	CONTAINS ! PROCEDURES INTERNAL TO SDIRK
415	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
416
417
418	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
419	SUBROUTINE SDIRK_Integrator( N,Tinitial,Tfinal,Y,Hmin,Hmax,Hstart, &
420	RelTol,AbsTol,ITOL, Max_no_steps, NewtonMaxit, &
421	Roundoff, FacMin, FacMax, FacRej, FacSafe, ThetaMin, &
422	NewtonTol, Qmin, Qmax, Hexit, Texit, IDID )
423	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
424	! CORE INTEGRATOR FOR SDIRK4
425	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
426
427	USE KPP_ROOT_Parameters
428	IMPLICIT NONE
429
430	!~~~> Arguments:
431	INTEGER :: N
432	KPP_REAL, INTENT(INOUT) :: Y(NVAR)
433	KPP_REAL, INTENT(IN) :: Tinitial, Tfinal, Hmin, Hmax, Hstart, &
434	RelTol(NVAR), AbsTol(NVAR), Roundoff, &
435	FacMin, FacMax, FacRej, FacSafe, ThetaMin, &
436	NewtonTol, Qmin, Qmax
437	KPP_REAL, INTENT(OUT) :: Hexit, Texit
438	INTEGER, INTENT(IN) :: ITOL, Max_no_steps, NewtonMaxit
439	INTEGER, INTENT(OUT) :: IDID
440
441	!~~~> Local variables:
442	KPP_REAL :: Z(NVAR,5), FV(NVAR,5), CONT(NVAR,4), &
443	NewtonFactor(5), SCAL(NVAR), RHS(NVAR), &
444	G(NVAR), Yhat(NVAR), TMP(NVAR), &
445	T, H, Hold, Theta, Hratio, Hmax1, W, &
446	HGammaInv, DYTH, QNEWT, ERR, Fac, Hnew, &
447	Tdirection, NewtonErr, NewtonErrOld
448	INTEGER :: i, j, ISING, istage, NewtonIter, IP(NVAR)
449	LOGICAL :: Reject, FIRST, NewtonReject, FreshJac, SkipJacUpdate, SkipLU
450
451	#ifdef FULL_ALGEBRA
452	KPP_REAL FJAC(NVAR,NVAR), E(NVAR,NVAR)
453	#else
454	KPP_REAL FJAC(LU_NONZERO), E(LU_NONZERO)
455	#endif
456	KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
457
458	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
459	! INITIALISATIONS
460	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
461	CALL SDIRK_Coefficients
462	T = Tinitial
463	Tdirection = SIGN(1.D0,Tfinal-Tinitial)
464	Hmax1=MIN(ABS(Hmax),ABS(Tfinal-Tinitial))
465	H = MAX(ABS(Hmin),ABS(Hstart))
466	IF (ABS(H) <= 10.D0*Roundoff) H=1.0D-6
467	H=MIN(ABS(H),Hmax1)
468	H=SIGN(H,Tdirection)
469	Hold=H
470	NewtonReject=.FALSE.
471	SkipLU =.FALSE.
472	FreshJac = .FALSE.
473	SkipJacUpdate = .FALSE.
474	Reject=.FALSE.
475	FIRST=.TRUE.
476	NewtonFactor(1:5)=ONE
477
478	CALL SDIRK_ErrorScale(ITOL, AbsTol, RelTol, Y, SCAL)
479
480	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
481	!~~~> Time loop begins
482	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
483	Tloop: DO WHILE ( (Tfinal-T)*Tdirection - Roundoff > ZERO )
484
485
486	!~~~> Compute E = 1/(h*gamma)-Jac and its LU decomposition
487	IF ( SkipLU ) THEN ! This time around skip the Jac update and LU
488	SkipLU = .FALSE.; FreshJac = .FALSE.; SkipJacUpdate = .FALSE.
489	ELSE
490	CALL SDIRK_PrepareMatrix ( H, T, Y, FJAC, &
491	FreshJac, SkipJacUpdate, E, IP, Reject, ISING )
492	IF (ISING /= 0) THEN
493	CALL SDIRK_ErrorMsg(-8,T,H,IDID); RETURN
494	END IF
495	END IF
496
497	IF (Nstp>Max_no_steps) THEN
498	CALL SDIRK_ErrorMsg(-6,T,H,IDID); RETURN
499	END IF
500	IF ( (T+0.1d0*H == T) .OR. (ABS(H) <= Roundoff) ) THEN
501	CALL SDIRK_ErrorMsg(-7,T,H,IDID); RETURN
502	END IF
503
504	HGammaInv = ONE/(H*rkGamma)
505
506	!~~~> NEWTON ITERATION
507	stages:DO istage=1,5
508
509	NewtonFactor(istage) = MAX(NewtonFactor(istage),Roundoff)**0.8d0
510
511	!~~~> STARTING VALUES FOR NEWTON ITERATION
512	CALL Set2zero(N,G)
513	CALL Set2zero(N,Z(1,istage))
514	IF (istage==1) THEN
515	IF (FIRST.OR.NewtonReject) THEN
516	CALL Set2zero(N,Z(1,istage))
517	ELSE
518	W=ONE+rkGamma*H/Hold
519	DO i=1,N
520	Z(i,istage)=W(CONT(i,1)+W(CONT(i,2)+W(CONT(i,3)+WCONT(i,4))))-Yhat(i)
521	END DO
522	END IF
523	ELSE
524	DO j = 1, istage-1
525	! Gj(:) = sum_j Beta(i,j)Zj(:) = H sum_j A(i,j)*Fun(Zj(:))
526	CALL WAXPY(N,rkBeta(istage,j),Z(1,j),1,G,1)
527	! CALL WAXPY(N,H*rkA(istage,j),FV(1,j),1,G,1)
528	! Zi(:) = sum_j Alpha(i,j)*Zj(:)
529	CALL WAXPY(N,rkAlpha(istage,j),Z(1,j),1,Z(1,istage),1)
530	END DO
531	IF (istage==5) CALL WCOPY(N,Z(1,istage),1,Yhat,1) ! Yhat(:) <- Z5(:)
532	END IF
533
534	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
535	! LOOP FOR THE SIMPLIFIED NEWTON ITERATION
536	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
537	NewtonIter=0
538	Theta=ABS(ThetaMin)
539	IF (Reject) Theta=2*ABS(ThetaMin)
540	NewtonErr = 1.0e+6 ! To force-enter Newton loop
541
542	Newton: DO WHILE (NewtonFactor(istage)*NewtonErr > NewtonTol)
543
544	IF (NewtonIter >= NewtonMaxit) THEN
545	H=H*0.5d0
546	Reject=.TRUE.
547	NewtonReject=.TRUE.
548	CYCLE Tloop
549	END IF
550	NewtonIter=NewtonIter+1
551
552	!~~~> COMPUTE THE RIGHT-HAND SIDE
553	TMP(1:N) = Y(1:N) + Z(1:N,istage)
554	CALL FUN_CHEM(T+rkC(istage)*H,TMP,RHS)
555	TMP(1:N) = G(1:N) - Z(1:N,istage)
556	CALL WAXPY(N,HGammaInv,TMP,1,RHS,1) ! RHS(:) <- RHS(:) + HGammaInv*(G(:)-Z(:))
557
558	!~~~> SOLVE THE LINEAR SYSTEMS
559	#ifdef FULL_ALGEBRA
560	CALL DGETRS( 'N', N, 1, E, N, IP, RHS, N, ISING )
561	#else
562	CALL KppSolve(E, RHS)
563	#endif
564	Nsol=Nsol+1
565
566	!~~~> CHECK CONVERGENCE OR IF NUMBER OF ITERATIONS TOO LARGE
567	CALL SDIRK_ErrorNorm(N, RHS, SCAL, NewtonErr)
568	IF ( (NewtonIter >= 2) .AND. (NewtonIter < NewtonMaxit) ) THEN
569	Theta = NewtonErr/NewtonErrOld
570	IF (Theta < 0.99d0) THEN
571	NewtonFactor(istage)=Theta/(ONE-Theta)
572	DYTH = NewtonFactor(istage)NewtonErr &
573	Theta**(NewtonMaxit-1-NewtonIter)
574	QNEWT = MAX(1.0d-4,MIN(16.0d0,DYTH/NewtonTol))
575	IF (QNEWT >= ONE) THEN
576	H=.8D0HQNEWT**(-ONE/(NewtonMaxit-NewtonIter))
577	Reject=.TRUE.
578	NewtonReject=.TRUE.
579	CYCLE Tloop ! go back to the beginning of DO step
580	END IF
581	ELSE
582	NewtonReject=.TRUE.
583	H=H*0.5d0
584	Reject=.TRUE.
585	CYCLE Tloop ! go back to the beginning of DO step
586	END IF
587	END IF
588	NewtonErrOld = NewtonErr
589	CALL WAXPY(N,ONE,RHS,1,Z(1,istage),1) ! Z(:) <-- Z(:)+RHS(:)
590
591	END DO Newton
592
593	!~~> END OF SIMPLIFIED NEWTON ITERATION
594	! Save function values
595	TMP(1:N) = Y(1:N) + Z(1:N,istage)
596	CALL FUN_CHEM(T+rkC(istage)*H,TMP,FV(1,istage))
597
598	END DO stages
599
600
601	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
602	! ERROR ESTIMATION
603	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
604	Nstp=Nstp+1
605	TMP(1:N)=HGammaInv*(Z(1:N,5)-Yhat(1:N))
606
607	#ifdef FULL_ALGEBRA
608	CALL DGETRS( 'N', N, 1, E, N, IP, TMP, N, ISING )
609	#else
610	CALL KppSolve(E, TMP)
611	#endif
612
613	CALL SDIRK_ErrorNorm(N, TMP, SCAL, ERR)
614
615	!~~~> COMPUTATION OF Hnew: WE REQUIRE FacMin <= Hnew/H <= FacMax
616	!Safe = FacSafeDBLE(1+2NewtonMaxit)/DBLE(NewtonIter+2*NewtonMaxit)
617	Fac = MAX(FacMin,MIN(FacMax,(ERR)*(-0.25d0)FacSafe))
618	Hnew = H*Fac
619
620	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
621	! ACCEPT/Reject STEP
622	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~
623	accept: IF ( ERR < ONE ) THEN !~~~> STEP IS ACCEPTED
624
625	FIRST=.FALSE.
626	Nacc=Nacc+1
627	Hold=H
628
629	!~~~> COEFFICIENTS FOR CONTINUOUS SOLUTION
630	CALL WAXPY(N,ONE,Z(1,5),1,Y,1) ! Y(:) <-- Y(:)+Z5(:)
631	CALL WCOPY(N,Z(1,5),1,Yhat,1) ! Yhat <-- Z5
632
633	DO i=1,4 ! CONTi <-- Sum_j rkD(i,j)*Zj
634	CALL Set2zero(N,CONT(1,i))
635	DO j = 1,5
636	CALL WAXPY(N,rkD(i,j),Z(1,j),1,CONT(1,i),1)
637	END DO
638	END DO
639
640	CALL SDIRK_ErrorScale(ITOL, AbsTol, RelTol, Y, SCAL)
641
642	T=T+H
643	FreshJac=.FALSE.
644
645	Hnew = Tdirection*MIN(ABS(Hnew),Hmax1)
646	Hexit = Hnew
647	IF (Reject) Hnew=Tdirection*MIN(ABS(Hnew),ABS(H))
648	Reject = .FALSE.
649	NewtonReject = .FALSE.
650	IF ((T+Hnew/Qmin-Tfinal)*Tdirection > 0.D0) THEN
651	H = Tfinal-T
652	ELSE
653	Hratio=Hnew/H
654	! If step not changed too much, keep it as is;
655	! do not update Jacobian and reuse LU
656	IF ( (Theta <= ThetaMin) .AND. (Hratio >= Qmin) &
657	.AND. (Hratio <= Qmax) ) THEN
658	SkipJacUpdate = .TRUE.
659	SkipLU = .TRUE.
660	ELSE
661	H = Hnew
662	END IF
663	END IF
664	! If convergence is fast enough, do not update Jacobian
665	IF (Theta <= ThetaMin) SkipJacUpdate = .TRUE.
666
667	ELSE accept !~~~> STEP IS REJECTED
668
669	Reject=.TRUE.
670	IF (FIRST) THEN
671	H=H*FacRej
672	ELSE
673	H=Hnew
674	END IF
675	IF (Nacc >= 1) Nrej=Nrej+1
676
677	END IF accept
678
679	END DO Tloop
680
681	! Successful return
682	Texit = T
683	IDID = 1
684
685	END SUBROUTINE SDIRK_Integrator
686
687
688	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
689	SUBROUTINE SDIRK_ErrorScale(ITOL, AbsTol, RelTol, Y, SCAL)
690	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
691	IMPLICIT NONE
692	INTEGER :: i, ITOL
693	KPP_REAL :: AbsTol(NVAR), RelTol(NVAR), &
694	Y(NVAR), SCAL(NVAR)
695	IF (ITOL == 0) THEN
696	DO i=1,NVAR
697	SCAL(i) = 1.0d0 / ( AbsTol(1)+RelTol(1)*ABS(Y(i)) )
698	END DO
699	ELSE
700	DO i=1,NVAR
701	SCAL(i) = 1.0d0 / ( AbsTol(i)+RelTol(i)*ABS(Y(i)) )
702	END DO
703	END IF
704	END SUBROUTINE SDIRK_ErrorScale
705
706
707	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
708	SUBROUTINE SDIRK_ErrorNorm(N, Y, SCAL, ERR)
709	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
710	!
711	INTEGER :: i, N
712	KPP_REAL :: Y(N), SCAL(N), ERR
713	ERR=0.0d0
714	DO i=1,N
715	ERR = ERR+(Y(i)SCAL(i))*2
716	END DO
717	ERR = MAX( SQRT(ERR/DBLE(N)), 1.0d-10 )
718	!
719	END SUBROUTINE SDIRK_ErrorNorm
720
721
722	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
723	SUBROUTINE SDIRK_ErrorMsg(Code,T,H,IERR)
724	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
725	! Handles all error messages
726	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
727
728	KPP_REAL, INTENT(IN) :: T, H
729	INTEGER, INTENT(IN) :: Code
730	INTEGER, INTENT(OUT) :: IERR
731
732	IERR = Code
733	PRINT * , &
734	'Forced exit from SDIRK due to the following error:'
735	IF ((Code>=-8).AND.(Code<=-1)) THEN
736	PRINT *, IERR_NAMES(Code)
737	ELSE
738	PRINT *, 'Unknown Error code: ', Code
739	ENDIF
740
741	PRINT *, "T=", T, "and H=", H
742
743	END SUBROUTINE SDIRK_ErrorMsg
744
745
746	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
747	SUBROUTINE SDIRK_PrepareMatrix ( H, T, Y, FJAC, &
748	FreshJac, SkipJacUpdate, E, IP, Reject, ISING )
749	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
750	!
751	IMPLICIT NONE
752
753	KPP_REAL, INTENT(INOUT) :: H
754	KPP_REAL, INTENT(IN) :: T, Y(NVAR)
755	LOGICAL, INTENT(INOUT) :: FreshJac, SkipJacUpdate
756	INTEGER, INTENT(OUT) :: ISING, IP(NVAR)
757	LOGICAL, INTENT(INOUT) :: Reject
758	#ifdef FULL_ALGEBRA
759	KPP_REAL, INTENT(INOUT) :: FJAC(NVAR,NVAR)
760	KPP_REAL, INTENT(OUT) :: E(NVAR,NVAR)
761	#else
762	KPP_REAL, INTENT(INOUT) :: FJAC(LU_NONZERO)
763	KPP_REAL, INTENT(OUT) :: E(LU_NONZERO)
764	#endif
765	KPP_REAL :: HGammaInv
766	INTEGER :: i, j, ConsecutiveSng
767	KPP_REAL, PARAMETER :: ONE = 1.0d0
768
769	20 CONTINUE
770
771	!~~~> COMPUTE THE JACOBIAN
772	IF (SkipJacUpdate) THEN
773	SkipJacUpdate = .FALSE.
774	ELSE IF ( .NOT.FreshJac ) THEN
775	CALL JAC_CHEM( T, Y, FJAC )
776	FreshJac = .TRUE.
777	END IF
778
779	!~~~> Compute the matrix E = 1/(HGAMMA)Jac, and its decomposition
780	ConsecutiveSng = 0
781	ISING = 1
782
783	Hloop: DO WHILE (ISING /= 0)
784
785	HGammaInv = ONE/(H*rkGamma)
786
787	#ifdef FULL_ALGEBRA
788	DO j=1,NVAR
789	DO i=1,NVAR
790	E(i,j)=-FJAC(i,j)
791	END DO
792	E(j,j)=E(j,j)+HGammaInv
793	END DO
794	CALL DGETRF( NVAR, NVAR, E, NVAR, IP, ISING )
795	#else
796	DO i = 1,LU_NONZERO
797	E(i) = -FJAC(i)
798	END DO
799	DO i = 1,NVAR
800	j = LU_DIAG(i); E(j) = E(j) + HGammaInv
801	END DO
802	CALL KppDecomp ( E, ISING)
803	IP(1) = 1
804	#endif
805	Ndec=Ndec+1
806
807	IF (ISING /= 0) THEN
808	WRITE (6,*) ' MATRIX IS SINGULAR, ISING=',ISING,' T=',T,' H=',H
809	Nsng = Nsng+1; ConsecutiveSng = ConsecutiveSng + 1
810	IF (ConsecutiveSng >= 6) RETURN ! Failure
811	H=H*0.5d0
812	Reject=.TRUE.
813	!~~~> Update Jacobian if not fresh
814	IF ( .NOT.FreshJac ) GOTO 20
815	END IF
816
817	END DO Hloop
818
819	END SUBROUTINE SDIRK_PrepareMatrix
820
821
822	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
823	SUBROUTINE SDIRK_Coefficients
824	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
825	rkGamma=4.0d0/15.0d0
826
827	rkA(1,1)= 4.d0/15.d0
828	rkA(2,1)= 1.d0/2.d0
829	rkA(2,2)= 4.d0/15.d0
830	rkA(3,1)= 51069.d0/144200.d0
831	rkA(3,2)=-7809.d0/144200.d0
832	rkA(3,3)= 4.d0/15.d0
833	rkA(4,1)= 12047244770625658.d0/141474406359725325.d0
834	rkA(4,2)=-3057890203562191.d0/47158135453241775.d0
835	rkA(4,3)= 2239631894905804.d0/28294881271945065.d0
836	rkA(4,4)= 4.d0/15.d0
837	rkA(5,1)= 181513.d0/86430.d0
838	rkA(5,2)=-89074.d0/116015.d0
839	rkA(5,3)= 83636.d0/34851.d0
840	rkA(5,4)=-69863904375173.d0/23297141763930.d0
841	rkA(5,5)= 4.d0/15.d0
842
843	rkB(1)= 181513.d0/86430.d0
844	rkB(2)=-89074.d0/116015.d0
845	rkB(3)= 83636.d0/34851.d0
846	rkB(4)=-69863904375173.d0/23297141763930.d0
847	rkB(5)= 4/15.d0
848
849	rkC(1)=4.d0/15.d0
850	rkC(2)=23.d0/30.d0
851	rkC(3)=17.d0/30.d0
852	rkC(4)=707.d0/1931.d0
853	rkC(5)=1.d0
854
855	rkBeta(2,1)=15.0d0/8.0d0
856	rkBeta(3,1)=1577061.0d0/922880.0d0
857	rkBeta(3,2)=-23427.0d0/115360.0d0
858	rkBeta(4,1)=647163682356923881.0d0/2414496535205978880.0d0
859	rkBeta(4,2)=-593512117011179.0d0/3245291041943520.0d0
860	rkBeta(4,3)=559907973726451.0d0/1886325418129671.0d0
861	rkBeta(5,1)=724545451.0d0/796538880.0d0
862	rkBeta(5,2)=-830832077.0d0/267298560.0d0
863	rkBeta(5,3)=30957577.0d0/2509272.0d0
864	rkBeta(5,4)=-69863904375173.0d0/6212571137048.0d0
865
866	rkAlpha(2,1)= 23.d0/8.d0
867	rkAlpha(3,1)= 0.9838473040915402d0
868	rkAlpha(3,2)= 0.3969226768377252d0
869	rkAlpha(4,1)= 0.6563374010466914d0
870	rkAlpha(4,2)= 0.0d0
871	rkAlpha(4,3)= 0.3372498196189311d0
872	rkAlpha(5,1)=7752107607.0d0/11393456128.0d0
873	rkAlpha(5,2)=-17881415427.0d0/11470078208.0d0
874	rkAlpha(5,3)=2433277665.0d0/179459416.0d0
875	rkAlpha(5,4)=-96203066666797.0d0/6212571137048.0d0
876
877	rkD(1,1)= 24.74416644927758d0
878	rkD(1,2)= -4.325375951824688d0
879	rkD(1,3)= 41.39683763286316d0
880	rkD(1,4)= -61.04144619901784d0
881	rkD(1,5)= -3.391332232917013d0
882	rkD(2,1)= -51.98245719616925d0
883	rkD(2,2)= 10.52501981094525d0
884	rkD(2,3)= -154.2067922191855d0
885	rkD(2,4)= 214.3082125319825d0
886	rkD(2,5)= 14.71166018088679d0
887	rkD(3,1)= 33.14347947522142d0
888	rkD(3,2)= -19.72986789558523d0
889	rkD(3,3)= 230.4878502285804d0
890	rkD(3,4)= -287.6629744338197d0
891	rkD(3,5)= -18.99932366302254d0
892	rkD(4,1)= -5.905188728329743d0
893	rkD(4,2)= 13.53022403646467d0
894	rkD(4,3)= -117.6778956422581d0
895	rkD(4,4)= 134.3962081008550d0
896	rkD(4,5)= 8.678995715052762d0
897
898	END SUBROUTINE SDIRK_Coefficients
899
900	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
901	END SUBROUTINE SDIRK ! AND ALL ITS INTERNAL PROCEDURES
902	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
903
904	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
905	SUBROUTINE FUN_CHEM( T, Y, P )
906	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
907
908	USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO
909	USE KPP_ROOT_Global, ONLY: TIME, FIX, RCONST
910	USE KPP_ROOT_Function
911	USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST
912
913	INTEGER N
914	KPP_REAL T, Told
915	KPP_REAL Y(NVAR), P(NVAR)
916
917	Told = TIME
918	TIME = T
919	CALL Update_SUN()
920	CALL Update_RCONST()
921
922	CALL Fun( Y, FIX, RCONST, P )
923
924	TIME = Told
925	Nfun=Nfun+1
926
927	END SUBROUTINE FUN_CHEM
928
929
930	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
931	SUBROUTINE JAC_CHEM( T, Y, JV )
932	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
933
934	USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO
935	USE KPP_ROOT_Global, ONLY: TIME, FIX, RCONST
936	USE KPP_ROOT_Jacobian
937	USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST
938
939	INTEGER N
940	KPP_REAL T, Told
941	KPP_REAL Y(NVAR)
942	#ifdef FULL_ALGEBRA
943	KPP_REAL :: JS(LU_NONZERO), JV(NVAR,NVAR)
944	INTEGER :: i, j
945	#else
946	KPP_REAL :: JV(LU_NONZERO)
947	#endif
948
949	Told = TIME
950	TIME = T
951	CALL Update_SUN()
952	CALL Update_RCONST()
953
954	#ifdef FULL_ALGEBRA
955	CALL Jac_SP(Y, FIX, RCONST, JS)
956	DO j=1,NVAR
957	DO j=1,NVAR
958	JV(i,j) = 0.0D0
959	END DO
960	END DO
961	DO i=1,LU_NONZERO
962	JV(LU_IROW(i),LU_ICOL(i)) = JS(i)
963	END DO
964	#else
965	CALL Jac_SP(Y, FIX, RCONST, JV)
966	#endif
967	TIME = Told
968	Njac = Njac+1
969
970	END SUBROUTINE JAC_CHEM
971
972	END MODULE KPP_ROOT_Integrator

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

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/kpp_sdirk4.f90 @ 4086

Download in other formats: