Home

Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

rosenbrock_adj.f90 @ 3512

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

Line
1	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
2	! Discrete adjoints of Rosenbrock, !
3	! for several Rosenbrock methods: !
4	! * Ros2 !
5	! * Ros3 !
6	! * Ros4 !
7	! * Rodas3 !
8	! * Rodas4 !
9	! By default the code employs the KPP sparse linear algebra routines !
10	! Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK) !
11	! !
12	! (C) Adrian Sandu, August 2004 !
13	! Virginia Polytechnic Institute and State University !
14	! Contact: sandu@cs.vt.edu !
15	! Revised by Philipp Miehe and Adrian Sandu, May 2006 !
16	! This implementation is part of KPP - the Kinetic PreProcessor !
17	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
18
19	MODULE KPP_ROOT_Integrator
20
21	USE KPP_ROOT_Precision
22	USE KPP_ROOT_Parameters
23	USE KPP_ROOT_Global
24	USE KPP_ROOT_LinearAlgebra
25	USE KPP_ROOT_Rates
26	USE KPP_ROOT_Function
27	USE KPP_ROOT_Jacobian
28	USE KPP_ROOT_Hessian
29	USE KPP_ROOT_Util
30
31	IMPLICIT NONE
32	PUBLIC
33	SAVE
34
35	!~~~> Statistics on the work performed by the Rosenbrock method
36	INTEGER, PARAMETER :: Nfun=1, Njac=2, Nstp=3, Nacc=4, &
37	Nrej=5, Ndec=6, Nsol=7, Nsng=8, &
38	Ntexit=1, Nhexit=2, Nhnew = 3
39
40
41	CONTAINS ! Routines in the module KPP_ROOT_Integrator
42
43
44	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
45	SUBROUTINE INTEGRATE_ADJ( NADJ, Y, Lambda, TIN, TOUT, &
46	ATOL_adj, RTOL_adj, &
47	ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U )
48	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
49
50	IMPLICIT NONE
51
52	!~~~> Y - Concentrations
53	KPP_REAL :: Y(NVAR)
54	!~~~> NADJ - No. of cost functionals for which adjoints
55	! are evaluated simultaneously
56	! If single cost functional is considered (like in
57	! most applications) simply set NADJ = 1
58	INTEGER NADJ
59	!~~~> Lambda - Sensitivities w.r.t. concentrations
60	! Note: Lambda (1:NVAR,j) contains sensitivities of
61	! the j-th cost functional w.r.t. Y(1:NVAR), j=1...NADJ
62	KPP_REAL :: Lambda(NVAR,NADJ)
63	KPP_REAL, INTENT(IN) :: TIN ! TIN - Start Time
64	KPP_REAL, INTENT(IN) :: TOUT ! TOUT - End Time
65	!~~~> Tolerances for adjoint calculations
66	! (used only for full continuous adjoint)
67	KPP_REAL, INTENT(IN) :: ATOL_adj(NVAR,NADJ), RTOL_adj(NVAR,NADJ)
68	!~~~> Optional input parameters and statistics
69	INTEGER, INTENT(IN), OPTIONAL :: ICNTRL_U(20)
70	KPP_REAL, INTENT(IN), OPTIONAL :: RCNTRL_U(20)
71	INTEGER, INTENT(OUT), OPTIONAL :: ISTATUS_U(20)
72	KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20)
73
74	KPP_REAL :: RCNTRL(20), RSTATUS(20)
75	INTEGER :: ICNTRL(20), ISTATUS(20), IERR
76
77	INTEGER, SAVE :: Ntotal
78
79
80	ICNTRL(1:20) = 0
81	RCNTRL(1:20) = 0.0_dp
82	ISTATUS(1:20) = 0
83	RSTATUS(1:20) = 0.0_dp
84
85
86	!~~~> fine-tune the integrator:
87	! ICNTRL(1) = 0 ! 0 = non-autonomous, 1 = autonomous
88	! ICNTRL(2) = 1 ! 0 = scalar, 1 = vector tolerances
89	! RCNTRL(3) = STEPMIN ! starting step
90	! ICNTRL(3) = 5 ! choice of the method for forward integration
91	! ICNTRL(6) = 1 ! choice of the method for continuous adjoint
92	! ICNTRL(7) = 2 ! 1=none, 2=discrete, 3=full continuous, 4=simplified continuous adjoint
93	! ICNTRL(8) = 1 ! Save fwd LU factorization: 0 = don't save, 1 = save
94
95
96	! if optional parameters are given, and if they are >=0, then they overwrite default settings
97	IF (PRESENT(ICNTRL_U)) THEN
98	WHERE(ICNTRL_U(:) >= 0) ICNTRL(:) = ICNTRL_U(:)
99	END IF
100	IF (PRESENT(RCNTRL_U)) THEN
101	WHERE(RCNTRL_U(:) >= 0) RCNTRL(:) = RCNTRL_U(:)
102	END IF
103
104
105	CALL RosenbrockADJ(Y, NADJ, Lambda, &
106	TIN, TOUT, &
107	ATOL, RTOL, ATOL_adj, RTOL_adj, &
108	RCNTRL, ICNTRL, RSTATUS, ISTATUS, IERR)
109
110
111	!~~~> Debug option: show number of steps
112	! Ntotal = Ntotal + ISTATUS(Nstp)
113	! WRITE(6,777) ISTATUS(Nstp),Ntotal,VAR(ind_O3),VAR(ind_NO2)
114	!777 FORMAT('NSTEPS=',I6,' (',I6,') O3=',E24.14,' NO2=',E24.14)
115
116	IF (IERR < 0) THEN
117	print *,'RosenbrockADJ: Unsucessful step at T=', &
118	TIN,' (IERR=',IERR,')'
119	END IF
120
121	STEPMIN = RSTATUS(Nhexit)
122	! if optional parameters are given for output
123	! copy to them to return information
124	IF (PRESENT(ISTATUS_U)) ISTATUS_U(:) = ISTATUS(:)
125	IF (PRESENT(RSTATUS_U)) RSTATUS_U(:) = RSTATUS(:)
126
127	END SUBROUTINE INTEGRATE_ADJ
128
129
130	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
131	SUBROUTINE RosenbrockADJ( Y, NADJ, Lambda, &
132	Tstart, Tend, &
133	AbsTol, RelTol, AbsTol_adj, RelTol_adj, &
134	RCNTRL, ICNTRL, RSTATUS, ISTATUS, IERR)
135	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
136	!
137	! ADJ = Adjoint of the Tangent Linear Model of a Rosenbrock Method
138	!
139	! Solves the system y'=F(t,y) using a RosenbrockADJ method defined by:
140	!
141	! G = 1/(H*gamma(1)) - Jac(t0,Y0)
142	! T_i = t0 + Alpha(i)*H
143	! Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
144	! G * K_i = Fun( T_i, Y_i ) + \sum_{j=1}^S C(i,j)/H * K_j +
145	! gamma(i)*dF/dT(t0, Y0)
146	! Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
147	!
148	! For details on RosenbrockADJ methods and their implementation consult:
149	! E. Hairer and G. Wanner
150	! "Solving ODEs II. Stiff and differential-algebraic problems".
151	! Springer series in computational mathematics, Springer-Verlag, 1996.
152	! The codes contained in the book inspired this implementation.
153	!
154	! (C) Adrian Sandu, August 2004
155	! Virginia Polytechnic Institute and State University
156	! Contact: sandu@cs.vt.edu
157	! Revised by Philipp Miehe and Adrian Sandu, May 2006
158	! This implementation is part of KPP - the Kinetic PreProcessor
159	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
160	!
161	!~~~> INPUT ARGUMENTS:
162	!
163	!- Y(NVAR) = vector of initial conditions (at T=Tstart)
164	! NADJ -> dimension of linearized system,
165	! i.e. the number of sensitivity coefficients
166	!- Lambda(NVAR,NADJ) -> vector of initial sensitivity conditions (at T=Tstart)
167	!- [Tstart,Tend] = time range of integration
168	! (if Tstart>Tend the integration is performed backwards in time)
169	!- RelTol, AbsTol = user precribed accuracy
170	!- SUBROUTINE Fun( T, Y, Ydot ) = ODE function,
171	! returns Ydot = Y' = F(T,Y)
172	!- SUBROUTINE Jac( T, Y, Jcb ) = Jacobian of the ODE function,
173	! returns Jcb = dF/dY
174	!- ICNTRL(1:10) = integer inputs parameters
175	!- RCNTRL(1:10) = real inputs parameters
176	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
177	!
178	!~~~> OUTPUT ARGUMENTS:
179	!
180	!- Y(NVAR) -> vector of final states (at T->Tend)
181	!- Lambda(NVAR,NADJ) -> vector of final sensitivities (at T=Tend)
182	!- ICNTRL(11:20) -> integer output parameters
183	!- RCNTRL(11:20) -> real output parameters
184	!- IERR -> job status upon return
185	! - succes (positive value) or failure (negative value) -
186	! = 1 : Success
187	! = -1 : Improper value for maximal no of steps
188	! = -2 : Selected RosenbrockADJ method not implemented
189	! = -3 : Hmin/Hmax/Hstart must be positive
190	! = -4 : FacMin/FacMax/FacRej must be positive
191	! = -5 : Improper tolerance values
192	! = -6 : No of steps exceeds maximum bound
193	! = -7 : Step size too small
194	! = -8 : Matrix is repeatedly singular
195	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
196	!
197	!~~~> INPUT PARAMETERS:
198	!
199	! Note: For input parameters equal to zero the default values of the
200	! corresponding variables are used.
201	!
202	! ICNTRL(1) = 1: F = F(y) Independent of T (AUTONOMOUS)
203	! = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
204	!
205	! ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors
206	! = 1: AbsTol, RelTol are scalars
207	!
208	! ICNTRL(3) -> selection of a particular Rosenbrock method
209	! = 0 : default method is Rodas3
210	! = 1 : method is Ros2
211	! = 2 : method is Ros3
212	! = 3 : method is Ros4
213	! = 4 : method is Rodas3
214	! = 5: method is Rodas4
215	!
216	! ICNTRL(4) -> maximum number of integration steps
217	! For ICNTRL(4)=0) the default value of BUFSIZE is used
218	!
219	! ICNTRL(6) -> selection of a particular Rosenbrock method for the
220	! continuous adjoint integration - for cts adjoint it
221	! can be different than the forward method ICNTRL(3)
222	! Note 1: to avoid interpolation errors (which can be huge!)
223	! it is recommended to use only ICNTRL(7) = 2 or 4
224	! Note 2: the performance of the full continuous adjoint
225	! strongly depends on the forward solution accuracy Abs/RelTol
226	!
227	! ICNTRL(7) -> Type of adjoint algorithm
228	! = 0 : default is discrete adjoint ( of method ICNTRL(3) )
229	! = 1 : no adjoint
230	! = 2 : discrete adjoint ( of method ICNTRL(3) )
231	! = 3 : fully adaptive continuous adjoint ( with method ICNTRL(6) )
232	! = 4 : simplified continuous adjoint ( with method ICNTRL(6) )
233	!
234	! ICNTRL(8) -> checkpointing the LU factorization at each step:
235	! ICNTRL(8)=0 : do not save LU factorization (the default)
236	! ICNTRL(8)=1 : save LU factorization
237	! Note: if ICNTRL(7)=1 the LU factorization is not saved
238	!
239	!~~~> Real input parameters:
240	!
241	! RCNTRL(1) -> Hmin, lower bound for the integration step size
242	! It is strongly recommended to keep Hmin = ZERO
243	!
244	! RCNTRL(2) -> Hmax, upper bound for the integration step size
245	!
246	! RCNTRL(3) -> Hstart, starting value for the integration step size
247	!
248	! RCNTRL(4) -> FacMin, lower bound on step decrease factor (default=0.2)
249	!
250	! RCNTRL(5) -> FacMax, upper bound on step increase factor (default=6)
251	!
252	! RCNTRL(6) -> FacRej, step decrease factor after multiple rejections
253	! (default=0.1)
254	!
255	! RCNTRL(7) -> FacSafe, by which the new step is slightly smaller
256	! than the predicted value (default=0.9)
257	!
258	! RCNTRL(8) -> ThetaMin. If Newton convergence rate smaller
259	! than ThetaMin the Jacobian is not recomputed;
260	! (default=0.001)
261	!
262	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
263	!
264	!~~~> OUTPUT PARAMETERS:
265	!
266	! Note: each call to RosenbrockADJ adds the corrent no. of fcn calls
267	! to previous value of ISTATUS(1), and similar for the other params.
268	! Set ISTATUS(1:10) = 0 before call to avoid this accumulation.
269	!
270	! ISTATUS(1) = No. of function calls
271	! ISTATUS(2) = No. of jacobian calls
272	! ISTATUS(3) = No. of steps
273	! ISTATUS(4) = No. of accepted steps
274	! ISTATUS(5) = No. of rejected steps (except at the beginning)
275	! ISTATUS(6) = No. of LU decompositions
276	! ISTATUS(7) = No. of forward/backward substitutions
277	! ISTATUS(8) = No. of singular matrix decompositions
278	!
279	! RSTATUS(1) -> Texit, the time corresponding to the
280	! computed Y upon return
281	! RSTATUS(2) -> Hexit, last accepted step before exit
282	! For multiple restarts, use Hexit as Hstart in the following run
283	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
284
285	IMPLICIT NONE
286
287	!~~~> Arguments
288	KPP_REAL, INTENT(INOUT) :: Y(NVAR)
289	INTEGER, INTENT(IN) :: NADJ
290	KPP_REAL, INTENT(INOUT) :: Lambda(NVAR,NADJ)
291	KPP_REAL, INTENT(IN) :: Tstart,Tend
292	KPP_REAL, INTENT(IN) :: AbsTol(NVAR),RelTol(NVAR)
293	KPP_REAL, INTENT(IN) :: AbsTol_adj(NVAR,NADJ), RelTol_adj(NVAR,NADJ)
294	INTEGER, INTENT(IN) :: ICNTRL(20)
295	KPP_REAL, INTENT(IN) :: RCNTRL(20)
296	INTEGER, INTENT(INOUT) :: ISTATUS(20)
297	KPP_REAL, INTENT(INOUT) :: RSTATUS(20)
298	INTEGER, INTENT(OUT) :: IERR
299	!~~~> Parameters of the Rosenbrock method, up to 6 stages
300	INTEGER :: ros_S, rosMethod
301	INTEGER, PARAMETER :: RS2=1, RS3=2, RS4=3, RD3=4, RD4=5
302	KPP_REAL :: ros_A(15), ros_C(15), ros_M(6), ros_E(6), &
303	ros_Alpha(6), ros_Gamma(6), ros_ELO
304	LOGICAL :: ros_NewF(6)
305	CHARACTER(LEN=12) :: ros_Name
306	!~~~> Types of Adjoints Implemented
307	INTEGER, PARAMETER :: Adj_none = 1, Adj_discrete = 2, &
308	Adj_continuous = 3, Adj_simple_continuous = 4
309	!~~~> Checkpoints in memory
310	INTEGER, PARAMETER :: bufsize = 200000
311	INTEGER :: stack_ptr = 0 ! last written entry
312	KPP_REAL, DIMENSION(:), POINTER :: chk_H, chk_T
313	KPP_REAL, DIMENSION(:,:), POINTER :: chk_Y, chk_K, chk_J
314	KPP_REAL, DIMENSION(:,:), POINTER :: chk_dY, chk_d2Y
315	!~~~> Local variables
316	KPP_REAL :: Roundoff, FacMin, FacMax, FacRej, FacSafe
317	KPP_REAL :: Hmin, Hmax, Hstart
318	KPP_REAL :: Texit
319	INTEGER :: i, UplimTol, Max_no_steps
320	INTEGER :: AdjointType, CadjMethod
321	LOGICAL :: Autonomous, VectorTol, SaveLU
322	!~~~> Parameters
323	KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
324	KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5
325
326	!~~~> Initialize statistics
327	ISTATUS(1:20) = 0
328	RSTATUS(1:20) = ZERO
329
330	!~~~> Autonomous or time dependent ODE. Default is time dependent.
331	Autonomous = .NOT.(ICNTRL(1) == 0)
332
333	!~~~> For Scalar tolerances (ICNTRL(2).NE.0) the code uses AbsTol(1) and RelTol(1)
334	! For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:NVAR) and RelTol(1:NVAR)
335	IF (ICNTRL(2) == 0) THEN
336	VectorTol = .TRUE.
337	UplimTol = NVAR
338	ELSE
339	VectorTol = .FALSE.
340	UplimTol = 1
341	END IF
342
343	!~~~> Initialize the particular Rosenbrock method selected
344	SELECT CASE (ICNTRL(3))
345	CASE (1)
346	CALL Ros2
347	CASE (2)
348	CALL Ros3
349	CASE (3)
350	CALL Ros4
351	CASE (0,4)
352	CALL Rodas3
353	CASE (5)
354	CALL Rodas4
355	CASE DEFAULT
356	PRINT * , 'Unknown Rosenbrock method: ICNTRL(3)=',ICNTRL(3)
357	CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
358	RETURN
359	END SELECT
360
361	!~~~> The maximum number of steps admitted
362	IF (ICNTRL(4) == 0) THEN
363	Max_no_steps = bufsize - 1
364	ELSEIF (Max_no_steps > 0) THEN
365	Max_no_steps=ICNTRL(4)
366	ELSE
367	PRINT * ,'User-selected max no. of steps: ICNTRL(4)=',ICNTRL(4)
368	CALL ros_ErrorMsg(-1,Tstart,ZERO,IERR)
369	RETURN
370	END IF
371
372	!~~~> The particular Rosenbrock method chosen for integrating the cts adjoint
373	IF (ICNTRL(6) == 0) THEN
374	CadjMethod = 4
375	ELSEIF ( (ICNTRL(6) >= 1).AND.(ICNTRL(6) <= 5) ) THEN
376	CadjMethod = ICNTRL(6)
377	ELSE
378	PRINT * , 'Unknown CADJ Rosenbrock method: ICNTRL(6)=', CadjMethod
379	CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
380	RETURN
381	END IF
382
383	!~~~> Discrete or continuous adjoint formulation
384	IF ( ICNTRL(7) == 0 ) THEN
385	AdjointType = Adj_discrete
386	ELSEIF ( (ICNTRL(7) >= 1).AND.(ICNTRL(7) <= 4) ) THEN
387	AdjointType = ICNTRL(7)
388	ELSE
389	PRINT * , 'User-selected adjoint type: ICNTRL(7)=', AdjointType
390	CALL ros_ErrorMsg(-9,Tstart,ZERO,IERR)
391	RETURN
392	END IF
393
394	!~~~> Save or not the forward LU factorization
395	SaveLU = (ICNTRL(8) /= 0)
396
397
398	!~~~> Unit roundoff (1+Roundoff>1)
399	Roundoff = WLAMCH('E')
400
401	!~~~> Lower bound on the step size: (positive value)
402	IF (RCNTRL(1) == ZERO) THEN
403	Hmin = ZERO
404	ELSEIF (RCNTRL(1) > ZERO) THEN
405	Hmin = RCNTRL(1)
406	ELSE
407	PRINT * , 'User-selected Hmin: RCNTRL(1)=', RCNTRL(1)
408	CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
409	RETURN
410	END IF
411	!~~~> Upper bound on the step size: (positive value)
412	IF (RCNTRL(2) == ZERO) THEN
413	Hmax = ABS(Tend-Tstart)
414	ELSEIF (RCNTRL(2) > ZERO) THEN
415	Hmax = MIN(ABS(RCNTRL(2)),ABS(Tend-Tstart))
416	ELSE
417	PRINT * , 'User-selected Hmax: RCNTRL(2)=', RCNTRL(2)
418	CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
419	RETURN
420	END IF
421	!~~~> Starting step size: (positive value)
422	IF (RCNTRL(3) == ZERO) THEN
423	Hstart = MAX(Hmin,DeltaMin)
424	ELSEIF (RCNTRL(3) > ZERO) THEN
425	Hstart = MIN(ABS(RCNTRL(3)),ABS(Tend-Tstart))
426	ELSE
427	PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3)
428	CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
429	RETURN
430	END IF
431	!~~~> Step size can be changed s.t. FacMin < Hnew/Hold < FacMax
432	IF (RCNTRL(4) == ZERO) THEN
433	FacMin = 0.2d0
434	ELSEIF (RCNTRL(4) > ZERO) THEN
435	FacMin = RCNTRL(4)
436	ELSE
437	PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4)
438	CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
439	RETURN
440	END IF
441	IF (RCNTRL(5) == ZERO) THEN
442	FacMax = 6.0d0
443	ELSEIF (RCNTRL(5) > ZERO) THEN
444	FacMax = RCNTRL(5)
445	ELSE
446	PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5)
447	CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
448	RETURN
449	END IF
450	!~~~> FacRej: Factor to decrease step after 2 succesive rejections
451	IF (RCNTRL(6) == ZERO) THEN
452	FacRej = 0.1d0
453	ELSEIF (RCNTRL(6) > ZERO) THEN
454	FacRej = RCNTRL(6)
455	ELSE
456	PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6)
457	CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
458	RETURN
459	END IF
460	!~~~> FacSafe: Safety Factor in the computation of new step size
461	IF (RCNTRL(7) == ZERO) THEN
462	FacSafe = 0.9d0
463	ELSEIF (RCNTRL(7) > ZERO) THEN
464	FacSafe = RCNTRL(7)
465	ELSE
466	PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7)
467	CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
468	RETURN
469	END IF
470	!~~~> Check if tolerances are reasonable
471	DO i=1,UplimTol
472	IF ( (AbsTol(i) <= ZERO) .OR. (RelTol(i) <= 10.d0*Roundoff) &
473	.OR. (RelTol(i) >= 1.0d0) ) THEN
474	PRINT * , ' AbsTol(',i,') = ',AbsTol(i)
475	PRINT * , ' RelTol(',i,') = ',RelTol(i)
476	CALL ros_ErrorMsg(-5,Tstart,ZERO,IERR)
477	RETURN
478	END IF
479	END DO
480
481
482	!~~~> Allocate checkpoint space or open checkpoint files
483	IF (AdjointType == Adj_discrete) THEN
484	CALL ros_AllocateDBuffers( ros_S )
485	ELSEIF ( (AdjointType == Adj_continuous).OR. &
486	(AdjointType == Adj_simple_continuous) ) THEN
487	CALL ros_AllocateCBuffers
488	END IF
489
490	!~~~> CALL Forward Rosenbrock method
491	CALL ros_FwdInt(Y,Tstart,Tend,Texit, &
492	AbsTol, RelTol, &
493	! Error indicator
494	IERR)
495
496	PRINT*,'FORWARD STATISTICS'
497	PRINT*,'Step=',Nstp,' Acc=',Nacc, &
498	' Rej=',Nrej, ' Singular=',Nsng
499
500	!~~~> If Forward integration failed return
501	IF (IERR<0) RETURN
502
503	!~~~> Initialize the particular Rosenbrock method for continuous adjoint
504	IF ( (AdjointType == Adj_continuous).OR. &
505	(AdjointType == Adj_simple_continuous) ) THEN
506	SELECT CASE (CadjMethod)
507	CASE (1)
508	CALL Ros2
509	CASE (2)
510	CALL Ros3
511	CASE (3)
512	CALL Ros4
513	CASE (4)
514	CALL Rodas3
515	CASE (5)
516	CALL Rodas4
517	CASE DEFAULT
518	PRINT * , 'Unknown Rosenbrock method: ICNTRL(3)=', ICNTRL(3)
519	CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
520	RETURN
521	END SELECT
522	END IF
523
524	SELECT CASE (AdjointType)
525	CASE (Adj_discrete)
526	CALL ros_DadjInt ( &
527	NADJ, Lambda, &
528	Tstart, Tend, Texit, &
529	IERR )
530	CASE (Adj_continuous)
531	CALL ros_CadjInt ( &
532	NADJ, Lambda, &
533	Tend, Tstart, Texit, &
534	AbsTol_adj, RelTol_adj, &
535	IERR )
536	CASE (Adj_simple_continuous)
537	CALL ros_SimpleCadjInt ( &
538	NADJ, Lambda, &
539	Tstart, Tend, Texit, &
540	IERR )
541	END SELECT ! AdjointType
542
543	PRINT*,'ADJOINT STATISTICS'
544	PRINT*,'Step=',Nstp,' Acc=',Nacc, &
545	' Rej=',Nrej, ' Singular=',Nsng
546
547	!~~~> Free checkpoint space or close checkpoint files
548	IF (AdjointType == Adj_discrete) THEN
549	CALL ros_FreeDBuffers
550	ELSEIF ( (AdjointType == Adj_continuous) .OR. &
551	(AdjointType == Adj_simple_continuous) ) THEN
552	CALL ros_FreeCBuffers
553	END IF
554
555
556	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
557	CONTAINS ! Procedures internal to RosenbrockADJ
558	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
559
560
561	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
562	SUBROUTINE ros_AllocateDBuffers( S )
563	!~~~> Allocate buffer space for discrete adjoint
564	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
565	INTEGER :: i, S
566
567	ALLOCATE( chk_H(bufsize), STAT=i )
568	IF (i/=0) THEN
569	PRINT*,'Failed allocation of buffer H'; STOP
570	END IF
571	ALLOCATE( chk_T(bufsize), STAT=i )
572	IF (i/=0) THEN
573	PRINT*,'Failed allocation of buffer T'; STOP
574	END IF
575	ALLOCATE( chk_Y(NVAR*S,bufsize), STAT=i )
576	IF (i/=0) THEN
577	PRINT*,'Failed allocation of buffer Y'; STOP
578	END IF
579	ALLOCATE( chk_K(NVAR*S,bufsize), STAT=i )
580	IF (i/=0) THEN
581	PRINT*,'Failed allocation of buffer K'; STOP
582	END IF
583	IF (SaveLU) THEN
584	#ifdef FULL_ALGEBRA
585	ALLOCATE( chk_J(NVAR*NVAR,bufsize), STAT=i )
586	#else
587	ALLOCATE( chk_J(LU_NONZERO,bufsize), STAT=i )
588	#endif
589	IF (i/=0) THEN
590	PRINT*,'Failed allocation of buffer J'; STOP
591	END IF
592	END IF
593
594	END SUBROUTINE ros_AllocateDBuffers
595
596
597	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
598	SUBROUTINE ros_FreeDBuffers
599	!~~~> Dallocate buffer space for discrete adjoint
600	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
601	INTEGER :: i
602
603	DEALLOCATE( chk_H, STAT=i )
604	IF (i/=0) THEN
605	PRINT*,'Failed deallocation of buffer H'; STOP
606	END IF
607	DEALLOCATE( chk_T, STAT=i )
608	IF (i/=0) THEN
609	PRINT*,'Failed deallocation of buffer T'; STOP
610	END IF
611	DEALLOCATE( chk_Y, STAT=i )
612	IF (i/=0) THEN
613	PRINT*,'Failed deallocation of buffer Y'; STOP
614	END IF
615	DEALLOCATE( chk_K, STAT=i )
616	IF (i/=0) THEN
617	PRINT*,'Failed deallocation of buffer K'; STOP
618	END IF
619	IF (SaveLU) THEN
620	DEALLOCATE( chk_J, STAT=i )
621	IF (i/=0) THEN
622	PRINT*,'Failed deallocation of buffer J'; STOP
623	END IF
624	END IF
625
626	END SUBROUTINE ros_FreeDBuffers
627
628
629	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
630	SUBROUTINE ros_AllocateCBuffers
631	!~~~> Allocate buffer space for continuous adjoint
632	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
633	INTEGER :: i
634
635	ALLOCATE( chk_H(bufsize), STAT=i )
636	IF (i/=0) THEN
637	PRINT*,'Failed allocation of buffer H'; STOP
638	END IF
639	ALLOCATE( chk_T(bufsize), STAT=i )
640	IF (i/=0) THEN
641	PRINT*,'Failed allocation of buffer T'; STOP
642	END IF
643	ALLOCATE( chk_Y(NVAR,bufsize), STAT=i )
644	IF (i/=0) THEN
645	PRINT*,'Failed allocation of buffer Y'; STOP
646	END IF
647	ALLOCATE( chk_dY(NVAR,bufsize), STAT=i )
648	IF (i/=0) THEN
649	PRINT*,'Failed allocation of buffer dY'; STOP
650	END IF
651	ALLOCATE( chk_d2Y(NVAR,bufsize), STAT=i )
652	IF (i/=0) THEN
653	PRINT*,'Failed allocation of buffer d2Y'; STOP
654	END IF
655
656	END SUBROUTINE ros_AllocateCBuffers
657
658
659	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
660	SUBROUTINE ros_FreeCBuffers
661	!~~~> Dallocate buffer space for continuous adjoint
662	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
663	INTEGER :: i
664
665	DEALLOCATE( chk_H, STAT=i )
666	IF (i/=0) THEN
667	PRINT*,'Failed deallocation of buffer H'; STOP
668	END IF
669	DEALLOCATE( chk_T, STAT=i )
670	IF (i/=0) THEN
671	PRINT*,'Failed deallocation of buffer T'; STOP
672	END IF
673	DEALLOCATE( chk_Y, STAT=i )
674	IF (i/=0) THEN
675	PRINT*,'Failed deallocation of buffer Y'; STOP
676	END IF
677	DEALLOCATE( chk_dY, STAT=i )
678	IF (i/=0) THEN
679	PRINT*,'Failed deallocation of buffer dY'; STOP
680	END IF
681	DEALLOCATE( chk_d2Y, STAT=i )
682	IF (i/=0) THEN
683	PRINT*,'Failed deallocation of buffer d2Y'; STOP
684	END IF
685
686	END SUBROUTINE ros_FreeCBuffers
687
688	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
689	SUBROUTINE ros_DPush( S, T, H, Ystage, K, E, P )
690	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
691	!~~~> Saves the next trajectory snapshot for discrete adjoints
692	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
693	INTEGER :: S ! no of stages
694	KPP_REAL :: T, H, Ystage(NVARS), K(NVARS)
695	INTEGER :: P(NVAR)
696	#ifdef FULL_ALGEBRA
697	KPP_REAL :: E(NVAR,NVAR)
698	#else
699	KPP_REAL :: E(LU_NONZERO)
700	#endif
701
702	stack_ptr = stack_ptr + 1
703	IF ( stack_ptr > bufsize ) THEN
704	PRINT*,'Push failed: buffer overflow'
705	STOP
706	END IF
707	chk_H( stack_ptr ) = H
708	chk_T( stack_ptr ) = T
709	!CALL WCOPY(NVAR*S,Ystage,1,chk_Y(1,stack_ptr),1)
710	!CALL WCOPY(NVAR*S,K,1,chk_K(1,stack_ptr),1)
711	chk_Y(1:NVARS,stack_ptr) = Ystage(1:NVARS)
712	chk_K(1:NVARS,stack_ptr) = K(1:NVARS)
713	IF (SaveLU) THEN
714	#ifdef FULL_ALGEBRA
715	chk_J(1:NVAR,1:NVAR,stack_ptr) = E(1:NVAR,1:NVAR)
716	chk_P(1:NVAR,stack_ptr) = P(1:NVAR)
717	#else
718	chk_J(1:LU_NONZERO,stack_ptr) = E(1:LU_NONZERO)
719	#endif
720	END IF
721
722	END SUBROUTINE ros_DPush
723
724
725	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
726	SUBROUTINE ros_DPop( S, T, H, Ystage, K, E, P )
727	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
728	!~~~> Retrieves the next trajectory snapshot for discrete adjoints
729	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
730
731	INTEGER :: S ! no of stages
732	KPP_REAL :: T, H, Ystage(NVARS), K(NVARS)
733	INTEGER :: P(NVAR)
734	#ifdef FULL_ALGEBRA
735	KPP_REAL :: E(NVAR,NVAR)
736	#else
737	KPP_REAL :: E(LU_NONZERO)
738	#endif
739
740	IF ( stack_ptr <= 0 ) THEN
741	PRINT*,'Pop failed: empty buffer'
742	STOP
743	END IF
744	H = chk_H( stack_ptr )
745	T = chk_T( stack_ptr )
746	!CALL WCOPY(NVAR*S,chk_Y(1,stack_ptr),1,Ystage,1)
747	!CALL WCOPY(NVAR*S,chk_K(1,stack_ptr),1,K,1)
748	Ystage(1:NVARS) = chk_Y(1:NVARS,stack_ptr)
749	K(1:NVARS) = chk_K(1:NVARS,stack_ptr)
750	!CALL WCOPY(LU_NONZERO,chk_J(1,stack_ptr),1,Jcb,1)
751	IF (SaveLU) THEN
752	#ifdef FULL_ALGEBRA
753	E(1:NVAR,1:NVAR) = chk_J(1:NVAR,1:NVAR,stack_ptr)
754	P(1:NVAR) = chk_P(1:NVAR,stack_ptr)
755	#else
756	E(1:LU_NONZERO) = chk_J(1:LU_NONZERO,stack_ptr)
757	#endif
758	END IF
759
760	stack_ptr = stack_ptr - 1
761
762	END SUBROUTINE ros_DPop
763
764	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
765	SUBROUTINE ros_CPush( T, H, Y, dY, d2Y )
766	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
767	!~~~> Saves the next trajectory snapshot for discrete adjoints
768	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
769
770	KPP_REAL :: T, H, Y(NVAR), dY(NVAR), d2Y(NVAR)
771
772	stack_ptr = stack_ptr + 1
773	IF ( stack_ptr > bufsize ) THEN
774	PRINT*,'Push failed: buffer overflow'
775	STOP
776	END IF
777	chk_H( stack_ptr ) = H
778	chk_T( stack_ptr ) = T
779	!CALL WCOPY(NVAR,Y,1,chk_Y(1,stack_ptr),1)
780	!CALL WCOPY(NVAR,dY,1,chk_dY(1,stack_ptr),1)
781	!CALL WCOPY(NVAR,d2Y,1,chk_d2Y(1,stack_ptr),1)
782	chk_Y(1:NVAR,stack_ptr) = Y(1:NVAR)
783	chk_dY(1:NVAR,stack_ptr) = dY(1:NVAR)
784	chk_d2Y(1:NVAR,stack_ptr) = d2Y(1:NVAR)
785	END SUBROUTINE ros_CPush
786
787
788	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
789	SUBROUTINE ros_CPop( T, H, Y, dY, d2Y )
790	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
791	!~~~> Retrieves the next trajectory snapshot for discrete adjoints
792	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
793
794	KPP_REAL :: T, H, Y(NVAR), dY(NVAR), d2Y(NVAR)
795
796	IF ( stack_ptr <= 0 ) THEN
797	PRINT*,'Pop failed: empty buffer'
798	STOP
799	END IF
800	H = chk_H( stack_ptr )
801	T = chk_T( stack_ptr )
802	!CALL WCOPY(NVAR,chk_Y(1,stack_ptr),1,Y,1)
803	!CALL WCOPY(NVAR,chk_dY(1,stack_ptr),1,dY,1)
804	!CALL WCOPY(NVAR,chk_d2Y(1,stack_ptr),1,d2Y,1)
805	Y(1:NVAR) = chk_Y(1:NVAR,stack_ptr)
806	dY(1:NVAR) = chk_dY(1:NVAR,stack_ptr)
807	d2Y(1:NVAR) = chk_d2Y(1:NVAR,stack_ptr)
808
809	stack_ptr = stack_ptr - 1
810
811	END SUBROUTINE ros_CPop
812
813
814
815	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
816	SUBROUTINE ros_ErrorMsg(Code,T,H,IERR)
817	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
818	! Handles all error messages
819	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
820
821	KPP_REAL, INTENT(IN) :: T, H
822	INTEGER, INTENT(IN) :: Code
823	INTEGER, INTENT(OUT) :: IERR
824
825	IERR = Code
826	PRINT * , &
827	'Forced exit from RosenbrockADJ due to the following error:'
828
829	SELECT CASE (Code)
830	CASE (-1)
831	PRINT * , '--> Improper value for maximal no of steps'
832	CASE (-2)
833	PRINT * , '--> Selected RosenbrockADJ method not implemented'
834	CASE (-3)
835	PRINT * , '--> Hmin/Hmax/Hstart must be positive'
836	CASE (-4)
837	PRINT * , '--> FacMin/FacMax/FacRej must be positive'
838	CASE (-5)
839	PRINT * , '--> Improper tolerance values'
840	CASE (-6)
841	PRINT * , '--> No of steps exceeds maximum buffer bound'
842	CASE (-7)
843	PRINT * , '--> Step size too small: T + 10*H = T', &
844	' or H < Roundoff'
845	CASE (-8)
846	PRINT * , '--> Matrix is repeatedly singular'
847	CASE (-9)
848	PRINT * , '--> Improper type of adjoint selected'
849	CASE DEFAULT
850	PRINT *, 'Unknown Error code: ', Code
851	END SELECT
852
853	PRINT *, "T=", T, "and H=", H
854
855	END SUBROUTINE ros_ErrorMsg
856
857
858
859	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
860	SUBROUTINE ros_FwdInt (Y, &
861	Tstart, Tend, T, &
862	AbsTol, RelTol, &
863	!~~~> Error indicator
864	IERR )
865	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
866	! Template for the implementation of a generic RosenbrockADJ method
867	! defined by ros_S (no of stages)
868	! and its coefficients ros_{A,C,M,E,Alpha,Gamma}
869	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
870
871	IMPLICIT NONE
872
873	!~~~> Input: the initial condition at Tstart; Output: the solution at T
874	KPP_REAL, INTENT(INOUT) :: Y(NVAR)
875	!~~~> Input: integration interval
876	KPP_REAL, INTENT(IN) :: Tstart,Tend
877	!~~~> Output: time at which the solution is returned (T=Tend if success)
878	KPP_REAL, INTENT(OUT) :: T
879	!~~~> Input: tolerances
880	KPP_REAL, INTENT(IN) :: AbsTol(NVAR), RelTol(NVAR)
881	!~~~> Output: Error indicator
882	INTEGER, INTENT(OUT) :: IERR
883	! ~~~~ Local variables
884	KPP_REAL :: Ynew(NVAR), Fcn0(NVAR), Fcn(NVAR)
885	KPP_REAL :: K(NVAR*ros_S), dFdT(NVAR)
886	KPP_REAL, DIMENSION(:), POINTER :: Ystage
887	#ifdef FULL_ALGEBRA
888	KPP_REAL :: Jac0(NVAR,NVAR), Ghimj(NVAR,NVAR)
889	#else
890	KPP_REAL :: Jac0(LU_NONZERO), Ghimj(LU_NONZERO)
891	#endif
892	KPP_REAL :: H, Hnew, HC, HG, Fac, Tau
893	KPP_REAL :: Err, Yerr(NVAR)
894	INTEGER :: Pivot(NVAR), Direction, ioffset, i, j, istage
895	LOGICAL :: RejectLastH, RejectMoreH, Singular
896	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
897
898	!~~~> Allocate stage vector buffer if needed
899	IF (AdjointType == Adj_discrete) THEN
900	ALLOCATE(Ystage(NVAR*ros_S), STAT=i)
901	! Uninitialized Ystage may lead to NaN on some compilers
902	Ystage = 0.0d0
903	IF (i/=0) THEN
904	PRINT*,'Allocation of Ystage failed'
905	STOP
906	END IF
907	END IF
908
909	!~~~> Initial preparations
910	T = Tstart
911	RSTATUS(Nhexit) = ZERO
912	H = MIN( MAX(ABS(Hmin),ABS(Hstart)) , ABS(Hmax) )
913	IF (ABS(H) <= 10.0_dp*Roundoff) H = DeltaMin
914
915	IF (Tend >= Tstart) THEN
916	Direction = +1
917	ELSE
918	Direction = -1
919	END IF
920	H = Direction*H
921
922	RejectLastH=.FALSE.
923	RejectMoreH=.FALSE.
924
925	!~~~> Time loop begins below
926
927	TimeLoop: DO WHILE ( (Direction > 0).AND.((T-Tend)+Roundoff <= ZERO) &
928	.OR. (Direction < 0).AND.((Tend-T)+Roundoff <= ZERO) )
929
930	IF ( ISTATUS(Nstp) > Max_no_steps ) THEN ! Too many steps
931	CALL ros_ErrorMsg(-6,T,H,IERR)
932	RETURN
933	END IF
934	IF ( ((T+0.1d0*H) == T).OR.(H <= Roundoff) ) THEN ! Step size too small
935	CALL ros_ErrorMsg(-7,T,H,IERR)
936	RETURN
937	END IF
938
939	!~~~> Limit H if necessary to avoid going beyond Tend
940	RSTATUS(Nhexit) = H
941	H = MIN(H,ABS(Tend-T))
942
943	!~~~> Compute the function at current time
944	CALL FunTemplate(T,Y,Fcn0)
945	ISTATUS(Nfun) = ISTATUS(Nfun) + 1
946
947	!~~~> Compute the function derivative with respect to T
948	IF (.NOT.Autonomous) THEN
949	CALL ros_FunTimeDerivative ( T, Roundoff, Y, &
950	Fcn0, dFdT )
951	END IF
952
953	!~~~> Compute the Jacobian at current time
954	CALL JacTemplate(T,Y,Jac0)
955	ISTATUS(Njac) = ISTATUS(Njac) + 1
956
957	!~~~> Repeat step calculation until current step accepted
958	UntilAccepted: DO
959
960	CALL ros_PrepareMatrix(H,Direction,ros_Gamma(1), &
961	Jac0,Ghimj,Pivot,Singular)
962	IF (Singular) THEN ! More than 5 consecutive failed decompositions
963	CALL ros_ErrorMsg(-8,T,H,IERR)
964	RETURN
965	END IF
966
967	!~~~> Compute the stages
968	Stage: DO istage = 1, ros_S
969
970	! Current istage offset. Current istage vector is K(ioffset+1:ioffset+NVAR)
971	ioffset = NVAR*(istage-1)
972
973	! For the 1st istage the function has been computed previously
974	IF ( istage == 1 ) THEN
975	CALL WCOPY(NVAR,Fcn0,1,Fcn,1)
976	IF (AdjointType == Adj_discrete) THEN ! Save stage solution
977	! CALL WCOPY(NVAR,Y,1,Ystage(1),1)
978	Ystage(1:NVAR) = Y(1:NVAR)
979	CALL WCOPY(NVAR,Y,1,Ynew,1)
980	END IF
981	! istage>1 and a new function evaluation is needed at the current istage
982	ELSEIF ( ros_NewF(istage) ) THEN
983	CALL WCOPY(NVAR,Y,1,Ynew,1)
984	DO j = 1, istage-1
985	CALL WAXPY(NVAR,ros_A((istage-1)*(istage-2)/2+j), &
986	K(NVAR*(j-1)+1),1,Ynew,1)
987	END DO
988	Tau = T + ros_Alpha(istage)DirectionH
989	CALL FunTemplate(Tau,Ynew,Fcn)
990	ISTATUS(Nfun) = ISTATUS(Nfun) + 1
991	END IF ! if istage == 1 elseif ros_NewF(istage)
992	! save stage solution every time even if ynew is not updated
993	IF ( ( istage > 1 ).AND.(AdjointType == Adj_discrete) ) THEN
994	! CALL WCOPY(NVAR,Ynew,1,Ystage(ioffset+1),1)
995	Ystage(ioffset+1:ioffset+NVAR) = Ynew(1:NVAR)
996	END IF
997	CALL WCOPY(NVAR,Fcn,1,K(ioffset+1),1)
998	DO j = 1, istage-1
999	HC = ros_C((istage-1)(istage-2)/2+j)/(DirectionH)
1000	CALL WAXPY(NVAR,HC,K(NVAR*(j-1)+1),1,K(ioffset+1),1)
1001	END DO
1002	IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN
1003	HG = DirectionHros_Gamma(istage)
1004	CALL WAXPY(NVAR,HG,dFdT,1,K(ioffset+1),1)
1005	END IF
1006	CALL ros_Solve('N', Ghimj, Pivot, K(ioffset+1))
1007
1008	END DO Stage
1009
1010
1011	!~~~> Compute the new solution
1012	CALL WCOPY(NVAR,Y,1,Ynew,1)
1013	DO j=1,ros_S
1014	CALL WAXPY(NVAR,ros_M(j),K(NVAR*(j-1)+1),1,Ynew,1)
1015	END DO
1016
1017	!~~~> Compute the error estimation
1018	CALL WSCAL(NVAR,ZERO,Yerr,1)
1019	DO j=1,ros_S
1020	CALL WAXPY(NVAR,ros_E(j),K(NVAR*(j-1)+1),1,Yerr,1)
1021	END DO
1022	Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol )
1023
1024	!~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax
1025	Fac = MIN(FacMax,MAX(FacMin,FacSafe/Err**(ONE/ros_ELO)))
1026	Hnew = H*Fac
1027
1028	!~~~> Check the error magnitude and adjust step size
1029	ISTATUS(Nstp) = ISTATUS(Nstp) + 1
1030	IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN !~~~> Accept step
1031	ISTATUS(Nacc) = ISTATUS(Nacc) + 1
1032	IF (AdjointType == Adj_discrete) THEN ! Save current state
1033	CALL ros_DPush( ros_S, T, H, Ystage, K, Ghimj, Pivot )
1034	ELSEIF ( (AdjointType == Adj_continuous) .OR. &
1035	(AdjointType == Adj_simple_continuous) ) THEN
1036	#ifdef FULL_ALGEBRA
1037	K = MATMUL(Jac0,Fcn0)
1038	#else
1039	CALL Jac_SP_Vec( Jac0, Fcn0, K(1) )
1040	#endif
1041	IF (.NOT. Autonomous) THEN
1042	CALL WAXPY(NVAR,ONE,dFdT,1,K(1),1)
1043	END IF
1044	CALL ros_CPush( T, H, Y, Fcn0, K(1) )
1045	END IF
1046	CALL WCOPY(NVAR,Ynew,1,Y,1)
1047	T = T + Direction*H
1048	Hnew = MAX(Hmin,MIN(Hnew,Hmax))
1049	IF (RejectLastH) THEN ! No step size increase after a rejected step
1050	Hnew = MIN(Hnew,H)
1051	END IF
1052	RSTATUS(Nhexit) = H
1053	RSTATUS(Nhnew) = Hnew
1054	RSTATUS(Ntexit) = T
1055	RejectLastH = .FALSE.
1056	RejectMoreH = .FALSE.
1057	H = Hnew
1058	EXIT UntilAccepted ! EXIT THE LOOP: WHILE STEP NOT ACCEPTED
1059	ELSE !~~~> Reject step
1060	IF (RejectMoreH) THEN
1061	Hnew = H*FacRej
1062	END IF
1063	RejectMoreH = RejectLastH
1064	RejectLastH = .TRUE.
1065	H = Hnew
1066	IF (ISTATUS(Nacc) >= 1) THEN
1067	ISTATUS(Nrej) = ISTATUS(Nrej) + 1
1068	END IF
1069	END IF ! Err <= 1
1070
1071	END DO UntilAccepted
1072
1073	END DO TimeLoop
1074
1075	!~~~> Save last state: only needed for continuous adjoint
1076	IF ( (AdjointType == Adj_continuous) .OR. &
1077	(AdjointType == Adj_simple_continuous) ) THEN
1078	CALL FunTemplate(T,Y,Fcn0)
1079	ISTATUS(Nfun) = ISTATUS(Nfun) + 1
1080	CALL JacTemplate(T,Y,Jac0)
1081	ISTATUS(Njac) = ISTATUS(Njac) + 1
1082	#ifdef FULL_ALGEBRA
1083	K = MATMUL(Jac0,Fcn0)
1084	#else
1085	CALL Jac_SP_Vec( Jac0, Fcn0, K(1) )
1086	#endif
1087	IF (.NOT. Autonomous) THEN
1088	CALL ros_FunTimeDerivative ( T, Roundoff, Y, &
1089	Fcn0, dFdT )
1090	CALL WAXPY(NVAR,ONE,dFdT,1,K(1),1)
1091	END IF
1092	CALL ros_CPush( T, H, Y, Fcn0, K(1) )
1093	!~~~> Deallocate stage buffer: only needed for discrete adjoint
1094	ELSEIF (AdjointType == Adj_discrete) THEN
1095	DEALLOCATE(Ystage, STAT=i)
1096	IF (i/=0) THEN
1097	PRINT*,'Deallocation of Ystage failed'
1098	STOP
1099	END IF
1100	END IF
1101
1102	!~~~> Succesful exit
1103	IERR = 1 !~~~> The integration was successful
1104
1105	END SUBROUTINE ros_FwdInt
1106
1107
1108	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1109	SUBROUTINE ros_DadjInt ( &
1110	NADJ, Lambda, &
1111	Tstart, Tend, T, &
1112	!~~~> Error indicator
1113	IERR )
1114	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1115	! Template for the implementation of a generic RosenbrockSOA method
1116	! defined by ros_S (no of stages)
1117	! and its coefficients ros_{A,C,M,E,Alpha,Gamma}
1118	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1119
1120	IMPLICIT NONE
1121
1122	!~~~> Input: the initial condition at Tstart; Output: the solution at T
1123	INTEGER, INTENT(IN) :: NADJ
1124	!~~~> First order adjoint
1125	KPP_REAL, INTENT(INOUT) :: Lambda(NVAR,NADJ)
1126	!!~~~> Input: integration interval
1127	KPP_REAL, INTENT(IN) :: Tstart,Tend
1128	!~~~> Output: time at which the solution is returned (T=Tend if success)
1129	KPP_REAL, INTENT(OUT) :: T
1130	!~~~> Output: Error indicator
1131	INTEGER, INTENT(OUT) :: IERR
1132	! ~~~~ Local variables
1133	KPP_REAL :: Ystage(NVARros_S), K(NVARros_S)
1134	KPP_REAL :: U(NVARros_S,NADJ), V(NVARros_S,NADJ)
1135	#ifdef FULL_ALGEBRA
1136	KPP_REAL, DIMENSION(NVAR,NVAR) :: Jac, dJdT, Ghimj
1137	#else
1138	KPP_REAL, DIMENSION(LU_NONZERO) :: Jac, dJdT, Ghimj
1139	#endif
1140	KPP_REAL :: Hes0(NHESS)
1141	KPP_REAL :: Tmp(NVAR), Tmp2(NVAR)
1142	KPP_REAL :: H, HC, HA, Tau
1143	INTEGER :: Pivot(NVAR), Direction
1144	INTEGER :: i, j, m, istage, istart, jstart
1145	!~~~> Local parameters
1146	KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
1147	KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5
1148	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1149
1150
1151
1152	IF (Tend >= Tstart) THEN
1153	Direction = +1
1154	ELSE
1155	Direction = -1
1156	END IF
1157
1158	!~~~> Time loop begins below
1159	TimeLoop: DO WHILE ( stack_ptr > 0 )
1160
1161	!~~~> Recover checkpoints for stage values and vectors
1162	CALL ros_DPop( ros_S, T, H, Ystage, K, Ghimj, Pivot )
1163
1164	! ISTATUS(Nstp) = ISTATUS(Nstp) + 1
1165
1166	!~~~> Compute LU decomposition
1167	IF (.NOT.SaveLU) THEN
1168	CALL JacTemplate(T,Ystage(1),Ghimj)
1169	ISTATUS(Njac) = ISTATUS(Njac) + 1
1170	Tau = ONE/(DirectionHros_Gamma(1))
1171	#ifdef FULL_ALGEBRA
1172	Ghimj(1:NVAR,1:NVAR) = -Ghimj(1:NVAR,1:NVAR)
1173	DO i=1,NVAR
1174	Ghimj(i,i) = Ghimj(i,i)+Tau
1175	END DO
1176	#else
1177	CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1)
1178	DO i=1,NVAR
1179	Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+Tau
1180	END DO
1181	#endif
1182	CALL ros_Decomp( Ghimj, Pivot, j )
1183	END IF
1184
1185	!~~~> Compute Hessian at the beginning of the interval
1186	CALL HessTemplate(T,Ystage(1),Hes0)
1187
1188	!~~~> Compute the stages
1189	Stage: DO istage = ros_S, 1, -1
1190
1191	!~~~> Current istage first entry
1192	istart = NVAR*(istage-1) + 1
1193
1194	!~~~> Compute U
1195	DO m = 1,NADJ
1196	CALL WCOPY(NVAR,Lambda(1,m),1,U(istart,m),1)
1197	CALL WSCAL(NVAR,ros_M(istage),U(istart,m),1)
1198	END DO ! m=1:NADJ
1199	DO j = istage+1, ros_S
1200	jstart = NVAR*(j-1) + 1
1201	HA = ros_A((j-1)*(j-2)/2+istage)
1202	HC = ros_C((j-1)(j-2)/2+istage)/(DirectionH)
1203	DO m = 1,NADJ
1204	CALL WAXPY(NVAR,HA,V(jstart,m),1,U(istart,m),1)
1205	CALL WAXPY(NVAR,HC,U(jstart,m),1,U(istart,m),1)
1206	END DO ! m=1:NADJ
1207	END DO
1208	DO m = 1,NADJ
1209	CALL ros_Solve('T', Ghimj, Pivot, U(istart,m))
1210	END DO ! m=1:NADJ
1211	!~~~> Compute V
1212	Tau = T + ros_Alpha(istage)DirectionH
1213	CALL JacTemplate(Tau,Ystage(istart),Jac)
1214	ISTATUS(Njac) = ISTATUS(Njac) + 1
1215	DO m = 1,NADJ
1216	#ifdef FULL_ALGEBRA
1217	V(istart:istart+NVAR-1,m) = MATMUL(TRANSPOSE(Jac),U(istart:istart+NVAR-1,m))
1218	#else
1219	CALL JacTR_SP_Vec(Jac,U(istart,m),V(istart,m))
1220	#endif
1221	END DO ! m=1:NADJ
1222
1223	END DO Stage
1224
1225	IF (.NOT.Autonomous) THEN
1226	!~~~> Compute the Jacobian derivative with respect to T.
1227	! Last "Jac" computed for stage 1
1228	CALL ros_JacTimeDerivative ( T, Roundoff, Ystage(1), Jac, dJdT )
1229	END IF
1230
1231	!~~~> Compute the new solution
1232
1233	!~~~> Compute Lambda
1234	DO istage=1,ros_S
1235	istart = NVAR*(istage-1) + 1
1236	DO m = 1,NADJ
1237	! Add V_i
1238	CALL WAXPY(NVAR,ONE,V(istart,m),1,Lambda(1,m),1)
1239	! Add (H0xK_i)^T * U_i
1240	CALL HessTR_Vec ( Hes0, U(istart,m), K(istart), Tmp )
1241	CALL WAXPY(NVAR,ONE,Tmp,1,Lambda(1,m),1)
1242	END DO ! m=1:NADJ
1243	END DO
1244	! Add H * dJac_dT_0^T * \sum(gamma_i U_i)
1245	! Tmp holds sum gamma_i U_i
1246	IF (.NOT.Autonomous) THEN
1247	DO m = 1,NADJ
1248	Tmp(1:NVAR) = ZERO
1249	DO istage = 1, ros_S
1250	istart = NVAR*(istage-1) + 1
1251	CALL WAXPY(NVAR,ros_Gamma(istage),U(istart,m),1,Tmp,1)
1252	END DO
1253	#ifdef FULL_ALGEBRA
1254	Tmp2 = MATMUL(TRANSPOSE(dJdT),Tmp)
1255	#else
1256	CALL JacTR_SP_Vec(dJdT,Tmp,Tmp2)
1257	#endif
1258	CALL WAXPY(NVAR,H,Tmp2,1,Lambda(1,m),1)
1259	END DO ! m=1:NADJ
1260	END IF ! .NOT.Autonomous
1261
1262
1263	END DO TimeLoop
1264
1265	!~~~> Save last state
1266
1267	!~~~> Succesful exit
1268	IERR = 1 !~~~> The integration was successful
1269
1270	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1271	END SUBROUTINE ros_DadjInt
1272	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1273
1274
1275
1276	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1277	SUBROUTINE ros_CadjInt ( &
1278	NADJ, Y, &
1279	Tstart, Tend, T, &
1280	AbsTol_adj, RelTol_adj, &
1281	!~~~> Error indicator
1282	IERR )
1283	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1284	! Template for the implementation of a generic RosenbrockADJ method
1285	! defined by ros_S (no of stages)
1286	! and its coefficients ros_{A,C,M,E,Alpha,Gamma}
1287	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1288
1289	IMPLICIT NONE
1290
1291	!~~~> Input: the initial condition at Tstart; Output: the solution at T
1292	INTEGER, INTENT(IN) :: NADJ
1293	KPP_REAL, INTENT(INOUT) :: Y(NVAR,NADJ)
1294	!~~~> Input: integration interval
1295	KPP_REAL, INTENT(IN) :: Tstart,Tend
1296	!~~~> Input: adjoint tolerances
1297	KPP_REAL, INTENT(IN) :: AbsTol_adj(NVAR,NADJ), RelTol_adj(NVAR,NADJ)
1298	!~~~> Output: time at which the solution is returned (T=Tend if success)
1299	KPP_REAL, INTENT(OUT) :: T
1300	!~~~> Output: Error indicator
1301	INTEGER, INTENT(OUT) :: IERR
1302	! ~~~~ Local variables
1303	KPP_REAL :: Y0(NVAR)
1304	KPP_REAL :: Ynew(NVAR,NADJ), Fcn0(NVAR,NADJ), Fcn(NVAR,NADJ)
1305	KPP_REAL :: K(NVAR*ros_S,NADJ), dFdT(NVAR,NADJ)
1306	#ifdef FULL_ALGEBRA
1307	KPP_REAL, DIMENSION(NVAR,NVAR) :: Jac0, Ghimj, Jac, dJdT
1308	#else
1309	KPP_REAL, DIMENSION(LU_NONZERO) :: Jac0, Ghimj, Jac, dJdT
1310	#endif
1311	KPP_REAL :: H, Hnew, HC, HG, Fac, Tau
1312	KPP_REAL :: Err, Yerr(NVAR,NADJ)
1313	INTEGER :: Pivot(NVAR), Direction, ioffset, j, istage, iadj
1314	LOGICAL :: RejectLastH, RejectMoreH, Singular
1315	!~~~> Local parameters
1316	KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
1317	KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5
1318	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1319
1320
1321	!~~~> Initial preparations
1322	T = Tstart
1323	RSTATUS(Nhexit) = 0.0_dp
1324	H = MIN( MAX(ABS(Hmin),ABS(Hstart)) , ABS(Hmax) )
1325	IF (ABS(H) <= 10.0_dp*Roundoff) H = DeltaMin
1326
1327	IF (Tend >= Tstart) THEN
1328	Direction = +1
1329	ELSE
1330	Direction = -1
1331	END IF
1332	H = Direction*H
1333
1334	RejectLastH=.FALSE.
1335	RejectMoreH=.FALSE.
1336
1337	!~~~> Time loop begins below
1338
1339	TimeLoop: DO WHILE ( (Direction > 0).AND.((T-Tend)+Roundoff <= ZERO) &
1340	.OR. (Direction < 0).AND.((Tend-T)+Roundoff <= ZERO) )
1341
1342	IF ( ISTATUS(Nstp) > Max_no_steps ) THEN ! Too many steps
1343	CALL ros_ErrorMsg(-6,T,H,IERR)
1344	RETURN
1345	END IF
1346	IF ( ((T+0.1d0*H) == T).OR.(H <= Roundoff) ) THEN ! Step size too small
1347	CALL ros_ErrorMsg(-7,T,H,IERR)
1348	RETURN
1349	END IF
1350
1351	!~~~> Limit H if necessary to avoid going beyond Tend
1352	RSTATUS(Nhexit) = H
1353	H = MIN(H,ABS(Tend-T))
1354
1355	!~~~> Interpolate forward solution
1356	CALL ros_cadj_Y( T, Y0 )
1357	!~~~> Compute the Jacobian at current time
1358	CALL JacTemplate(T, Y0, Jac0)
1359	ISTATUS(Njac) = ISTATUS(Njac) + 1
1360
1361	!~~~> Compute the function derivative with respect to T
1362	IF (.NOT.Autonomous) THEN
1363	CALL ros_JacTimeDerivative ( T, Roundoff, Y0, &
1364	Jac0, dJdT )
1365	DO iadj = 1, NADJ
1366	#ifdef FULL_ALGEBRA
1367	dFdT(1:NVAR,iadj) = MATMUL(TRANSPOSE(dJdT),Y(1:NVAR,iadj))
1368	#else
1369	CALL JacTR_SP_Vec(dJdT,Y(1,iadj),dFdT(1,iadj))
1370	#endif
1371	CALL WSCAL(NVAR,(-ONE),dFdT(1,iadj),1)
1372	END DO
1373	END IF
1374
1375	!~~~> Ydot = -J^T*Y
1376	#ifdef FULL_ALGEBRA
1377	Jac0(1:NVAR,1:NVAR) = -Jac0(1:NVAR,1:NVAR)
1378	#else
1379	CALL WSCAL(LU_NONZERO,(-ONE),Jac0,1)
1380	#endif
1381	DO iadj = 1, NADJ
1382	#ifdef FULL_ALGEBRA
1383	Fcn0(1:NVAR,iadj) = MATMUL(TRANSPOSE(Jac0),Y(1:NVAR,iadj))
1384	#else
1385	CALL JacTR_SP_Vec(Jac0,Y(1,iadj),Fcn0(1,iadj))
1386	#endif
1387	END DO
1388
1389	!~~~> Repeat step calculation until current step accepted
1390	UntilAccepted: DO
1391
1392	CALL ros_PrepareMatrix(H,Direction,ros_Gamma(1), &
1393	Jac0,Ghimj,Pivot,Singular)
1394	IF (Singular) THEN ! More than 5 consecutive failed decompositions
1395	CALL ros_ErrorMsg(-8,T,H,IERR)
1396	RETURN
1397	END IF
1398
1399	!~~~> Compute the stages
1400	Stage: DO istage = 1, ros_S
1401
1402	! Current istage offset. Current istage vector is K(ioffset+1:ioffset+NVAR)
1403	ioffset = NVAR*(istage-1)
1404
1405	! For the 1st istage the function has been computed previously
1406	IF ( istage == 1 ) THEN
1407	DO iadj = 1, NADJ
1408	CALL WCOPY(NVAR,Fcn0(1,iadj),1,Fcn(1,iadj),1)
1409	END DO
1410	! istage>1 and a new function evaluation is needed at the current istage
1411	ELSEIF ( ros_NewF(istage) ) THEN
1412	CALL WCOPY(NVAR*NADJ,Y,1,Ynew,1)
1413	DO j = 1, istage-1
1414	DO iadj = 1, NADJ
1415	CALL WAXPY(NVAR,ros_A((istage-1)*(istage-2)/2+j), &
1416	K(NVAR*(j-1)+1,iadj),1,Ynew(1,iadj),1)
1417	END DO
1418	END DO
1419	Tau = T + ros_Alpha(istage)DirectionH
1420	! CALL FunTemplate(Tau,Ynew,Fcn)
1421	! ISTATUS(Nfun) = ISTATUS(Nfun) + 1
1422	CALL ros_cadj_Y( Tau, Y0 )
1423	CALL JacTemplate(Tau, Y0, Jac)
1424	ISTATUS(Njac) = ISTATUS(Njac) + 1
1425	#ifdef FULL_ALGEBRA
1426	Jac(1:NVAR,1:NVAR) = -Jac(1:NVAR,1:NVAR)
1427	#else
1428	CALL WSCAL(LU_NONZERO,(-ONE),Jac,1)
1429	#endif
1430	DO iadj = 1, NADJ
1431	#ifdef FULL_ALGEBRA
1432	Fcn(1:NVAR,iadj) = MATMUL(TRANSPOSE(Jac),Ynew(1:NVAR,iadj))
1433	#else
1434	CALL JacTR_SP_Vec(Jac,Ynew(1,iadj),Fcn(1,iadj))
1435	#endif
1436	!CALL WSCAL(NVAR,(-ONE),Fcn(1,iadj),1)
1437	END DO
1438	END IF ! if istage == 1 elseif ros_NewF(istage)
1439
1440	DO iadj = 1, NADJ
1441	CALL WCOPY(NVAR,Fcn(1,iadj),1,K(ioffset+1,iadj),1)
1442	END DO
1443	DO j = 1, istage-1
1444	HC = ros_C((istage-1)(istage-2)/2+j)/(DirectionH)
1445	DO iadj = 1, NADJ
1446	CALL WAXPY(NVAR,HC,K(NVAR*(j-1)+1,iadj),1, &
1447	K(ioffset+1,iadj),1)
1448	END DO
1449	END DO
1450	IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN
1451	HG = DirectionHros_Gamma(istage)
1452	DO iadj = 1, NADJ
1453	CALL WAXPY(NVAR,HG,dFdT(1,iadj),1,K(ioffset+1,iadj),1)
1454	END DO
1455	END IF
1456	DO iadj = 1, NADJ
1457	CALL ros_Solve('T', Ghimj, Pivot, K(ioffset+1,iadj))
1458	END DO
1459
1460	END DO Stage
1461
1462
1463	!~~~> Compute the new solution
1464	DO iadj = 1, NADJ
1465	CALL WCOPY(NVAR,Y(1,iadj),1,Ynew(1,iadj),1)
1466	DO j=1,ros_S
1467	CALL WAXPY(NVAR,ros_M(j),K(NVAR*(j-1)+1,iadj),1,Ynew(1,iadj),1)
1468	END DO
1469	END DO
1470
1471	!~~~> Compute the error estimation
1472	CALL WSCAL(NVAR*NADJ,ZERO,Yerr,1)
1473	DO j=1,ros_S
1474	DO iadj = 1, NADJ
1475	CALL WAXPY(NVAR,ros_E(j),K(NVAR*(j-1)+1,iadj),1,Yerr(1,iadj),1)
1476	END DO
1477	END DO
1478	!~~~> Max error among all adjoint components
1479	iadj = 1
1480	Err = ros_ErrorNorm ( Y(1,iadj), Ynew(1,iadj), Yerr(1,iadj), &
1481	AbsTol_adj(1,iadj), RelTol_adj(1,iadj), VectorTol )
1482
1483	!~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax
1484	Fac = MIN(FacMax,MAX(FacMin,FacSafe/Err**(ONE/ros_ELO)))
1485	Hnew = H*Fac
1486
1487	!~~~> Check the error magnitude and adjust step size
1488	! ISTATUS(Nstp) = ISTATUS(Nstp) + 1
1489	IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN !~~~> Accept step
1490	ISTATUS(Nacc) = ISTATUS(Nacc) + 1
1491	CALL WCOPY(NVAR*NADJ,Ynew,1,Y,1)
1492	T = T + Direction*H
1493	Hnew = MAX(Hmin,MIN(Hnew,Hmax))
1494	IF (RejectLastH) THEN ! No step size increase after a rejected step
1495	Hnew = MIN(Hnew,H)
1496	END IF
1497	RSTATUS(Nhexit) = H
1498	RSTATUS(Nhnew) = Hnew
1499	RSTATUS(Ntexit) = T
1500	RejectLastH = .FALSE.
1501	RejectMoreH = .FALSE.
1502	H = Hnew
1503	EXIT UntilAccepted ! EXIT THE LOOP: WHILE STEP NOT ACCEPTED
1504	ELSE !~~~> Reject step
1505	IF (RejectMoreH) THEN
1506	Hnew = H*FacRej
1507	END IF
1508	RejectMoreH = RejectLastH
1509	RejectLastH = .TRUE.
1510	H = Hnew
1511	IF (ISTATUS(Nacc) >= 1) THEN
1512	ISTATUS(Nrej) = ISTATUS(Nrej) + 1
1513	END IF
1514	END IF ! Err <= 1
1515
1516	END DO UntilAccepted
1517
1518	END DO TimeLoop
1519
1520	!~~~> Succesful exit
1521	IERR = 1 !~~~> The integration was successful
1522
1523	END SUBROUTINE ros_CadjInt
1524
1525
1526	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1527	SUBROUTINE ros_SimpleCadjInt ( &
1528	NADJ, Y, &
1529	Tstart, Tend, T, &
1530	!~~~> Error indicator
1531	IERR )
1532	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1533	! Template for the implementation of a generic RosenbrockADJ method
1534	! defined by ros_S (no of stages)
1535	! and its coefficients ros_{A,C,M,E,Alpha,Gamma}
1536	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1537
1538	IMPLICIT NONE
1539
1540	!~~~> Input: the initial condition at Tstart; Output: the solution at T
1541	INTEGER, INTENT(IN) :: NADJ
1542	KPP_REAL, INTENT(INOUT) :: Y(NVAR,NADJ)
1543	!~~~> Input: integration interval
1544	KPP_REAL, INTENT(IN) :: Tstart,Tend
1545	!~~~> Output: time at which the solution is returned (T=Tend if success)
1546	KPP_REAL, INTENT(OUT) :: T
1547	!~~~> Output: Error indicator
1548	INTEGER, INTENT(OUT) :: IERR
1549	! ~~~~ Local variables
1550	KPP_REAL :: Y0(NVAR)
1551	KPP_REAL :: Ynew(NVAR,NADJ), Fcn0(NVAR,NADJ), Fcn(NVAR,NADJ)
1552	KPP_REAL :: K(NVAR*ros_S,NADJ), dFdT(NVAR,NADJ)
1553	#ifdef FULL_ALGEBRA
1554	KPP_REAL,DIMENSION(NVAR,NVAR) :: Jac0, Ghimj, Jac, dJdT
1555	#else
1556	KPP_REAL,DIMENSION(LU_NONZERO) :: Jac0, Ghimj, Jac, dJdT
1557	#endif
1558	KPP_REAL :: H, HC, HG, Tau
1559	KPP_REAL :: ghinv
1560	INTEGER :: Pivot(NVAR), Direction, ioffset, i, j, istage, iadj
1561	INTEGER :: istack
1562	!~~~> Local parameters
1563	KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
1564	KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5
1565	!~~~> Locally called functions
1566	! KPP_REAL WLAMCH
1567	! EXTERNAL WLAMCH
1568	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1569
1570
1571	!~~~> INITIAL PREPARATIONS
1572
1573	IF (Tend >= Tstart) THEN
1574	Direction = -1
1575	ELSE
1576	Direction = +1
1577	END IF
1578
1579	!~~~> Time loop begins below
1580	TimeLoop: DO istack = stack_ptr,2,-1
1581
1582	T = chk_T(istack)
1583	H = chk_H(istack-1)
1584	!CALL WCOPY(NVAR,chk_Y(1,istack),1,Y0,1)
1585	Y0(1:NVAR) = chk_Y(1:NVAR,istack)
1586
1587	!~~~> Compute the Jacobian at current time
1588	CALL JacTemplate(T, Y0, Jac0)
1589	ISTATUS(Njac) = ISTATUS(Njac) + 1
1590
1591	!~~~> Compute the function derivative with respect to T
1592	IF (.NOT.Autonomous) THEN
1593	CALL ros_JacTimeDerivative ( T, Roundoff, Y0, &
1594	Jac0, dJdT )
1595	DO iadj = 1, NADJ
1596	#ifdef FULL_ALGEBRA
1597	dFdT(1:NVAR,iadj) = MATMUL(TRANSPOSE(dJdT),Y(1:NVAR,iadj))
1598	#else
1599	CALL JacTR_SP_Vec(dJdT,Y(1,iadj),dFdT(1,iadj))
1600	#endif
1601	CALL WSCAL(NVAR,(-ONE),dFdT(1,iadj),1)
1602	END DO
1603	END IF
1604
1605	!~~~> Ydot = -J^T*Y
1606	#ifdef FULL_ALGEBRA
1607	Jac0(1:NVAR,1:NVAR) = -Jac0(1:NVAR,1:NVAR)
1608	#else
1609	CALL WSCAL(LU_NONZERO,(-ONE),Jac0,1)
1610	#endif
1611	DO iadj = 1, NADJ
1612	#ifdef FULL_ALGEBRA
1613	Fcn0(1:NVAR,iadj) = MATMUL(TRANSPOSE(Jac0),Y(1:NVAR,iadj))
1614	#else
1615	CALL JacTR_SP_Vec(Jac0,Y(1,iadj),Fcn0(1,iadj))
1616	#endif
1617	END DO
1618
1619	!~~~> Construct Ghimj = 1/(H*ham) - Jac0
1620	ghinv = ONE/(DirectionHros_Gamma(1))
1621	#ifdef FULL_ALGEBRA
1622	Ghimj(1:NVAR,1:NVAR) = -Jac0(1:NVAR,1:NVAR)
1623	DO i=1,NVAR
1624	Ghimj(i,i) = Ghimj(i,i)+ghinv
1625	END DO
1626	#else
1627	CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1)
1628	CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1)
1629	DO i=1,NVAR
1630	Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+ghinv
1631	END DO
1632	#endif
1633	!~~~> Compute LU decomposition
1634	CALL ros_Decomp( Ghimj, Pivot, j )
1635	IF (j /= 0) THEN
1636	CALL ros_ErrorMsg(-8,T,H,IERR)
1637	PRINT*,' The matrix is singular !'
1638	STOP
1639	END IF
1640
1641	!~~~> Compute the stages
1642	Stage: DO istage = 1, ros_S
1643
1644	! Current istage offset. Current istage vector is K(ioffset+1:ioffset+NVAR)
1645	ioffset = NVAR*(istage-1)
1646
1647	! For the 1st istage the function has been computed previously
1648	IF ( istage == 1 ) THEN
1649	DO iadj = 1, NADJ
1650	CALL WCOPY(NVAR,Fcn0(1,iadj),1,Fcn(1,iadj),1)
1651	END DO
1652	! istage>1 and a new function evaluation is needed at the current istage
1653	ELSEIF ( ros_NewF(istage) ) THEN
1654	CALL WCOPY(NVAR*NADJ,Y,1,Ynew,1)
1655	DO j = 1, istage-1
1656	DO iadj = 1, NADJ
1657	CALL WAXPY(NVAR,ros_A((istage-1)*(istage-2)/2+j), &
1658	K(NVAR*(j-1)+1,iadj),1,Ynew(1,iadj),1)
1659	END DO
1660	END DO
1661	Tau = T + ros_Alpha(istage)DirectionH
1662	CALL ros_Hermite3( chk_T(istack-1), chk_T(istack), Tau, &
1663	chk_Y(1:NVAR,istack-1), chk_Y(1:NVAR,istack), &
1664	chk_dY(1:NVAR,istack-1), chk_dY(1:NVAR,istack), Y0 )
1665	CALL JacTemplate(Tau, Y0, Jac)
1666	ISTATUS(Njac) = ISTATUS(Njac) + 1
1667	#ifdef FULL_ALGEBRA
1668	Jac(1:NVAR,1:NVAR) = -Jac(1:NVAR,1:NVAR)
1669	#else
1670	CALL WSCAL(LU_NONZERO,(-ONE),Jac,1)
1671	#endif
1672	DO iadj = 1, NADJ
1673	#ifdef FULL_ALGEBRA
1674	Fcn(1:NVAR,iadj) = MATMUL(TRANSPOSE(Jac),Ynew(1:NVAR,iadj))
1675	#else
1676	CALL JacTR_SP_Vec(Jac,Ynew(1,iadj),Fcn(1,iadj))
1677	#endif
1678	END DO
1679	END IF ! if istage == 1 elseif ros_NewF(istage)
1680
1681	DO iadj = 1, NADJ
1682	CALL WCOPY(NVAR,Fcn(1,iadj),1,K(ioffset+1,iadj),1)
1683	END DO
1684	DO j = 1, istage-1
1685	HC = ros_C((istage-1)(istage-2)/2+j)/(DirectionH)
1686	DO iadj = 1, NADJ
1687	CALL WAXPY(NVAR,HC,K(NVAR*(j-1)+1,iadj),1, &
1688	K(ioffset+1,iadj),1)
1689	END DO
1690	END DO
1691	IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN
1692	HG = DirectionHros_Gamma(istage)
1693	DO iadj = 1, NADJ
1694	CALL WAXPY(NVAR,HG,dFdT(1,iadj),1,K(ioffset+1,iadj),1)
1695	END DO
1696	END IF
1697	DO iadj = 1, NADJ
1698	CALL ros_Solve('T', Ghimj, Pivot, K(ioffset+1,iadj))
1699	END DO
1700
1701	END DO Stage
1702
1703
1704	!~~~> Compute the new solution
1705	DO iadj = 1, NADJ
1706	DO j=1,ros_S
1707	CALL WAXPY(NVAR,ros_M(j),K(NVAR*(j-1)+1,iadj),1,Y(1,iadj),1)
1708	END DO
1709	END DO
1710
1711	END DO TimeLoop
1712
1713	!~~~> Succesful exit
1714	IERR = 1 !~~~> The integration was successful
1715
1716	END SUBROUTINE ros_SimpleCadjInt
1717
1718	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1719	KPP_REAL FUNCTION ros_ErrorNorm ( Y, Ynew, Yerr, &
1720	AbsTol, RelTol, VectorTol )
1721	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1722	!~~~> Computes the "scaled norm" of the error vector Yerr
1723	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1724	IMPLICIT NONE
1725
1726	! Input arguments
1727	KPP_REAL, INTENT(IN) :: Y(NVAR), Ynew(NVAR), &
1728	Yerr(NVAR), AbsTol(NVAR), RelTol(NVAR)
1729	LOGICAL, INTENT(IN) :: VectorTol
1730	! Local variables
1731	KPP_REAL :: Err, Scale, Ymax
1732	INTEGER :: i
1733
1734	Err = ZERO
1735	DO i=1,NVAR
1736	Ymax = MAX(ABS(Y(i)),ABS(Ynew(i)))
1737	IF (VectorTol) THEN
1738	Scale = AbsTol(i)+RelTol(i)*Ymax
1739	ELSE
1740	Scale = AbsTol(1)+RelTol(1)*Ymax
1741	END IF
1742	Err = Err+(Yerr(i)/Scale)**2
1743	END DO
1744	Err = SQRT(Err/NVAR)
1745
1746	ros_ErrorNorm = MAX(Err,1.0d-10)
1747
1748	END FUNCTION ros_ErrorNorm
1749
1750
1751	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1752	SUBROUTINE ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, dFdT )
1753	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1754	!~~~> The time partial derivative of the function by finite differences
1755	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1756	IMPLICIT NONE
1757
1758	!~~~> Input arguments
1759	KPP_REAL, INTENT(IN) :: T, Roundoff, Y(NVAR), Fcn0(NVAR)
1760	!~~~> Output arguments
1761	KPP_REAL, INTENT(OUT) :: dFdT(NVAR)
1762	!~~~> Local variables
1763	KPP_REAL :: Delta
1764	KPP_REAL, PARAMETER :: ONE = 1.0d0, DeltaMin = 1.0d-6
1765
1766	Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T))
1767	CALL FunTemplate(T+Delta,Y,dFdT)
1768	ISTATUS(Nfun) = ISTATUS(Nfun) + 1
1769	CALL WAXPY(NVAR,(-ONE),Fcn0,1,dFdT,1)
1770	CALL WSCAL(NVAR,(ONE/Delta),dFdT,1)
1771
1772	END SUBROUTINE ros_FunTimeDerivative
1773
1774
1775	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1776	SUBROUTINE ros_JacTimeDerivative ( T, Roundoff, Y, &
1777	Jac0, dJdT )
1778	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1779	!~~~> The time partial derivative of the Jacobian by finite differences
1780	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1781	IMPLICIT NONE
1782
1783	!~~~> Arguments
1784	KPP_REAL, INTENT(IN) :: T, Roundoff, Y(NVAR)
1785	#ifdef FULL_ALGEBRA
1786	KPP_REAL, INTENT(IN) :: Jac0(NVAR,NVAR)
1787	KPP_REAL, INTENT(OUT) :: dJdT(NVAR,NVAR)
1788	#else
1789	KPP_REAL, INTENT(IN) :: Jac0(LU_NONZERO)
1790	KPP_REAL, INTENT(OUT) :: dJdT(LU_NONZERO)
1791	#endif
1792	!~~~> Local variables
1793	KPP_REAL :: Delta
1794
1795	Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T))
1796	CALL JacTemplate(T+Delta,Y,dJdT)
1797	ISTATUS(Njac) = ISTATUS(Njac) + 1
1798	#ifdef FULL_ALGEBRA
1799	CALL WAXPY(NVAR*NVAR,(-ONE),Jac0,1,dJdT,1)
1800	CALL WSCAL(NVAR*NVAR,(ONE/Delta),dJdT,1)
1801	#else
1802	CALL WAXPY(LU_NONZERO,(-ONE),Jac0,1,dJdT,1)
1803	CALL WSCAL(LU_NONZERO,(ONE/Delta),dJdT,1)
1804	#endif
1805
1806	END SUBROUTINE ros_JacTimeDerivative
1807
1808
1809
1810	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1811	SUBROUTINE ros_PrepareMatrix ( H, Direction, gam, &
1812	Jac0, Ghimj, Pivot, Singular )
1813	! --- --- --- --- --- --- --- --- --- --- --- --- ---
1814	! Prepares the LHS matrix for stage calculations
1815	! 1. Construct Ghimj = 1/(H*gam) - Jac0
1816	! "(Gamma H) Inverse Minus Jacobian"
1817	! 2. Repeat LU decomposition of Ghimj until successful.
1818	! -half the step size if LU decomposition fails and retry
1819	! -exit after 5 consecutive fails
1820	! --- --- --- --- --- --- --- --- --- --- --- --- ---
1821	IMPLICIT NONE
1822
1823	!~~~> Input arguments
1824	#ifdef FULL_ALGEBRA
1825	KPP_REAL, INTENT(IN) :: Jac0(NVAR,NVAR)
1826	#else
1827	KPP_REAL, INTENT(IN) :: Jac0(LU_NONZERO)
1828	#endif
1829	KPP_REAL, INTENT(IN) :: gam
1830	INTEGER, INTENT(IN) :: Direction
1831	!~~~> Output arguments
1832	#ifdef FULL_ALGEBRA
1833	KPP_REAL, INTENT(OUT) :: Ghimj(NVAR,NVAR)
1834	#else
1835	KPP_REAL, INTENT(OUT) :: Ghimj(LU_NONZERO)
1836	#endif
1837	LOGICAL, INTENT(OUT) :: Singular
1838	INTEGER, INTENT(OUT) :: Pivot(NVAR)
1839	!~~~> Inout arguments
1840	KPP_REAL, INTENT(INOUT) :: H ! step size is decreased when LU fails
1841	!~~~> Local variables
1842	INTEGER :: i, ISING, Nconsecutive
1843	KPP_REAL :: ghinv
1844	KPP_REAL, PARAMETER :: ONE = 1.0_dp, HALF = 0.5_dp
1845
1846	Nconsecutive = 0
1847	Singular = .TRUE.
1848
1849	DO WHILE (Singular)
1850
1851	!~~~> Construct Ghimj = 1/(H*gam) - Jac0
1852	#ifdef FULL_ALGEBRA
1853	CALL WCOPY(NVAR*NVAR,Jac0,1,Ghimj,1)
1854	CALL WSCAL(NVAR*NVAR,(-ONE),Ghimj,1)
1855	ghinv = ONE/(DirectionHgam)
1856	DO i=1,NVAR
1857	Ghimj(i,i) = Ghimj(i,i)+ghinv
1858	END DO
1859	#else
1860	CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1)
1861	CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1)
1862	ghinv = ONE/(DirectionHgam)
1863	DO i=1,NVAR
1864	Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+ghinv
1865	END DO
1866	#endif
1867	!~~~> Compute LU decomposition
1868	CALL ros_Decomp( Ghimj, Pivot, ISING )
1869	IF (ISING == 0) THEN
1870	!~~~> If successful done
1871	Singular = .FALSE.
1872	ELSE ! ISING .ne. 0
1873	!~~~> If unsuccessful half the step size; if 5 consecutive fails then return
1874	ISTATUS(Nsng) = ISTATUS(Nsng) + 1
1875	Nconsecutive = Nconsecutive+1
1876	Singular = .TRUE.
1877	PRINT*,'Warning: LU Decomposition returned ISING = ',ISING
1878	IF (Nconsecutive <= 5) THEN ! Less than 5 consecutive failed decompositions
1879	H = H*HALF
1880	ELSE ! More than 5 consecutive failed decompositions
1881	RETURN
1882	END IF ! Nconsecutive
1883	END IF ! ISING
1884
1885	END DO ! WHILE Singular
1886
1887	END SUBROUTINE ros_PrepareMatrix
1888
1889
1890	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1891	SUBROUTINE ros_Decomp( A, Pivot, ISING )
1892	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1893	! Template for the LU decomposition
1894	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1895	IMPLICIT NONE
1896	!~~~> Inout variables
1897	#ifdef FULL_ALGEBRA
1898	KPP_REAL, INTENT(INOUT) :: A(NVAR,NVAR)
1899	#else
1900	KPP_REAL, INTENT(INOUT) :: A(LU_NONZERO)
1901	#endif
1902	!~~~> Output variables
1903	INTEGER, INTENT(OUT) :: Pivot(NVAR), ISING
1904
1905	#ifdef FULL_ALGEBRA
1906	CALL DGETRF( NVAR, NVAR, A, NVAR, Pivot, ISING )
1907	#else
1908	CALL KppDecomp ( A, ISING )
1909	Pivot(1) = 1
1910	#endif
1911	ISTATUS(Ndec) = ISTATUS(Ndec) + 1
1912
1913	END SUBROUTINE ros_Decomp
1914
1915
1916	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1917	SUBROUTINE ros_Solve( How, A, Pivot, b )
1918	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1919	! Template for the forward/backward substitution (using pre-computed LU decomposition)
1920	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1921	IMPLICIT NONE
1922	!~~~> Input variables
1923	CHARACTER, INTENT(IN) :: How
1924	#ifdef FULL_ALGEBRA
1925	KPP_REAL, INTENT(IN) :: A(NVAR,NVAR)
1926	INTEGER :: ISING
1927	#else
1928	KPP_REAL, INTENT(IN) :: A(LU_NONZERO)
1929	#endif
1930	INTEGER, INTENT(IN) :: Pivot(NVAR)
1931	!~~~> InOut variables
1932	KPP_REAL, INTENT(INOUT) :: b(NVAR)
1933
1934	SELECT CASE (How)
1935	CASE ('N')
1936	#ifdef FULL_ALGEBRA
1937	CALL DGETRS( 'N', NVAR , 1, A, NVAR, Pivot, b, NVAR, ISING )
1938	#else
1939	CALL KppSolve( A, b )
1940	#endif
1941	CASE ('T')
1942	#ifdef FULL_ALGEBRA
1943	CALL DGETRS( 'T', NVAR , 1, A, NVAR, Pivot, b, NVAR, ISING )
1944	#else
1945	CALL KppSolveTR( A, b, b )
1946	#endif
1947	CASE DEFAULT
1948	PRINT*,'Error: unknown argument in ros_Solve: How=',How
1949	STOP
1950	END SELECT
1951	ISTATUS(Nsol) = ISTATUS(Nsol) + 1
1952
1953	END SUBROUTINE ros_Solve
1954
1955
1956
1957	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1958	SUBROUTINE ros_cadj_Y( T, Y )
1959	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1960	! Finds the solution Y at T by interpolating the stored forward trajectory
1961	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1962	IMPLICIT NONE
1963	!~~~> Input variables
1964	KPP_REAL, INTENT(IN) :: T
1965	!~~~> Output variables
1966	KPP_REAL, INTENT(OUT) :: Y(NVAR)
1967	!~~~> Local variables
1968	INTEGER :: i
1969	KPP_REAL, PARAMETER :: ONE = 1.0d0
1970
1971	! chk_H, chk_T, chk_Y, chk_dY, chk_d2Y
1972
1973	IF( (T < chk_T(1)).OR.(T> chk_T(stack_ptr)) ) THEN
1974	PRINT*,'Cannot locate solution at T = ',T
1975	PRINT*,'Stored trajectory is between Tstart = ',chk_T(1)
1976	PRINT*,' and Tend = ',chk_T(stack_ptr)
1977	STOP
1978	END IF
1979	DO i = 1, stack_ptr-1
1980	IF( (T>= chk_T(i)).AND.(T<= chk_T(i+1)) ) EXIT
1981	END DO
1982
1983
1984	! IF (.FALSE.) THEN
1985	!
1986	! CALL ros_Hermite5( chk_T(i), chk_T(i+1), T, &
1987	! chk_Y(1,i), chk_Y(1,i+1), &
1988	! chk_dY(1,i), chk_dY(1,i+1), &
1989	! chk_d2Y(1,i), chk_d2Y(1,i+1), Y )
1990	!
1991	! ELSE
1992
1993	CALL ros_Hermite3( chk_T(i), chk_T(i+1), T, &
1994	chk_Y(1:NVAR,i), chk_Y(1:NVAR,i+1), &
1995	chk_dY(1:NVAR,i), chk_dY(1:NVAR,i+1), &
1996	Y )
1997
1998	!
1999	! END IF
2000
2001	END SUBROUTINE ros_cadj_Y
2002
2003
2004	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2005	SUBROUTINE ros_Hermite3( a, b, T, Ya, Yb, Ja, Jb, Y )
2006	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2007	! Template for Hermite interpolation of order 5 on the interval [a,b]
2008	! P = c(1) + c(2)(x-a) + ... + c(4)(x-a)^3
2009	! P[a,b] = [Ya,Yb], P'[a,b] = [Ja,Jb]
2010	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2011	IMPLICIT NONE
2012	!~~~> Input variables
2013	KPP_REAL, INTENT(IN) :: a, b, T, Ya(NVAR), Yb(NVAR)
2014	KPP_REAL, INTENT(IN) :: Ja(NVAR), Jb(NVAR)
2015	!~~~> Output variables
2016	KPP_REAL, INTENT(OUT) :: Y(NVAR)
2017	!~~~> Local variables
2018	KPP_REAL :: Tau, amb(3), C(NVAR,4)
2019	KPP_REAL, PARAMETER :: ZERO = 0.0d0
2020	INTEGER :: i, j
2021
2022	amb(1) = 1.0d0/(a-b)
2023	DO i=2,3
2024	amb(i) = amb(i-1)*amb(1)
2025	END DO
2026
2027
2028	! c(1) = ya;
2029	CALL WCOPY(NVAR,Ya,1,C(1,1),1)
2030	! c(2) = ja;
2031	CALL WCOPY(NVAR,Ja,1,C(1,2),1)
2032	! c(3) = 2/(a-b)ja + 1/(a-b)jb - 3/(a - b)^2ya + 3/(a - b)^2yb ;
2033	CALL WCOPY(NVAR,Ya,1,C(1,3),1)
2034	CALL WSCAL(NVAR,-3.0*amb(2),C(1,3),1)
2035	CALL WAXPY(NVAR,3.0*amb(2),Yb,1,C(1,3),1)
2036	CALL WAXPY(NVAR,2.0*amb(1),Ja,1,C(1,3),1)
2037	CALL WAXPY(NVAR,amb(1),Jb,1,C(1,3),1)
2038	! c(4) = 1/(a-b)^2ja + 1/(a-b)^2jb - 2/(a-b)^3ya + 2/(a-b)^3yb ;
2039	CALL WCOPY(NVAR,Ya,1,C(1,4),1)
2040	CALL WSCAL(NVAR,-2.0*amb(3),C(1,4),1)
2041	CALL WAXPY(NVAR,2.0*amb(3),Yb,1,C(1,4),1)
2042	CALL WAXPY(NVAR,amb(2),Ja,1,C(1,4),1)
2043	CALL WAXPY(NVAR,amb(2),Jb,1,C(1,4),1)
2044
2045	Tau = T - a
2046	CALL WCOPY(NVAR,C(1,4),1,Y,1)
2047	CALL WSCAL(NVAR,Tau**3,Y,1)
2048	DO j = 3,1,-1
2049	CALL WAXPY(NVAR,TAU**(j-1),C(1,j),1,Y,1)
2050	END DO
2051
2052	END SUBROUTINE ros_Hermite3
2053
2054	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2055	SUBROUTINE ros_Hermite5( a, b, T, Ya, Yb, Ja, Jb, Ha, Hb, Y )
2056	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2057	! Template for Hermite interpolation of order 5 on the interval [a,b]
2058	! P = c(1) + c(2)(x-a) + ... + c(6)(x-a)^5
2059	! P[a,b] = [Ya,Yb], P'[a,b] = [Ja,Jb], P"[a,b] = [Ha,Hb]
2060	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2061	IMPLICIT NONE
2062	!~~~> Input variables
2063	KPP_REAL, INTENT(IN) :: a, b, T, Ya(NVAR), Yb(NVAR)
2064	KPP_REAL, INTENT(IN) :: Ja(NVAR), Jb(NVAR), Ha(NVAR), Hb(NVAR)
2065	!~~~> Output variables
2066	KPP_REAL, INTENT(OUT) :: Y(NVAR)
2067	!~~~> Local variables
2068	KPP_REAL :: Tau, amb(5), C(NVAR,6)
2069	KPP_REAL, PARAMETER :: ZERO = 0.0d0, HALF = 0.5d0
2070	INTEGER :: i, j
2071
2072	amb(1) = 1.0d0/(a-b)
2073	DO i=2,5
2074	amb(i) = amb(i-1)*amb(1)
2075	END DO
2076
2077	! c(1) = ya;
2078	CALL WCOPY(NVAR,Ya,1,C(1,1),1)
2079	! c(2) = ja;
2080	CALL WCOPY(NVAR,Ja,1,C(1,2),1)
2081	! c(3) = ha/2;
2082	CALL WCOPY(NVAR,Ha,1,C(1,3),1)
2083	CALL WSCAL(NVAR,HALF,C(1,3),1)
2084
2085	! c(4) = 10amb(3)ya - 10amb(3)yb - 6amb(2)ja - 4amb(2)jb + 1.5amb(1)ha - 0.5amb(1)hb ;
2086	CALL WCOPY(NVAR,Ya,1,C(1,4),1)
2087	CALL WSCAL(NVAR,10.0*amb(3),C(1,4),1)
2088	CALL WAXPY(NVAR,-10.0*amb(3),Yb,1,C(1,4),1)
2089	CALL WAXPY(NVAR,-6.0*amb(2),Ja,1,C(1,4),1)
2090	CALL WAXPY(NVAR,-4.0*amb(2),Jb,1,C(1,4),1)
2091	CALL WAXPY(NVAR, 1.5*amb(1),Ha,1,C(1,4),1)
2092	CALL WAXPY(NVAR,-0.5*amb(1),Hb,1,C(1,4),1)
2093
2094	! c(5) = 15amb(4)ya - 15amb(4)yb - 8.amb(3)ja - 7amb(3)jb + 1.5amb(2)ha - 1amb(2)hb ;
2095	CALL WCOPY(NVAR,Ya,1,C(1,5),1)
2096	CALL WSCAL(NVAR, 15.0*amb(4),C(1,5),1)
2097	CALL WAXPY(NVAR,-15.0*amb(4),Yb,1,C(1,5),1)
2098	CALL WAXPY(NVAR,-8.0*amb(3),Ja,1,C(1,5),1)
2099	CALL WAXPY(NVAR,-7.0*amb(3),Jb,1,C(1,5),1)
2100	CALL WAXPY(NVAR,1.5*amb(2),Ha,1,C(1,5),1)
2101	CALL WAXPY(NVAR,-amb(2),Hb,1,C(1,5),1)
2102
2103	! c(6) = 6amb(5)ya - 6amb(5)yb - 3.amb(4)ja - 3.amb(4)jb + 0.5amb(3)ha -0.5amb(3)hb ;
2104	CALL WCOPY(NVAR,Ya,1,C(1,6),1)
2105	CALL WSCAL(NVAR, 6.0*amb(5),C(1,6),1)
2106	CALL WAXPY(NVAR,-6.0*amb(5),Yb,1,C(1,6),1)
2107	CALL WAXPY(NVAR,-3.0*amb(4),Ja,1,C(1,6),1)
2108	CALL WAXPY(NVAR,-3.0*amb(4),Jb,1,C(1,6),1)
2109	CALL WAXPY(NVAR, 0.5*amb(3),Ha,1,C(1,6),1)
2110	CALL WAXPY(NVAR,-0.5*amb(3),Hb,1,C(1,6),1)
2111
2112	Tau = T - a
2113	CALL WCOPY(NVAR,C(1,6),1,Y,1)
2114	DO j = 5,1,-1
2115	CALL WSCAL(NVAR,Tau,Y,1)
2116	CALL WAXPY(NVAR,ONE,C(1,j),1,Y,1)
2117	END DO
2118
2119	END SUBROUTINE ros_Hermite5
2120
2121	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2122	SUBROUTINE Ros2
2123	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2124	! --- AN L-STABLE METHOD, 2 stages, order 2
2125	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2126
2127	IMPLICIT NONE
2128	DOUBLE PRECISION g
2129
2130	g = 1.0d0 + 1.0d0/SQRT(2.0d0)
2131
2132	rosMethod = RS2
2133	!~~~> Name of the method
2134	ros_Name = 'ROS-2'
2135	!~~~> Number of stages
2136	ros_S = 2
2137
2138	!~~~> The coefficient matrices A and C are strictly lower triangular.
2139	! The lower triangular (subdiagonal) elements are stored in row-wise order:
2140	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
2141	! The general mapping formula is:
2142	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
2143	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
2144
2145	ros_A(1) = (1.d0)/g
2146	ros_C(1) = (-2.d0)/g
2147	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
2148	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
2149	ros_NewF(1) = .TRUE.
2150	ros_NewF(2) = .TRUE.
2151	!~~~> M_i = Coefficients for new step solution
2152	ros_M(1)= (3.d0)/(2.d0*g)
2153	ros_M(2)= (1.d0)/(2.d0*g)
2154	! E_i = Coefficients for error estimator
2155	ros_E(1) = 1.d0/(2.d0*g)
2156	ros_E(2) = 1.d0/(2.d0*g)
2157	!~~~> ros_ELO = estimator of local order - the minimum between the
2158	! main and the embedded scheme orders plus one
2159	ros_ELO = 2.0d0
2160	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
2161	ros_Alpha(1) = 0.0d0
2162	ros_Alpha(2) = 1.0d0
2163	!~~~> Gamma_i = \sum_j gamma_{i,j}
2164	ros_Gamma(1) = g
2165	ros_Gamma(2) =-g
2166
2167	END SUBROUTINE Ros2
2168
2169
2170	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2171	SUBROUTINE Ros3
2172	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2173	! --- AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations
2174	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2175
2176	IMPLICIT NONE
2177
2178	rosMethod = RS3
2179	!~~~> Name of the method
2180	ros_Name = 'ROS-3'
2181	!~~~> Number of stages
2182	ros_S = 3
2183
2184	!~~~> The coefficient matrices A and C are strictly lower triangular.
2185	! The lower triangular (subdiagonal) elements are stored in row-wise order:
2186	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
2187	! The general mapping formula is:
2188	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
2189	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
2190
2191	ros_A(1)= 1.d0
2192	ros_A(2)= 1.d0
2193	ros_A(3)= 0.d0
2194
2195	ros_C(1) = -0.10156171083877702091975600115545d+01
2196	ros_C(2) = 0.40759956452537699824805835358067d+01
2197	ros_C(3) = 0.92076794298330791242156818474003d+01
2198	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
2199	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
2200	ros_NewF(1) = .TRUE.
2201	ros_NewF(2) = .TRUE.
2202	ros_NewF(3) = .FALSE.
2203	!~~~> M_i = Coefficients for new step solution
2204	ros_M(1) = 0.1d+01
2205	ros_M(2) = 0.61697947043828245592553615689730d+01
2206	ros_M(3) = -0.42772256543218573326238373806514d+00
2207	! E_i = Coefficients for error estimator
2208	ros_E(1) = 0.5d+00
2209	ros_E(2) = -0.29079558716805469821718236208017d+01
2210	ros_E(3) = 0.22354069897811569627360909276199d+00
2211	!~~~> ros_ELO = estimator of local order - the minimum between the
2212	! main and the embedded scheme orders plus 1
2213	ros_ELO = 3.0d0
2214	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
2215	ros_Alpha(1)= 0.0d+00
2216	ros_Alpha(2)= 0.43586652150845899941601945119356d+00
2217	ros_Alpha(3)= 0.43586652150845899941601945119356d+00
2218	!~~~> Gamma_i = \sum_j gamma_{i,j}
2219	ros_Gamma(1)= 0.43586652150845899941601945119356d+00
2220	ros_Gamma(2)= 0.24291996454816804366592249683314d+00
2221	ros_Gamma(3)= 0.21851380027664058511513169485832d+01
2222
2223	END SUBROUTINE Ros3
2224
2225	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2226
2227
2228	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2229	SUBROUTINE Ros4
2230	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2231	! L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES
2232	! L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
2233	!
2234	! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
2235	! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
2236	! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
2237	! SPRINGER-VERLAG (1990)
2238	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2239
2240	IMPLICIT NONE
2241
2242	rosMethod = RS4
2243	!~~~> Name of the method
2244	ros_Name = 'ROS-4'
2245	!~~~> Number of stages
2246	ros_S = 4
2247
2248	!~~~> The coefficient matrices A and C are strictly lower triangular.
2249	! The lower triangular (subdiagonal) elements are stored in row-wise order:
2250	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
2251	! The general mapping formula is:
2252	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
2253	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
2254
2255	ros_A(1) = 0.2000000000000000d+01
2256	ros_A(2) = 0.1867943637803922d+01
2257	ros_A(3) = 0.2344449711399156d+00
2258	ros_A(4) = ros_A(2)
2259	ros_A(5) = ros_A(3)
2260	ros_A(6) = 0.0D0
2261
2262	ros_C(1) =-0.7137615036412310d+01
2263	ros_C(2) = 0.2580708087951457d+01
2264	ros_C(3) = 0.6515950076447975d+00
2265	ros_C(4) =-0.2137148994382534d+01
2266	ros_C(5) =-0.3214669691237626d+00
2267	ros_C(6) =-0.6949742501781779d+00
2268	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
2269	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
2270	ros_NewF(1) = .TRUE.
2271	ros_NewF(2) = .TRUE.
2272	ros_NewF(3) = .TRUE.
2273	ros_NewF(4) = .FALSE.
2274	!~~~> M_i = Coefficients for new step solution
2275	ros_M(1) = 0.2255570073418735d+01
2276	ros_M(2) = 0.2870493262186792d+00
2277	ros_M(3) = 0.4353179431840180d+00
2278	ros_M(4) = 0.1093502252409163d+01
2279	!~~~> E_i = Coefficients for error estimator
2280	ros_E(1) =-0.2815431932141155d+00
2281	ros_E(2) =-0.7276199124938920d-01
2282	ros_E(3) =-0.1082196201495311d+00
2283	ros_E(4) =-0.1093502252409163d+01
2284	!~~~> ros_ELO = estimator of local order - the minimum between the
2285	! main and the embedded scheme orders plus 1
2286	ros_ELO = 4.0d0
2287	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
2288	ros_Alpha(1) = 0.D0
2289	ros_Alpha(2) = 0.1145640000000000d+01
2290	ros_Alpha(3) = 0.6552168638155900d+00
2291	ros_Alpha(4) = ros_Alpha(3)
2292	!~~~> Gamma_i = \sum_j gamma_{i,j}
2293	ros_Gamma(1) = 0.5728200000000000d+00
2294	ros_Gamma(2) =-0.1769193891319233d+01
2295	ros_Gamma(3) = 0.7592633437920482d+00
2296	ros_Gamma(4) =-0.1049021087100450d+00
2297
2298	END SUBROUTINE Ros4
2299
2300	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2301	SUBROUTINE Rodas3
2302	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2303	! --- A STIFFLY-STABLE METHOD, 4 stages, order 3
2304	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2305
2306	IMPLICIT NONE
2307
2308	rosMethod = RD3
2309	!~~~> Name of the method
2310	ros_Name = 'RODAS-3'
2311	!~~~> Number of stages
2312	ros_S = 4
2313
2314	!~~~> The coefficient matrices A and C are strictly lower triangular.
2315	! The lower triangular (subdiagonal) elements are stored in row-wise order:
2316	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
2317	! The general mapping formula is:
2318	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
2319	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
2320
2321	ros_A(1) = 0.0d+00
2322	ros_A(2) = 2.0d+00
2323	ros_A(3) = 0.0d+00
2324	ros_A(4) = 2.0d+00
2325	ros_A(5) = 0.0d+00
2326	ros_A(6) = 1.0d+00
2327
2328	ros_C(1) = 4.0d+00
2329	ros_C(2) = 1.0d+00
2330	ros_C(3) =-1.0d+00
2331	ros_C(4) = 1.0d+00
2332	ros_C(5) =-1.0d+00
2333	ros_C(6) =-(8.0d+00/3.0d+00)
2334
2335	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
2336	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
2337	ros_NewF(1) = .TRUE.
2338	ros_NewF(2) = .FALSE.
2339	ros_NewF(3) = .TRUE.
2340	ros_NewF(4) = .TRUE.
2341	!~~~> M_i = Coefficients for new step solution
2342	ros_M(1) = 2.0d+00
2343	ros_M(2) = 0.0d+00
2344	ros_M(3) = 1.0d+00
2345	ros_M(4) = 1.0d+00
2346	!~~~> E_i = Coefficients for error estimator
2347	ros_E(1) = 0.0d+00
2348	ros_E(2) = 0.0d+00
2349	ros_E(3) = 0.0d+00
2350	ros_E(4) = 1.0d+00
2351	!~~~> ros_ELO = estimator of local order - the minimum between the
2352	! main and the embedded scheme orders plus 1
2353	ros_ELO = 3.0d+00
2354	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
2355	ros_Alpha(1) = 0.0d+00
2356	ros_Alpha(2) = 0.0d+00
2357	ros_Alpha(3) = 1.0d+00
2358	ros_Alpha(4) = 1.0d+00
2359	!~~~> Gamma_i = \sum_j gamma_{i,j}
2360	ros_Gamma(1) = 0.5d+00
2361	ros_Gamma(2) = 1.5d+00
2362	ros_Gamma(3) = 0.0d+00
2363	ros_Gamma(4) = 0.0d+00
2364
2365	END SUBROUTINE Rodas3
2366
2367	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2368	SUBROUTINE Rodas4
2369	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2370	! STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES
2371	!
2372	! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
2373	! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
2374	! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
2375	! SPRINGER-VERLAG (1996)
2376	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2377
2378	IMPLICIT NONE
2379
2380	rosMethod = RD4
2381	!~~~> Name of the method
2382	ros_Name = 'RODAS-4'
2383	!~~~> Number of stages
2384	ros_S = 6
2385
2386	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
2387	ros_Alpha(1) = 0.000d0
2388	ros_Alpha(2) = 0.386d0
2389	ros_Alpha(3) = 0.210d0
2390	ros_Alpha(4) = 0.630d0
2391	ros_Alpha(5) = 1.000d0
2392	ros_Alpha(6) = 1.000d0
2393
2394	!~~~> Gamma_i = \sum_j gamma_{i,j}
2395	ros_Gamma(1) = 0.2500000000000000d+00
2396	ros_Gamma(2) =-0.1043000000000000d+00
2397	ros_Gamma(3) = 0.1035000000000000d+00
2398	ros_Gamma(4) =-0.3620000000000023d-01
2399	ros_Gamma(5) = 0.0d0
2400	ros_Gamma(6) = 0.0d0
2401
2402	!~~~> The coefficient matrices A and C are strictly lower triangular.
2403	! The lower triangular (subdiagonal) elements are stored in row-wise order:
2404	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
2405	! The general mapping formula is: A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
2406	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
2407
2408	ros_A(1) = 0.1544000000000000d+01
2409	ros_A(2) = 0.9466785280815826d+00
2410	ros_A(3) = 0.2557011698983284d+00
2411	ros_A(4) = 0.3314825187068521d+01
2412	ros_A(5) = 0.2896124015972201d+01
2413	ros_A(6) = 0.9986419139977817d+00
2414	ros_A(7) = 0.1221224509226641d+01
2415	ros_A(8) = 0.6019134481288629d+01
2416	ros_A(9) = 0.1253708332932087d+02
2417	ros_A(10) =-0.6878860361058950d+00
2418	ros_A(11) = ros_A(7)
2419	ros_A(12) = ros_A(8)
2420	ros_A(13) = ros_A(9)
2421	ros_A(14) = ros_A(10)
2422	ros_A(15) = 1.0d+00
2423
2424	ros_C(1) =-0.5668800000000000d+01
2425	ros_C(2) =-0.2430093356833875d+01
2426	ros_C(3) =-0.2063599157091915d+00
2427	ros_C(4) =-0.1073529058151375d+00
2428	ros_C(5) =-0.9594562251023355d+01
2429	ros_C(6) =-0.2047028614809616d+02
2430	ros_C(7) = 0.7496443313967647d+01
2431	ros_C(8) =-0.1024680431464352d+02
2432	ros_C(9) =-0.3399990352819905d+02
2433	ros_C(10) = 0.1170890893206160d+02
2434	ros_C(11) = 0.8083246795921522d+01
2435	ros_C(12) =-0.7981132988064893d+01
2436	ros_C(13) =-0.3152159432874371d+02
2437	ros_C(14) = 0.1631930543123136d+02
2438	ros_C(15) =-0.6058818238834054d+01
2439
2440	!~~~> M_i = Coefficients for new step solution
2441	ros_M(1) = ros_A(7)
2442	ros_M(2) = ros_A(8)
2443	ros_M(3) = ros_A(9)
2444	ros_M(4) = ros_A(10)
2445	ros_M(5) = 1.0d+00
2446	ros_M(6) = 1.0d+00
2447
2448	!~~~> E_i = Coefficients for error estimator
2449	ros_E(1) = 0.0d+00
2450	ros_E(2) = 0.0d+00
2451	ros_E(3) = 0.0d+00
2452	ros_E(4) = 0.0d+00
2453	ros_E(5) = 0.0d+00
2454	ros_E(6) = 1.0d+00
2455
2456	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
2457	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
2458	ros_NewF(1) = .TRUE.
2459	ros_NewF(2) = .TRUE.
2460	ros_NewF(3) = .TRUE.
2461	ros_NewF(4) = .TRUE.
2462	ros_NewF(5) = .TRUE.
2463	ros_NewF(6) = .TRUE.
2464
2465	!~~~> ros_ELO = estimator of local order - the minimum between the
2466	! main and the embedded scheme orders plus 1
2467	ros_ELO = 4.0d0
2468
2469	END SUBROUTINE Rodas4
2470
2471
2472	END SUBROUTINE RosenbrockADJ ! and its internal procedures
2473
2474
2475	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2476	SUBROUTINE FunTemplate( T, Y, Ydot )
2477	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2478	! Template for the ODE function call.
2479	! Updates the rate coefficients (and possibly the fixed species) at each call
2480	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2481
2482	!~~~> Input variables
2483	KPP_REAL, INTENT(IN) :: T, Y(NVAR)
2484	!~~~> Output variables
2485	KPP_REAL, INTENT(OUT) :: Ydot(NVAR)
2486	!~~~> Local variables
2487	KPP_REAL :: Told
2488
2489	Told = TIME
2490	TIME = T
2491	CALL Update_SUN()
2492	CALL Update_RCONST()
2493	CALL Fun( Y, FIX, RCONST, Ydot )
2494	TIME = Told
2495
2496	END SUBROUTINE FunTemplate
2497
2498
2499	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2500	SUBROUTINE JacTemplate( T, Y, Jcb )
2501	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2502	! Template for the ODE Jacobian call.
2503	! Updates the rate coefficients (and possibly the fixed species) at each call
2504	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2505
2506	!~~~> Input variables
2507	KPP_REAL :: T, Y(NVAR)
2508	!~~~> Output variables
2509	#ifdef FULL_ALGEBRA
2510	KPP_REAL :: JV(LU_NONZERO), Jcb(NVAR,NVAR)
2511	#else
2512	KPP_REAL :: Jcb(LU_NONZERO)
2513	#endif
2514	!~~~> Local variables
2515	KPP_REAL :: Told
2516	#ifdef FULL_ALGEBRA
2517	INTEGER :: i, j
2518	#endif
2519
2520	Told = TIME
2521	TIME = T
2522	CALL Update_SUN()
2523	CALL Update_RCONST()
2524	#ifdef FULL_ALGEBRA
2525	CALL Jac_SP(Y, FIX, RCONST, JV)
2526	DO j=1,NVAR
2527	DO i=1,NVAR
2528	Jcb(i,j) = 0.0_dp
2529	END DO
2530	END DO
2531	DO i=1,LU_NONZERO
2532	Jcb(LU_IROW(i),LU_ICOL(i)) = JV(i)
2533	END DO
2534	#else
2535	CALL Jac_SP( Y, FIX, RCONST, Jcb )
2536	#endif
2537	TIME = Told
2538
2539	END SUBROUTINE JacTemplate
2540
2541
2542
2543	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2544	SUBROUTINE HessTemplate( T, Y, Hes )
2545	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2546	! Template for the ODE Hessian call.
2547	! Updates the rate coefficients (and possibly the fixed species) at each call
2548	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2549
2550	!~~~> Input variables
2551	KPP_REAL, INTENT(IN) :: T, Y(NVAR)
2552	!~~~> Output variables
2553	KPP_REAL, INTENT(OUT) :: Hes(NHESS)
2554	!~~~> Local variables
2555	KPP_REAL :: Told
2556
2557	Told = TIME
2558	TIME = T
2559	CALL Update_SUN()
2560	CALL Update_RCONST()
2561	CALL Hessian( Y, FIX, RCONST, Hes )
2562	TIME = Told
2563
2564	END SUBROUTINE HessTemplate
2565
2566	END MODULE KPP_ROOT_Integrator
2567
2568
2569
2570

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

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/rosenbrock_adj.f90 @ 3512

Download in other formats: