Home

Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

rosenbrock_adj.f90 @ 4901

Last change on this file since 4901 was 2696, checked in by kanani, 6 years ago
Merge of branch palm4u into trunk
File size: 85.3 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 Fun_Template(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 Jac_Template(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:NVARj),1,Ynew,1)
987	END DO
988	Tau = T + ros_Alpha(istage)DirectionH
989	CALL Fun_Template(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:ioffset+NVAR),1)
995	Ystage(ioffset+1:ioffset+NVAR) = Ynew(1:NVAR)
996	END IF
997	!CALL WCOPY(NVAR,Fcn,1,K(ioffset+1:ioffset+NVAR),1)
998	K(ioffset+1:ioffset+NVAR) = Fcn(1:NVAR)
999	DO j = 1, istage-1
1000	HC = ros_C((istage-1)(istage-2)/2+j)/(DirectionH)
1001	CALL WAXPY(NVAR,HC,K(NVAR(j-1)+1:NVARj),1,K(ioffset+1:ioffset+NVAR),1)
1002	END DO
1003	IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN
1004	HG = DirectionHros_Gamma(istage)
1005	CALL WAXPY(NVAR,HG,dFdT,1,K(ioffset+1:ioffset+NVAR),1)
1006	END IF
1007	CALL ros_Solve('N', Ghimj, Pivot, K(ioffset+1:ioffset+NVAR))
1008
1009	END DO Stage
1010
1011
1012	!~~~> Compute the new solution
1013	CALL WCOPY(NVAR,Y,1,Ynew,1)
1014	DO j=1,ros_S
1015	CALL WAXPY(NVAR,ros_M(j),K(NVAR(j-1)+1:NVARj),1,Ynew,1)
1016	END DO
1017
1018	!~~~> Compute the error estimation
1019	CALL WSCAL(NVAR,ZERO,Yerr,1)
1020	DO j=1,ros_S
1021	CALL WAXPY(NVAR,ros_E(j),K(NVAR(j-1)+1:NVARj),1,Yerr,1)
1022	END DO
1023	Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol )
1024
1025	!~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax
1026	Fac = MIN(FacMax,MAX(FacMin,FacSafe/Err**(ONE/ros_ELO)))
1027	Hnew = H*Fac
1028
1029	!~~~> Check the error magnitude and adjust step size
1030	ISTATUS(Nstp) = ISTATUS(Nstp) + 1
1031	IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN !~~~> Accept step
1032	ISTATUS(Nacc) = ISTATUS(Nacc) + 1
1033	IF (AdjointType == Adj_discrete) THEN ! Save current state
1034	CALL ros_DPush( ros_S, T, H, Ystage, K, Ghimj, Pivot )
1035	ELSEIF ( (AdjointType == Adj_continuous) .OR. &
1036	(AdjointType == Adj_simple_continuous) ) THEN
1037	#ifdef FULL_ALGEBRA
1038	K = MATMUL(Jac0,Fcn0)
1039	#else
1040	CALL Jac_SP_Vec( Jac0, Fcn0, K(1:NVAR) )
1041	#endif
1042	IF (.NOT. Autonomous) THEN
1043	CALL WAXPY(NVAR,ONE,dFdT,1,K(1:NVAR),1)
1044	END IF
1045	CALL ros_CPush( T, H, Y, Fcn0, K(1:NVAR) )
1046	END IF
1047	CALL WCOPY(NVAR,Ynew,1,Y,1)
1048	T = T + Direction*H
1049	Hnew = MAX(Hmin,MIN(Hnew,Hmax))
1050	IF (RejectLastH) THEN ! No step size increase after a rejected step
1051	Hnew = MIN(Hnew,H)
1052	END IF
1053	RSTATUS(Nhexit) = H
1054	RSTATUS(Nhnew) = Hnew
1055	RSTATUS(Ntexit) = T
1056	RejectLastH = .FALSE.
1057	RejectMoreH = .FALSE.
1058	H = Hnew
1059	EXIT UntilAccepted ! EXIT THE LOOP: WHILE STEP NOT ACCEPTED
1060	ELSE !~~~> Reject step
1061	IF (RejectMoreH) THEN
1062	Hnew = H*FacRej
1063	END IF
1064	RejectMoreH = RejectLastH
1065	RejectLastH = .TRUE.
1066	H = Hnew
1067	IF (ISTATUS(Nacc) >= 1) THEN
1068	ISTATUS(Nrej) = ISTATUS(Nrej) + 1
1069	END IF
1070	END IF ! Err <= 1
1071
1072	END DO UntilAccepted
1073
1074	END DO TimeLoop
1075
1076	!~~~> Save last state: only needed for continuous adjoint
1077	IF ( (AdjointType == Adj_continuous) .OR. &
1078	(AdjointType == Adj_simple_continuous) ) THEN
1079	CALL Fun_Template(T,Y,Fcn0)
1080	ISTATUS(Nfun) = ISTATUS(Nfun) + 1
1081	CALL Jac_Template(T,Y,Jac0)
1082	ISTATUS(Njac) = ISTATUS(Njac) + 1
1083	#ifdef FULL_ALGEBRA
1084	K = MATMUL(Jac0,Fcn0)
1085	#else
1086	CALL Jac_SP_Vec( Jac0, Fcn0, K(1:NVAR) )
1087	#endif
1088	IF (.NOT. Autonomous) THEN
1089	CALL ros_FunTimeDerivative ( T, Roundoff, Y, &
1090	Fcn0, dFdT )
1091	CALL WAXPY(NVAR,ONE,dFdT,1,K(1:NVAR),1)
1092	END IF
1093	CALL ros_CPush( T, H, Y, Fcn0, K(1:NVAR) )
1094	!~~~> Deallocate stage buffer: only needed for discrete adjoint
1095	ELSEIF (AdjointType == Adj_discrete) THEN
1096	DEALLOCATE(Ystage, STAT=i)
1097	IF (i/=0) THEN
1098	PRINT*,'Deallocation of Ystage failed'
1099	STOP
1100	END IF
1101	END IF
1102
1103	!~~~> Succesful exit
1104	IERR = 1 !~~~> The integration was successful
1105
1106	END SUBROUTINE ros_FwdInt
1107
1108
1109	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1110	SUBROUTINE ros_DadjInt ( &
1111	NADJ, Lambda, &
1112	Tstart, Tend, T, &
1113	!~~~> Error indicator
1114	IERR )
1115	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1116	! Template for the implementation of a generic RosenbrockSOA method
1117	! defined by ros_S (no of stages)
1118	! and its coefficients ros_{A,C,M,E,Alpha,Gamma}
1119	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1120
1121	IMPLICIT NONE
1122
1123	!~~~> Input: the initial condition at Tstart; Output: the solution at T
1124	INTEGER, INTENT(IN) :: NADJ
1125	!~~~> First order adjoint
1126	KPP_REAL, INTENT(INOUT) :: Lambda(NVAR,NADJ)
1127	!!~~~> Input: integration interval
1128	KPP_REAL, INTENT(IN) :: Tstart,Tend
1129	!~~~> Output: time at which the solution is returned (T=Tend if success)
1130	KPP_REAL, INTENT(OUT) :: T
1131	!~~~> Output: Error indicator
1132	INTEGER, INTENT(OUT) :: IERR
1133	! ~~~~ Local variables
1134	KPP_REAL :: Ystage(NVARros_S), K(NVARros_S)
1135	KPP_REAL :: U(NVARros_S,NADJ), V(NVARros_S,NADJ)
1136	#ifdef FULL_ALGEBRA
1137	KPP_REAL, DIMENSION(NVAR,NVAR) :: Jac, dJdT, Ghimj
1138	#else
1139	KPP_REAL, DIMENSION(LU_NONZERO) :: Jac, dJdT, Ghimj
1140	#endif
1141	KPP_REAL :: Hes0(NHESS)
1142	KPP_REAL :: Tmp(NVAR), Tmp2(NVAR)
1143	KPP_REAL :: H, HC, HA, Tau
1144	INTEGER :: Pivot(NVAR), Direction
1145	INTEGER :: i, j, m, istage, ioffset, joffset
1146	!~~~> Local parameters
1147	KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
1148	KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5
1149	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1150
1151
1152
1153	IF (Tend >= Tstart) THEN
1154	Direction = +1
1155	ELSE
1156	Direction = -1
1157	END IF
1158
1159	!~~~> Time loop begins below
1160	TimeLoop: DO WHILE ( stack_ptr > 0 )
1161
1162	!~~~> Recover checkpoints for stage values and vectors
1163	CALL ros_DPop( ros_S, T, H, Ystage, K, Ghimj, Pivot )
1164
1165	! ISTATUS(Nstp) = ISTATUS(Nstp) + 1
1166
1167	!~~~> Compute LU decomposition
1168	IF (.NOT.SaveLU) THEN
1169	CALL Jac_Template(T,Ystage(1:NVAR),Ghimj)
1170	ISTATUS(Njac) = ISTATUS(Njac) + 1
1171	Tau = ONE/(DirectionHros_Gamma(1))
1172	#ifdef FULL_ALGEBRA
1173	Ghimj(1:NVAR,1:NVAR) = -Ghimj(1:NVAR,1:NVAR)
1174	DO i=1,NVAR
1175	Ghimj(i,i) = Ghimj(i,i)+Tau
1176	END DO
1177	#else
1178	CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1)
1179	DO i=1,NVAR
1180	Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+Tau
1181	END DO
1182	#endif
1183	CALL ros_Decomp( Ghimj, Pivot, j )
1184	END IF
1185
1186	!~~~> Compute Hessian at the beginning of the interval
1187	CALL Hess_Template(T,Ystage(1:NVAR),Hes0)
1188
1189	!~~~> Compute the stages
1190	Stage: DO istage = ros_S, 1, -1
1191
1192	!~~~> Current istage first entry
1193	ioffset = NVAR*(istage-1)
1194
1195	!~~~> Compute U
1196	DO m = 1,NADJ
1197	!CALL WCOPY(NVAR,Lambda(1:NVAR,m),1,U(ioffset+1:ioffset+NVAR,m),1)
1198	U(ioffset+1:ioffset+NVAR,m) = Lambda(1:NVAR,m)
1199	CALL WSCAL(NVAR,ros_M(istage),U(ioffset+1:ioffset+NVAR,m),1)
1200	END DO ! m=1:NADJ
1201	DO j = istage+1, ros_S
1202	joffset = NVAR*(j-1)
1203	HA = ros_A((j-1)*(j-2)/2+istage)
1204	HC = ros_C((j-1)(j-2)/2+istage)/(DirectionH)
1205	DO m = 1,NADJ
1206	CALL WAXPY(NVAR,HA,V(joffset+1:joffset+NVAR,m),1,U(ioffset+1:ioffset+NVAR,m),1)
1207	CALL WAXPY(NVAR,HC,U(joffset+1:joffset+NVAR,m),1,U(ioffset+1:ioffset+NVAR,m),1)
1208	END DO ! m=1:NADJ
1209	END DO
1210	DO m = 1,NADJ
1211	CALL ros_Solve('T', Ghimj, Pivot, U(ioffset+1:ioffset+NVAR,m))
1212	END DO ! m=1:NADJ
1213	!~~~> Compute V
1214	Tau = T + ros_Alpha(istage)DirectionH
1215	CALL Jac_Template(Tau,Ystage(ioffset+1:ioffset+NVAR),Jac)
1216	ISTATUS(Njac) = ISTATUS(Njac) + 1
1217	DO m = 1,NADJ
1218	#ifdef FULL_ALGEBRA
1219	V(ioffset+1:ioffset+NVAR,m) = MATMUL(TRANSPOSE(Jac),U(ioffset+1:ioffset+NVAR,m))
1220	#else
1221	CALL JacTR_SP_Vec(Jac,U(ioffset+1:ioffset+NVAR,m),V(ioffset+1:ioffset+NVAR,m))
1222	#endif
1223	END DO ! m=1:NADJ
1224
1225	END DO Stage
1226
1227	IF (.NOT.Autonomous) THEN
1228	!~~~> Compute the Jacobian derivative with respect to T.
1229	! Last "Jac" computed for stage 1
1230	CALL ros_JacTimeDerivative ( T, Roundoff, Ystage(1), Jac, dJdT )
1231	END IF
1232
1233	!~~~> Compute the new solution
1234
1235	!~~~> Compute Lambda
1236	DO istage=1,ros_S
1237	ioffset = NVAR*(istage-1)
1238	DO m = 1,NADJ
1239	! Add V_i
1240	CALL WAXPY(NVAR,ONE,V(ioffset+1:ioffset+NVAR,m),1,Lambda(1,m),1)
1241	! Add (H0xK_i)^T * U_i
1242	CALL HessTR_Vec ( Hes0, U(ioffset+1:ioffset+NVAR,m), &
1243	K(ioffset+1:ioffset+NVAR), Tmp )
1244	CALL WAXPY(NVAR,ONE,Tmp,1,Lambda(1:NVAR,m),1)
1245	END DO ! m=1:NADJ
1246	END DO
1247	! Add H * dJac_dT_0^T * \sum(gamma_i U_i)
1248	! Tmp holds sum gamma_i U_i
1249	IF (.NOT.Autonomous) THEN
1250	DO m = 1,NADJ
1251	Tmp(1:NVAR) = ZERO
1252	DO istage = 1, ros_S
1253	ioffset = NVAR*(istage-1)
1254	CALL WAXPY(NVAR,ros_Gamma(istage),U(ioffset+1:ioffset+NVAR,m),1,Tmp,1)
1255	END DO
1256	#ifdef FULL_ALGEBRA
1257	Tmp2 = MATMUL(TRANSPOSE(dJdT),Tmp)
1258	#else
1259	CALL JacTR_SP_Vec(dJdT,Tmp,Tmp2)
1260	#endif
1261	CALL WAXPY(NVAR,H,Tmp2,1,Lambda(1:NVAR,m),1)
1262	END DO ! m=1:NADJ
1263	END IF ! .NOT.Autonomous
1264
1265
1266	END DO TimeLoop
1267
1268	!~~~> Save last state
1269
1270	!~~~> Succesful exit
1271	IERR = 1 !~~~> The integration was successful
1272
1273	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1274	END SUBROUTINE ros_DadjInt
1275	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1276
1277
1278
1279	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1280	SUBROUTINE ros_CadjInt ( &
1281	NADJ, Y, &
1282	Tstart, Tend, T, &
1283	AbsTol_adj, RelTol_adj, &
1284	!~~~> Error indicator
1285	IERR )
1286	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1287	! Template for the implementation of a generic RosenbrockADJ method
1288	! defined by ros_S (no of stages)
1289	! and its coefficients ros_{A,C,M,E,Alpha,Gamma}
1290	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1291
1292	IMPLICIT NONE
1293
1294	!~~~> Input: the initial condition at Tstart; Output: the solution at T
1295	INTEGER, INTENT(IN) :: NADJ
1296	KPP_REAL, INTENT(INOUT) :: Y(NVAR,NADJ)
1297	!~~~> Input: integration interval
1298	KPP_REAL, INTENT(IN) :: Tstart,Tend
1299	!~~~> Input: adjoint tolerances
1300	KPP_REAL, INTENT(IN) :: AbsTol_adj(NVAR,NADJ), RelTol_adj(NVAR,NADJ)
1301	!~~~> Output: time at which the solution is returned (T=Tend if success)
1302	KPP_REAL, INTENT(OUT) :: T
1303	!~~~> Output: Error indicator
1304	INTEGER, INTENT(OUT) :: IERR
1305	! ~~~~ Local variables
1306	KPP_REAL :: Y0(NVAR)
1307	KPP_REAL :: Ynew(NVAR,NADJ), Fcn0(NVAR,NADJ), Fcn(NVAR,NADJ)
1308	KPP_REAL :: K(NVAR*ros_S,NADJ), dFdT(NVAR,NADJ)
1309	#ifdef FULL_ALGEBRA
1310	KPP_REAL, DIMENSION(NVAR,NVAR) :: Jac0, Ghimj, Jac, dJdT
1311	#else
1312	KPP_REAL, DIMENSION(LU_NONZERO) :: Jac0, Ghimj, Jac, dJdT
1313	#endif
1314	KPP_REAL :: H, Hnew, HC, HG, Fac, Tau
1315	KPP_REAL :: Err, Yerr(NVAR,NADJ)
1316	INTEGER :: Pivot(NVAR), Direction, ioffset, j, istage, iadj
1317	LOGICAL :: RejectLastH, RejectMoreH, Singular
1318	!~~~> Local parameters
1319	KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
1320	KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5
1321	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1322
1323
1324	!~~~> Initial preparations
1325	T = Tstart
1326	RSTATUS(Nhexit) = 0.0_dp
1327	H = MIN( MAX(ABS(Hmin),ABS(Hstart)) , ABS(Hmax) )
1328	IF (ABS(H) <= 10.0_dp*Roundoff) H = DeltaMin
1329
1330	IF (Tend >= Tstart) THEN
1331	Direction = +1
1332	ELSE
1333	Direction = -1
1334	END IF
1335	H = Direction*H
1336
1337	RejectLastH=.FALSE.
1338	RejectMoreH=.FALSE.
1339
1340	!~~~> Time loop begins below
1341
1342	TimeLoop: DO WHILE ( (Direction > 0).AND.((T-Tend)+Roundoff <= ZERO) &
1343	.OR. (Direction < 0).AND.((Tend-T)+Roundoff <= ZERO) )
1344
1345	IF ( ISTATUS(Nstp) > Max_no_steps ) THEN ! Too many steps
1346	CALL ros_ErrorMsg(-6,T,H,IERR)
1347	RETURN
1348	END IF
1349	IF ( ((T+0.1d0*H) == T).OR.(H <= Roundoff) ) THEN ! Step size too small
1350	CALL ros_ErrorMsg(-7,T,H,IERR)
1351	RETURN
1352	END IF
1353
1354	!~~~> Limit H if necessary to avoid going beyond Tend
1355	RSTATUS(Nhexit) = H
1356	H = MIN(H,ABS(Tend-T))
1357
1358	!~~~> Interpolate forward solution
1359	CALL ros_cadj_Y( T, Y0 )
1360	!~~~> Compute the Jacobian at current time
1361	CALL Jac_Template(T, Y0, Jac0)
1362	ISTATUS(Njac) = ISTATUS(Njac) + 1
1363
1364	!~~~> Compute the function derivative with respect to T
1365	IF (.NOT.Autonomous) THEN
1366	CALL ros_JacTimeDerivative ( T, Roundoff, Y0, &
1367	Jac0, dJdT )
1368	DO iadj = 1, NADJ
1369	#ifdef FULL_ALGEBRA
1370	dFdT(1:NVAR,iadj) = MATMUL(TRANSPOSE(dJdT),Y(1:NVAR,iadj))
1371	#else
1372	CALL JacTR_SP_Vec(dJdT,Y(1:NVAR,iadj),dFdT(1:NVAR,iadj))
1373	#endif
1374	CALL WSCAL(NVAR,(-ONE),dFdT(1:NVAR,iadj),1)
1375	END DO
1376	END IF
1377
1378	!~~~> Ydot = -J^T*Y
1379	#ifdef FULL_ALGEBRA
1380	Jac0(1:NVAR,1:NVAR) = -Jac0(1:NVAR,1:NVAR)
1381	#else
1382	CALL WSCAL(LU_NONZERO,(-ONE),Jac0,1)
1383	#endif
1384	DO iadj = 1, NADJ
1385	#ifdef FULL_ALGEBRA
1386	Fcn0(1:NVAR,iadj) = MATMUL(TRANSPOSE(Jac0),Y(1:NVAR,iadj))
1387	#else
1388	CALL JacTR_SP_Vec(Jac0,Y(1:NVAR,iadj),Fcn0(1:NVAR,iadj))
1389	#endif
1390	END DO
1391
1392	!~~~> Repeat step calculation until current step accepted
1393	UntilAccepted: DO
1394
1395	CALL ros_PrepareMatrix(H,Direction,ros_Gamma(1), &
1396	Jac0,Ghimj,Pivot,Singular)
1397	IF (Singular) THEN ! More than 5 consecutive failed decompositions
1398	CALL ros_ErrorMsg(-8,T,H,IERR)
1399	RETURN
1400	END IF
1401
1402	!~~~> Compute the stages
1403	Stage: DO istage = 1, ros_S
1404
1405	! Current istage offset. Current istage vector is K(ioffset+1:ioffset+NVAR)
1406	ioffset = NVAR*(istage-1)
1407
1408	! For the 1st istage the function has been computed previously
1409	IF ( istage == 1 ) THEN
1410	DO iadj = 1, NADJ
1411	!CALL WCOPY(NVAR,Fcn0(1,iadj),1,Fcn(1,iadj),1)
1412	Fcn(1:NVAR,iadj) = Fcn0(1:NVAR,iadj)
1413	END DO
1414	! istage>1 and a new function evaluation is needed at the current istage
1415	ELSEIF ( ros_NewF(istage) ) THEN
1416	CALL WCOPY(NVAR*NADJ,Y,1,Ynew,1)
1417	DO j = 1, istage-1
1418	DO iadj = 1, NADJ
1419	CALL WAXPY(NVAR,ros_A((istage-1)*(istage-2)/2+j), &
1420	K(NVAR(j-1)+1:NVARj,iadj),1,Ynew(1,iadj),1)
1421	END DO
1422	END DO
1423	Tau = T + ros_Alpha(istage)DirectionH
1424	! CALL Fun_Template(Tau,Ynew,Fcn)
1425	! ISTATUS(Nfun) = ISTATUS(Nfun) + 1
1426	CALL ros_cadj_Y( Tau, Y0 )
1427	CALL Jac_Template(Tau, Y0, Jac)
1428	ISTATUS(Njac) = ISTATUS(Njac) + 1
1429	#ifdef FULL_ALGEBRA
1430	Jac(1:NVAR,1:NVAR) = -Jac(1:NVAR,1:NVAR)
1431	#else
1432	CALL WSCAL(LU_NONZERO,(-ONE),Jac,1)
1433	#endif
1434	DO iadj = 1, NADJ
1435	#ifdef FULL_ALGEBRA
1436	Fcn(1:NVAR,iadj) = MATMUL(TRANSPOSE(Jac),Ynew(1:NVAR,iadj))
1437	#else
1438	CALL JacTR_SP_Vec(Jac,Ynew(1:NVAR,iadj),Fcn(1:NVAR,iadj))
1439	#endif
1440	!CALL WSCAL(NVAR,(-ONE),Fcn(1,iadj),1)
1441	END DO
1442	END IF ! if istage == 1 elseif ros_NewF(istage)
1443
1444	DO iadj = 1, NADJ
1445	!CALL WCOPY(NVAR,Fcn(1,iadj),1,K(ioffset+1,iadj),1)
1446	K(ioffset+1:ioffset+NVAR,iadj) = Fcn(1:NVAR,iadj)
1447	END DO
1448	DO j = 1, istage-1
1449	HC = ros_C((istage-1)(istage-2)/2+j)/(DirectionH)
1450	DO iadj = 1, NADJ
1451	CALL WAXPY(NVAR,HC,K(NVAR(j-1)+1:NVARj,iadj),1, &
1452	K(ioffset+1:ioffset+NVAR,iadj),1)
1453	END DO
1454	END DO
1455	IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN
1456	HG = DirectionHros_Gamma(istage)
1457	DO iadj = 1, NADJ
1458	CALL WAXPY(NVAR,HG,dFdT(1:NVAR,iadj),1,K(ioffset+1:ioffset+NVAR,iadj),1)
1459	END DO
1460	END IF
1461	DO iadj = 1, NADJ
1462	CALL ros_Solve('T', Ghimj, Pivot, K(ioffset+1:ioffset+NVAR,iadj))
1463	END DO
1464
1465	END DO Stage
1466
1467
1468	!~~~> Compute the new solution
1469	DO iadj = 1, NADJ
1470	!CALL WCOPY(NVAR,Y(1:NVAR,iadj),1,Ynew(1,iadj),1)
1471	Ynew(1:NVAR,iadj) = Y(1:NVAR,iadj)
1472	DO j=1,ros_S
1473	CALL WAXPY(NVAR,ros_M(j),K(NVAR(j-1)+1:NVARj,iadj),1,Ynew(1:NVAR,iadj),1)
1474	END DO
1475	END DO
1476
1477	!~~~> Compute the error estimation
1478	CALL WSCAL(NVAR*NADJ,ZERO,Yerr,1)
1479	DO j=1,ros_S
1480	DO iadj = 1, NADJ
1481	CALL WAXPY(NVAR,ros_E(j),K(NVAR(j-1)+1:NVARj,iadj),1,Yerr(1:NVAR,iadj),1)
1482	END DO
1483	END DO
1484	!~~~> Max error among all adjoint components
1485	iadj = 1
1486	Err = ros_ErrorNorm ( Y(1:NVAR,iadj), Ynew(1,iadj), Yerr(1,iadj), &
1487	AbsTol_adj(1,iadj), RelTol_adj(1,iadj), VectorTol )
1488
1489	!~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax
1490	Fac = MIN(FacMax,MAX(FacMin,FacSafe/Err**(ONE/ros_ELO)))
1491	Hnew = H*Fac
1492
1493	!~~~> Check the error magnitude and adjust step size
1494	! ISTATUS(Nstp) = ISTATUS(Nstp) + 1
1495	IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN !~~~> Accept step
1496	ISTATUS(Nacc) = ISTATUS(Nacc) + 1
1497	CALL WCOPY(NVAR*NADJ,Ynew,1,Y,1)
1498	T = T + Direction*H
1499	Hnew = MAX(Hmin,MIN(Hnew,Hmax))
1500	IF (RejectLastH) THEN ! No step size increase after a rejected step
1501	Hnew = MIN(Hnew,H)
1502	END IF
1503	RSTATUS(Nhexit) = H
1504	RSTATUS(Nhnew) = Hnew
1505	RSTATUS(Ntexit) = T
1506	RejectLastH = .FALSE.
1507	RejectMoreH = .FALSE.
1508	H = Hnew
1509	EXIT UntilAccepted ! EXIT THE LOOP: WHILE STEP NOT ACCEPTED
1510	ELSE !~~~> Reject step
1511	IF (RejectMoreH) THEN
1512	Hnew = H*FacRej
1513	END IF
1514	RejectMoreH = RejectLastH
1515	RejectLastH = .TRUE.
1516	H = Hnew
1517	IF (ISTATUS(Nacc) >= 1) THEN
1518	ISTATUS(Nrej) = ISTATUS(Nrej) + 1
1519	END IF
1520	END IF ! Err <= 1
1521
1522	END DO UntilAccepted
1523
1524	END DO TimeLoop
1525
1526	!~~~> Succesful exit
1527	IERR = 1 !~~~> The integration was successful
1528
1529	END SUBROUTINE ros_CadjInt
1530
1531
1532	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1533	SUBROUTINE ros_SimpleCadjInt ( &
1534	NADJ, Y, &
1535	Tstart, Tend, T, &
1536	!~~~> Error indicator
1537	IERR )
1538	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1539	! Template for the implementation of a generic RosenbrockADJ method
1540	! defined by ros_S (no of stages)
1541	! and its coefficients ros_{A,C,M,E,Alpha,Gamma}
1542	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1543
1544	IMPLICIT NONE
1545
1546	!~~~> Input: the initial condition at Tstart; Output: the solution at T
1547	INTEGER, INTENT(IN) :: NADJ
1548	KPP_REAL, INTENT(INOUT) :: Y(NVAR,NADJ)
1549	!~~~> Input: integration interval
1550	KPP_REAL, INTENT(IN) :: Tstart,Tend
1551	!~~~> Output: time at which the solution is returned (T=Tend if success)
1552	KPP_REAL, INTENT(OUT) :: T
1553	!~~~> Output: Error indicator
1554	INTEGER, INTENT(OUT) :: IERR
1555	! ~~~~ Local variables
1556	KPP_REAL :: Y0(NVAR)
1557	KPP_REAL :: Ynew(NVAR,NADJ), Fcn0(NVAR,NADJ), Fcn(NVAR,NADJ)
1558	KPP_REAL :: K(NVAR*ros_S,NADJ), dFdT(NVAR,NADJ)
1559	#ifdef FULL_ALGEBRA
1560	KPP_REAL,DIMENSION(NVAR,NVAR) :: Jac0, Ghimj, Jac, dJdT
1561	#else
1562	KPP_REAL,DIMENSION(LU_NONZERO) :: Jac0, Ghimj, Jac, dJdT
1563	#endif
1564	KPP_REAL :: H, HC, HG, Tau
1565	KPP_REAL :: ghinv
1566	INTEGER :: Pivot(NVAR), Direction, ioffset, i, j, istage, iadj
1567	INTEGER :: istack
1568	!~~~> Local parameters
1569	KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
1570	KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5
1571	!~~~> Locally called functions
1572	! KPP_REAL WLAMCH
1573	! EXTERNAL WLAMCH
1574	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1575
1576
1577	!~~~> INITIAL PREPARATIONS
1578
1579	IF (Tend >= Tstart) THEN
1580	Direction = -1
1581	ELSE
1582	Direction = +1
1583	END IF
1584
1585	!~~~> Time loop begins below
1586	TimeLoop: DO istack = stack_ptr,2,-1
1587
1588	T = chk_T(istack)
1589	H = chk_H(istack-1)
1590	!CALL WCOPY(NVAR,chk_Y(1,istack),1,Y0,1)
1591	Y0(1:NVAR) = chk_Y(1:NVAR,istack)
1592
1593	!~~~> Compute the Jacobian at current time
1594	CALL Jac_Template(T, Y0, Jac0)
1595	ISTATUS(Njac) = ISTATUS(Njac) + 1
1596
1597	!~~~> Compute the function derivative with respect to T
1598	IF (.NOT.Autonomous) THEN
1599	CALL ros_JacTimeDerivative ( T, Roundoff, Y0, &
1600	Jac0, dJdT )
1601	DO iadj = 1, NADJ
1602	#ifdef FULL_ALGEBRA
1603	dFdT(1:NVAR,iadj) = MATMUL(TRANSPOSE(dJdT),Y(1:NVAR,iadj))
1604	#else
1605	CALL JacTR_SP_Vec(dJdT,Y(1:NVAR,iadj),dFdT(1:NVAR,iadj))
1606	#endif
1607	CALL WSCAL(NVAR,(-ONE),dFdT(1:NVAR,iadj),1)
1608	END DO
1609	END IF
1610
1611	!~~~> Ydot = -J^T*Y
1612	#ifdef FULL_ALGEBRA
1613	Jac0(1:NVAR,1:NVAR) = -Jac0(1:NVAR,1:NVAR)
1614	#else
1615	CALL WSCAL(LU_NONZERO,(-ONE),Jac0,1)
1616	#endif
1617	DO iadj = 1, NADJ
1618	#ifdef FULL_ALGEBRA
1619	Fcn0(1:NVAR,iadj) = MATMUL(TRANSPOSE(Jac0),Y(1:NVAR,iadj))
1620	#else
1621	CALL JacTR_SP_Vec(Jac0,Y(1:NVAR,iadj),Fcn0(1:NVAR,iadj))
1622	#endif
1623	END DO
1624
1625	!~~~> Construct Ghimj = 1/(H*ham) - Jac0
1626	ghinv = ONE/(DirectionHros_Gamma(1))
1627	#ifdef FULL_ALGEBRA
1628	Ghimj(1:NVAR,1:NVAR) = -Jac0(1:NVAR,1:NVAR)
1629	DO i=1,NVAR
1630	Ghimj(i,i) = Ghimj(i,i)+ghinv
1631	END DO
1632	#else
1633	CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1)
1634	CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1)
1635	DO i=1,NVAR
1636	Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+ghinv
1637	END DO
1638	#endif
1639	!~~~> Compute LU decomposition
1640	CALL ros_Decomp( Ghimj, Pivot, j )
1641	IF (j /= 0) THEN
1642	CALL ros_ErrorMsg(-8,T,H,IERR)
1643	PRINT*,' The matrix is singular !'
1644	STOP
1645	END IF
1646
1647	!~~~> Compute the stages
1648	Stage: DO istage = 1, ros_S
1649
1650	! Current istage offset. Current istage vector is K(ioffset+1:ioffset+NVAR)
1651	ioffset = NVAR*(istage-1)
1652
1653	! For the 1st istage the function has been computed previously
1654	IF ( istage == 1 ) THEN
1655	DO iadj = 1, NADJ
1656	CALL WCOPY(NVAR,Fcn0(1,iadj),1,Fcn(1,iadj),1)
1657	END DO
1658	! istage>1 and a new function evaluation is needed at the current istage
1659	ELSEIF ( ros_NewF(istage) ) THEN
1660	CALL WCOPY(NVAR*NADJ,Y,1,Ynew,1)
1661	DO j = 1, istage-1
1662	DO iadj = 1, NADJ
1663	CALL WAXPY(NVAR,ros_A((istage-1)*(istage-2)/2+j), &
1664	K(NVAR(j-1)+1:NVARj,iadj),1,Ynew(1:NVAR,iadj),1)
1665	END DO
1666	END DO
1667	Tau = T + ros_Alpha(istage)DirectionH
1668	CALL ros_Hermite3( chk_T(istack-1), chk_T(istack), Tau, &
1669	chk_Y(1:NVAR,istack-1), chk_Y(1:NVAR,istack), &
1670	chk_dY(1:NVAR,istack-1), chk_dY(1:NVAR,istack), Y0 )
1671	CALL Jac_Template(Tau, Y0, Jac)
1672	ISTATUS(Njac) = ISTATUS(Njac) + 1
1673	#ifdef FULL_ALGEBRA
1674	Jac(1:NVAR,1:NVAR) = -Jac(1:NVAR,1:NVAR)
1675	#else
1676	CALL WSCAL(LU_NONZERO,(-ONE),Jac,1)
1677	#endif
1678	DO iadj = 1, NADJ
1679	#ifdef FULL_ALGEBRA
1680	Fcn(1:NVAR,iadj) = MATMUL(TRANSPOSE(Jac),Ynew(1:NVAR,iadj))
1681	#else
1682	CALL JacTR_SP_Vec(Jac,Ynew(1:NVAR,iadj),Fcn(1:NVAR,iadj))
1683	#endif
1684	END DO
1685	END IF ! if istage == 1 elseif ros_NewF(istage)
1686
1687	DO iadj = 1, NADJ
1688	CALL WCOPY(NVAR,Fcn(1,iadj),1,K(ioffset+1,iadj),1)
1689	END DO
1690	DO j = 1, istage-1
1691	HC = ros_C((istage-1)(istage-2)/2+j)/(DirectionH)
1692	DO iadj = 1, NADJ
1693	CALL WAXPY(NVAR,HC,K(NVAR(j-1)+1:NVARj,iadj),1, &
1694	K(ioffset+1:ioffset+NVAR,iadj),1)
1695	END DO
1696	END DO
1697	IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN
1698	HG = DirectionHros_Gamma(istage)
1699	DO iadj = 1, NADJ
1700	CALL WAXPY(NVAR,HG,dFdT(1:NVAR,iadj),1,K(ioffset+1:ioffset+NVAR,iadj),1)
1701	END DO
1702	END IF
1703	DO iadj = 1, NADJ
1704	CALL ros_Solve('T', Ghimj, Pivot, K(ioffset+1:ioffset+NVAR,iadj))
1705	END DO
1706
1707	END DO Stage
1708
1709
1710	!~~~> Compute the new solution
1711	DO iadj = 1, NADJ
1712	DO j=1,ros_S
1713	CALL WAXPY(NVAR,ros_M(j),K(NVAR(j-1)+1:NVARj,iadj),1,Y(1:NVAR,iadj),1)
1714	END DO
1715	END DO
1716
1717	END DO TimeLoop
1718
1719	!~~~> Succesful exit
1720	IERR = 1 !~~~> The integration was successful
1721
1722	END SUBROUTINE ros_SimpleCadjInt
1723
1724	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1725	KPP_REAL FUNCTION ros_ErrorNorm ( Y, Ynew, Yerr, &
1726	AbsTol, RelTol, VectorTol )
1727	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1728	!~~~> Computes the "scaled norm" of the error vector Yerr
1729	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1730	IMPLICIT NONE
1731
1732	! Input arguments
1733	KPP_REAL, INTENT(IN) :: Y(NVAR), Ynew(NVAR), &
1734	Yerr(NVAR), AbsTol(NVAR), RelTol(NVAR)
1735	LOGICAL, INTENT(IN) :: VectorTol
1736	! Local variables
1737	KPP_REAL :: Err, Scale, Ymax
1738	INTEGER :: i
1739
1740	Err = ZERO
1741	DO i=1,NVAR
1742	Ymax = MAX(ABS(Y(i)),ABS(Ynew(i)))
1743	IF (VectorTol) THEN
1744	Scale = AbsTol(i)+RelTol(i)*Ymax
1745	ELSE
1746	Scale = AbsTol(1)+RelTol(1)*Ymax
1747	END IF
1748	Err = Err+(Yerr(i)/Scale)**2
1749	END DO
1750	Err = SQRT(Err/NVAR)
1751
1752	ros_ErrorNorm = MAX(Err,1.0d-10)
1753
1754	END FUNCTION ros_ErrorNorm
1755
1756
1757	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1758	SUBROUTINE ros_FunTimeDerivative ( T, Roundoff, Y, Fcn0, dFdT )
1759	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1760	!~~~> The time partial derivative of the function by finite differences
1761	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1762	IMPLICIT NONE
1763
1764	!~~~> Input arguments
1765	KPP_REAL, INTENT(IN) :: T, Roundoff, Y(NVAR), Fcn0(NVAR)
1766	!~~~> Output arguments
1767	KPP_REAL, INTENT(OUT) :: dFdT(NVAR)
1768	!~~~> Local variables
1769	KPP_REAL :: Delta
1770	KPP_REAL, PARAMETER :: ONE = 1.0d0, DeltaMin = 1.0d-6
1771
1772	Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T))
1773	CALL Fun_Template(T+Delta,Y,dFdT)
1774	ISTATUS(Nfun) = ISTATUS(Nfun) + 1
1775	CALL WAXPY(NVAR,(-ONE),Fcn0,1,dFdT,1)
1776	CALL WSCAL(NVAR,(ONE/Delta),dFdT,1)
1777
1778	END SUBROUTINE ros_FunTimeDerivative
1779
1780
1781	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1782	SUBROUTINE ros_JacTimeDerivative ( T, Roundoff, Y, &
1783	Jac0, dJdT )
1784	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1785	!~~~> The time partial derivative of the Jacobian by finite differences
1786	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1787	IMPLICIT NONE
1788
1789	!~~~> Arguments
1790	KPP_REAL, INTENT(IN) :: T, Roundoff, Y(NVAR)
1791	#ifdef FULL_ALGEBRA
1792	KPP_REAL, INTENT(IN) :: Jac0(NVAR,NVAR)
1793	KPP_REAL, INTENT(OUT) :: dJdT(NVAR,NVAR)
1794	#else
1795	KPP_REAL, INTENT(IN) :: Jac0(LU_NONZERO)
1796	KPP_REAL, INTENT(OUT) :: dJdT(LU_NONZERO)
1797	#endif
1798	!~~~> Local variables
1799	KPP_REAL :: Delta
1800
1801	Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T))
1802	CALL Jac_Template(T+Delta,Y,dJdT)
1803	ISTATUS(Njac) = ISTATUS(Njac) + 1
1804	#ifdef FULL_ALGEBRA
1805	CALL WAXPY(NVAR*NVAR,(-ONE),Jac0,1,dJdT,1)
1806	CALL WSCAL(NVAR*NVAR,(ONE/Delta),dJdT,1)
1807	#else
1808	CALL WAXPY(LU_NONZERO,(-ONE),Jac0,1,dJdT,1)
1809	CALL WSCAL(LU_NONZERO,(ONE/Delta),dJdT,1)
1810	#endif
1811
1812	END SUBROUTINE ros_JacTimeDerivative
1813
1814
1815
1816	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1817	SUBROUTINE ros_PrepareMatrix ( H, Direction, gam, &
1818	Jac0, Ghimj, Pivot, Singular )
1819	! --- --- --- --- --- --- --- --- --- --- --- --- ---
1820	! Prepares the LHS matrix for stage calculations
1821	! 1. Construct Ghimj = 1/(H*gam) - Jac0
1822	! "(Gamma H) Inverse Minus Jacobian"
1823	! 2. Repeat LU decomposition of Ghimj until successful.
1824	! -half the step size if LU decomposition fails and retry
1825	! -exit after 5 consecutive fails
1826	! --- --- --- --- --- --- --- --- --- --- --- --- ---
1827	IMPLICIT NONE
1828
1829	!~~~> Input arguments
1830	#ifdef FULL_ALGEBRA
1831	KPP_REAL, INTENT(IN) :: Jac0(NVAR,NVAR)
1832	#else
1833	KPP_REAL, INTENT(IN) :: Jac0(LU_NONZERO)
1834	#endif
1835	KPP_REAL, INTENT(IN) :: gam
1836	INTEGER, INTENT(IN) :: Direction
1837	!~~~> Output arguments
1838	#ifdef FULL_ALGEBRA
1839	KPP_REAL, INTENT(OUT) :: Ghimj(NVAR,NVAR)
1840	#else
1841	KPP_REAL, INTENT(OUT) :: Ghimj(LU_NONZERO)
1842	#endif
1843	LOGICAL, INTENT(OUT) :: Singular
1844	INTEGER, INTENT(OUT) :: Pivot(NVAR)
1845	!~~~> Inout arguments
1846	KPP_REAL, INTENT(INOUT) :: H ! step size is decreased when LU fails
1847	!~~~> Local variables
1848	INTEGER :: i, ising, Nconsecutive
1849	KPP_REAL :: ghinv
1850	KPP_REAL, PARAMETER :: ONE = 1.0_dp, HALF = 0.5_dp
1851
1852	Nconsecutive = 0
1853	Singular = .TRUE.
1854
1855	DO WHILE (Singular)
1856
1857	!~~~> Construct Ghimj = 1/(H*gam) - Jac0
1858	#ifdef FULL_ALGEBRA
1859	CALL WCOPY(NVAR*NVAR,Jac0,1,Ghimj,1)
1860	CALL WSCAL(NVAR*NVAR,(-ONE),Ghimj,1)
1861	ghinv = ONE/(DirectionHgam)
1862	DO i=1,NVAR
1863	Ghimj(i,i) = Ghimj(i,i)+ghinv
1864	END DO
1865	#else
1866	CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1)
1867	CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1)
1868	ghinv = ONE/(DirectionHgam)
1869	DO i=1,NVAR
1870	Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+ghinv
1871	END DO
1872	#endif
1873	!~~~> Compute LU decomposition
1874	CALL ros_Decomp( Ghimj, Pivot, ising )
1875	IF (ising == 0) THEN
1876	!~~~> If successful done
1877	Singular = .FALSE.
1878	ELSE ! ising .ne. 0
1879	!~~~> If unsuccessful half the step size; if 5 consecutive fails then return
1880	ISTATUS(Nsng) = ISTATUS(Nsng) + 1
1881	Nconsecutive = Nconsecutive+1
1882	Singular = .TRUE.
1883	PRINT*,'Warning: LU Decomposition returned ising = ',ising
1884	IF (Nconsecutive <= 5) THEN ! Less than 5 consecutive failed decompositions
1885	H = H*HALF
1886	ELSE ! More than 5 consecutive failed decompositions
1887	RETURN
1888	END IF ! Nconsecutive
1889	END IF ! ising
1890
1891	END DO ! WHILE Singular
1892
1893	END SUBROUTINE ros_PrepareMatrix
1894
1895
1896	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1897	SUBROUTINE ros_Decomp( A, Pivot, ising )
1898	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1899	! Template for the LU decomposition
1900	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1901	IMPLICIT NONE
1902	!~~~> Inout variables
1903	#ifdef FULL_ALGEBRA
1904	KPP_REAL, INTENT(INOUT) :: A(NVAR,NVAR)
1905	#else
1906	KPP_REAL, INTENT(INOUT) :: A(LU_NONZERO)
1907	#endif
1908	!~~~> Output variables
1909	INTEGER, INTENT(OUT) :: Pivot(NVAR), ising
1910
1911	#ifdef FULL_ALGEBRA
1912	CALL DGETRF( NVAR, NVAR, A, NVAR, Pivot, ising )
1913	#else
1914	CALL KppDecomp ( A, ising )
1915	Pivot(1) = 1
1916	#endif
1917	ISTATUS(Ndec) = ISTATUS(Ndec) + 1
1918
1919	END SUBROUTINE ros_Decomp
1920
1921
1922	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1923	SUBROUTINE ros_Solve( How, A, Pivot, b )
1924	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1925	! Template for the forward/backward substitution (using pre-computed LU decomposition)
1926	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1927	IMPLICIT NONE
1928	!~~~> Input variables
1929	CHARACTER, INTENT(IN) :: How
1930	#ifdef FULL_ALGEBRA
1931	KPP_REAL, INTENT(IN) :: A(NVAR,NVAR)
1932	#else
1933	KPP_REAL, INTENT(IN) :: A(LU_NONZERO)
1934	#endif
1935	INTEGER, INTENT(IN) :: Pivot(NVAR)
1936	!~~~> InOut variables
1937	KPP_REAL, INTENT(INOUT) :: b(NVAR)
1938
1939	SELECT CASE (How)
1940	CASE ('N')
1941	#ifdef FULL_ALGEBRA
1942	CALL DGETRS( 'N', NVAR , 1, A, NVAR, Pivot, b, NVAR, 0 )
1943	#else
1944	CALL KppSolve( A, b )
1945	#endif
1946	CASE ('T')
1947	#ifdef FULL_ALGEBRA
1948	CALL DGETRS( 'T', NVAR , 1, A, NVAR, Pivot, b, NVAR, 0 )
1949	#else
1950	CALL KppSolveTR( A, b, b )
1951	#endif
1952	CASE DEFAULT
1953	PRINT*,'Error: unknown argument in ros_Solve: How=',How
1954	STOP
1955	END SELECT
1956	ISTATUS(Nsol) = ISTATUS(Nsol) + 1
1957
1958	END SUBROUTINE ros_Solve
1959
1960
1961
1962	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1963	SUBROUTINE ros_cadj_Y( T, Y )
1964	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1965	! Finds the solution Y at T by interpolating the stored forward trajectory
1966	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1967	IMPLICIT NONE
1968	!~~~> Input variables
1969	KPP_REAL, INTENT(IN) :: T
1970	!~~~> Output variables
1971	KPP_REAL, INTENT(OUT) :: Y(NVAR)
1972	!~~~> Local variables
1973	INTEGER :: i
1974	KPP_REAL, PARAMETER :: ONE = 1.0d0
1975
1976	! chk_H, chk_T, chk_Y, chk_dY, chk_d2Y
1977
1978	IF( (T < chk_T(1)).OR.(T> chk_T(stack_ptr)) ) THEN
1979	PRINT*,'Cannot locate solution at T = ',T
1980	PRINT*,'Stored trajectory is between Tstart = ',chk_T(1)
1981	PRINT*,' and Tend = ',chk_T(stack_ptr)
1982	STOP
1983	END IF
1984	DO i = 1, stack_ptr-1
1985	IF( (T>= chk_T(i)).AND.(T<= chk_T(i+1)) ) EXIT
1986	END DO
1987
1988
1989	! IF (.FALSE.) THEN
1990	!
1991	! CALL ros_Hermite5( chk_T(i), chk_T(i+1), T, &
1992	! chk_Y(1,i), chk_Y(1,i+1), &
1993	! chk_dY(1,i), chk_dY(1,i+1), &
1994	! chk_d2Y(1,i), chk_d2Y(1,i+1), Y )
1995	!
1996	! ELSE
1997
1998	CALL ros_Hermite3( chk_T(i), chk_T(i+1), T, &
1999	chk_Y(1:NVAR,i), chk_Y(1:NVAR,i+1), &
2000	chk_dY(1:NVAR,i), chk_dY(1:NVAR,i+1), &
2001	Y )
2002
2003	!
2004	! END IF
2005
2006	END SUBROUTINE ros_cadj_Y
2007
2008
2009	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2010	SUBROUTINE ros_Hermite3( a, b, T, Ya, Yb, Ja, Jb, Y )
2011	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2012	! Template for Hermite interpolation of order 5 on the interval [a,b]
2013	! P = c(1) + c(2)(x-a) + ... + c(4)(x-a)^3
2014	! P[a,b] = [Ya,Yb], P'[a,b] = [Ja,Jb]
2015	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2016	IMPLICIT NONE
2017	!~~~> Input variables
2018	KPP_REAL, INTENT(IN) :: a, b, T, Ya(NVAR), Yb(NVAR)
2019	KPP_REAL, INTENT(IN) :: Ja(NVAR), Jb(NVAR)
2020	!~~~> Output variables
2021	KPP_REAL, INTENT(OUT) :: Y(NVAR)
2022	!~~~> Local variables
2023	KPP_REAL :: Tau, amb(3), C(NVAR,4)
2024	KPP_REAL, PARAMETER :: ZERO = 0.0d0
2025	INTEGER :: i, j
2026
2027	amb(1) = 1.0d0/(a-b)
2028	DO i=2,3
2029	amb(i) = amb(i-1)*amb(1)
2030	END DO
2031
2032
2033	! c(1) = ya;
2034	CALL WCOPY(NVAR,Ya,1,C(1:NVAR,1),1)
2035	! c(2) = ja;
2036	CALL WCOPY(NVAR,Ja,1,C(1:NVAR,2),1)
2037	! c(3) = 2/(a-b)ja + 1/(a-b)jb - 3/(a - b)^2ya + 3/(a - b)^2yb ;
2038	CALL WCOPY(NVAR,Ya,1,C(1:NVAR,3),1)
2039	CALL WSCAL(NVAR,-3.0*amb(2),C(1:NVAR,3),1)
2040	CALL WAXPY(NVAR,3.0*amb(2),Yb,1,C(1:NVAR,3),1)
2041	CALL WAXPY(NVAR,2.0*amb(1),Ja,1,C(1:NVAR,3),1)
2042	CALL WAXPY(NVAR,amb(1),Jb,1,C(1:NVAR,3),1)
2043	! c(4) = 1/(a-b)^2ja + 1/(a-b)^2jb - 2/(a-b)^3ya + 2/(a-b)^3yb ;
2044	CALL WCOPY(NVAR,Ya,1,C(1:NVAR,4),1)
2045	CALL WSCAL(NVAR,-2.0*amb(3),C(1:NVAR,4),1)
2046	CALL WAXPY(NVAR,2.0*amb(3),Yb,1,C(1:NVAR,4),1)
2047	CALL WAXPY(NVAR,amb(2),Ja,1,C(1:NVAR,4),1)
2048	CALL WAXPY(NVAR,amb(2),Jb,1,C(1:NVAR,4),1)
2049
2050	Tau = T - a
2051	CALL WCOPY(NVAR,C(1:NVAR,4),1,Y,1)
2052	CALL WSCAL(NVAR,Tau**3,Y,1)
2053	DO j = 3,1,-1
2054	CALL WAXPY(NVAR,TAU**(j-1),C(1:NVAR,j),1,Y,1)
2055	END DO
2056
2057	END SUBROUTINE ros_Hermite3
2058
2059	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2060	SUBROUTINE ros_Hermite5( a, b, T, Ya, Yb, Ja, Jb, Ha, Hb, Y )
2061	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2062	! Template for Hermite interpolation of order 5 on the interval [a,b]
2063	! P = c(1) + c(2)(x-a) + ... + c(6)(x-a)^5
2064	! P[a,b] = [Ya,Yb], P'[a,b] = [Ja,Jb], P"[a,b] = [Ha,Hb]
2065	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2066	IMPLICIT NONE
2067	!~~~> Input variables
2068	KPP_REAL, INTENT(IN) :: a, b, T, Ya(NVAR), Yb(NVAR)
2069	KPP_REAL, INTENT(IN) :: Ja(NVAR), Jb(NVAR), Ha(NVAR), Hb(NVAR)
2070	!~~~> Output variables
2071	KPP_REAL, INTENT(OUT) :: Y(NVAR)
2072	!~~~> Local variables
2073	KPP_REAL :: Tau, amb(5), C(NVAR,6)
2074	KPP_REAL, PARAMETER :: ZERO = 0.0d0, HALF = 0.5d0
2075	INTEGER :: i, j
2076
2077	amb(1) = 1.0d0/(a-b)
2078	DO i=2,5
2079	amb(i) = amb(i-1)*amb(1)
2080	END DO
2081
2082	! c(1) = ya;
2083	CALL WCOPY(NVAR,Ya,1,C(1:NVAR,1),1)
2084	! c(2) = ja;
2085	CALL WCOPY(NVAR,Ja,1,C(1:NVAR,2),1)
2086	! c(3) = ha/2;
2087	CALL WCOPY(NVAR,Ha,1,C(1:NVAR,3),1)
2088	CALL WSCAL(NVAR,HALF,C(1:NVAR,3),1)
2089
2090	! c(4) = 10amb(3)ya - 10amb(3)yb - 6amb(2)ja - 4amb(2)jb + 1.5amb(1)ha - 0.5amb(1)hb ;
2091	CALL WCOPY(NVAR,Ya,1,C(1:NVAR,4),1)
2092	CALL WSCAL(NVAR,10.0*amb(3),C(1:NVAR,4),1)
2093	CALL WAXPY(NVAR,-10.0*amb(3),Yb,1,C(1:NVAR,4),1)
2094	CALL WAXPY(NVAR,-6.0*amb(2),Ja,1,C(1:NVAR,4),1)
2095	CALL WAXPY(NVAR,-4.0*amb(2),Jb,1,C(1:NVAR,4),1)
2096	CALL WAXPY(NVAR, 1.5*amb(1),Ha,1,C(1:NVAR,4),1)
2097	CALL WAXPY(NVAR,-0.5*amb(1),Hb,1,C(1:NVAR,4),1)
2098
2099	! c(5) = 15amb(4)ya - 15amb(4)yb - 8.amb(3)ja - 7amb(3)jb + 1.5amb(2)ha - 1amb(2)hb ;
2100	CALL WCOPY(NVAR,Ya,1,C(1:NVAR,5),1)
2101	CALL WSCAL(NVAR, 15.0*amb(4),C(1:NVAR,5),1)
2102	CALL WAXPY(NVAR,-15.0*amb(4),Yb,1,C(1:NVAR,5),1)
2103	CALL WAXPY(NVAR,-8.0*amb(3),Ja,1,C(1:NVAR,5),1)
2104	CALL WAXPY(NVAR,-7.0*amb(3),Jb,1,C(1:NVAR,5),1)
2105	CALL WAXPY(NVAR,1.5*amb(2),Ha,1,C(1:NVAR,5),1)
2106	CALL WAXPY(NVAR,-amb(2),Hb,1,C(1:NVAR,5),1)
2107
2108	! c(6) = 6amb(5)ya - 6amb(5)yb - 3.amb(4)ja - 3.amb(4)jb + 0.5amb(3)ha -0.5amb(3)hb ;
2109	CALL WCOPY(NVAR,Ya,1,C(1:NVAR,6),1)
2110	CALL WSCAL(NVAR, 6.0*amb(5),C(1:NVAR,6),1)
2111	CALL WAXPY(NVAR,-6.0*amb(5),Yb,1,C(1:NVAR,6),1)
2112	CALL WAXPY(NVAR,-3.0*amb(4),Ja,1,C(1:NVAR,6),1)
2113	CALL WAXPY(NVAR,-3.0*amb(4),Jb,1,C(1:NVAR,6),1)
2114	CALL WAXPY(NVAR, 0.5*amb(3),Ha,1,C(1:NVAR,6),1)
2115	CALL WAXPY(NVAR,-0.5*amb(3),Hb,1,C(1:NVAR,6),1)
2116
2117	Tau = T - a
2118	CALL WCOPY(NVAR,C(1:NVAR,6),1,Y,1)
2119	DO j = 5,1,-1
2120	CALL WSCAL(NVAR,Tau,Y,1)
2121	CALL WAXPY(NVAR,ONE,C(1:NVAR,j),1,Y,1)
2122	END DO
2123
2124	END SUBROUTINE ros_Hermite5
2125
2126	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2127	SUBROUTINE Ros2
2128	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2129	! --- AN L-STABLE METHOD, 2 stages, order 2
2130	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2131
2132	IMPLICIT NONE
2133	DOUBLE PRECISION g
2134
2135	g = 1.0d0 + 1.0d0/SQRT(2.0d0)
2136
2137	rosMethod = RS2
2138	!~~~> Name of the method
2139	ros_Name = 'ROS-2'
2140	!~~~> Number of stages
2141	ros_S = 2
2142
2143	!~~~> The coefficient matrices A and C are strictly lower triangular.
2144	! The lower triangular (subdiagonal) elements are stored in row-wise order:
2145	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
2146	! The general mapping formula is:
2147	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
2148	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
2149
2150	ros_A(1) = (1.d0)/g
2151	ros_C(1) = (-2.d0)/g
2152	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
2153	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
2154	ros_NewF(1) = .TRUE.
2155	ros_NewF(2) = .TRUE.
2156	!~~~> M_i = Coefficients for new step solution
2157	ros_M(1)= (3.d0)/(2.d0*g)
2158	ros_M(2)= (1.d0)/(2.d0*g)
2159	! E_i = Coefficients for error estimator
2160	ros_E(1) = 1.d0/(2.d0*g)
2161	ros_E(2) = 1.d0/(2.d0*g)
2162	!~~~> ros_ELO = estimator of local order - the minimum between the
2163	! main and the embedded scheme orders plus one
2164	ros_ELO = 2.0d0
2165	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
2166	ros_Alpha(1) = 0.0d0
2167	ros_Alpha(2) = 1.0d0
2168	!~~~> Gamma_i = \sum_j gamma_{i,j}
2169	ros_Gamma(1) = g
2170	ros_Gamma(2) =-g
2171
2172	END SUBROUTINE Ros2
2173
2174
2175	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2176	SUBROUTINE Ros3
2177	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2178	! --- AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations
2179	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2180
2181	IMPLICIT NONE
2182
2183	rosMethod = RS3
2184	!~~~> Name of the method
2185	ros_Name = 'ROS-3'
2186	!~~~> Number of stages
2187	ros_S = 3
2188
2189	!~~~> The coefficient matrices A and C are strictly lower triangular.
2190	! The lower triangular (subdiagonal) elements are stored in row-wise order:
2191	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
2192	! The general mapping formula is:
2193	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
2194	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
2195
2196	ros_A(1)= 1.d0
2197	ros_A(2)= 1.d0
2198	ros_A(3)= 0.d0
2199
2200	ros_C(1) = -0.10156171083877702091975600115545d+01
2201	ros_C(2) = 0.40759956452537699824805835358067d+01
2202	ros_C(3) = 0.92076794298330791242156818474003d+01
2203	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
2204	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
2205	ros_NewF(1) = .TRUE.
2206	ros_NewF(2) = .TRUE.
2207	ros_NewF(3) = .FALSE.
2208	!~~~> M_i = Coefficients for new step solution
2209	ros_M(1) = 0.1d+01
2210	ros_M(2) = 0.61697947043828245592553615689730d+01
2211	ros_M(3) = -0.42772256543218573326238373806514d+00
2212	! E_i = Coefficients for error estimator
2213	ros_E(1) = 0.5d+00
2214	ros_E(2) = -0.29079558716805469821718236208017d+01
2215	ros_E(3) = 0.22354069897811569627360909276199d+00
2216	!~~~> ros_ELO = estimator of local order - the minimum between the
2217	! main and the embedded scheme orders plus 1
2218	ros_ELO = 3.0d0
2219	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
2220	ros_Alpha(1)= 0.0d+00
2221	ros_Alpha(2)= 0.43586652150845899941601945119356d+00
2222	ros_Alpha(3)= 0.43586652150845899941601945119356d+00
2223	!~~~> Gamma_i = \sum_j gamma_{i,j}
2224	ros_Gamma(1)= 0.43586652150845899941601945119356d+00
2225	ros_Gamma(2)= 0.24291996454816804366592249683314d+00
2226	ros_Gamma(3)= 0.21851380027664058511513169485832d+01
2227
2228	END SUBROUTINE Ros3
2229
2230	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2231
2232
2233	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2234	SUBROUTINE Ros4
2235	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2236	! L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES
2237	! L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
2238	!
2239	! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
2240	! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
2241	! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
2242	! SPRINGER-VERLAG (1990)
2243	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2244
2245	IMPLICIT NONE
2246
2247	rosMethod = RS4
2248	!~~~> Name of the method
2249	ros_Name = 'ROS-4'
2250	!~~~> Number of stages
2251	ros_S = 4
2252
2253	!~~~> The coefficient matrices A and C are strictly lower triangular.
2254	! The lower triangular (subdiagonal) elements are stored in row-wise order:
2255	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
2256	! The general mapping formula is:
2257	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
2258	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
2259
2260	ros_A(1) = 0.2000000000000000d+01
2261	ros_A(2) = 0.1867943637803922d+01
2262	ros_A(3) = 0.2344449711399156d+00
2263	ros_A(4) = ros_A(2)
2264	ros_A(5) = ros_A(3)
2265	ros_A(6) = 0.0D0
2266
2267	ros_C(1) =-0.7137615036412310d+01
2268	ros_C(2) = 0.2580708087951457d+01
2269	ros_C(3) = 0.6515950076447975d+00
2270	ros_C(4) =-0.2137148994382534d+01
2271	ros_C(5) =-0.3214669691237626d+00
2272	ros_C(6) =-0.6949742501781779d+00
2273	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
2274	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
2275	ros_NewF(1) = .TRUE.
2276	ros_NewF(2) = .TRUE.
2277	ros_NewF(3) = .TRUE.
2278	ros_NewF(4) = .FALSE.
2279	!~~~> M_i = Coefficients for new step solution
2280	ros_M(1) = 0.2255570073418735d+01
2281	ros_M(2) = 0.2870493262186792d+00
2282	ros_M(3) = 0.4353179431840180d+00
2283	ros_M(4) = 0.1093502252409163d+01
2284	!~~~> E_i = Coefficients for error estimator
2285	ros_E(1) =-0.2815431932141155d+00
2286	ros_E(2) =-0.7276199124938920d-01
2287	ros_E(3) =-0.1082196201495311d+00
2288	ros_E(4) =-0.1093502252409163d+01
2289	!~~~> ros_ELO = estimator of local order - the minimum between the
2290	! main and the embedded scheme orders plus 1
2291	ros_ELO = 4.0d0
2292	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
2293	ros_Alpha(1) = 0.D0
2294	ros_Alpha(2) = 0.1145640000000000d+01
2295	ros_Alpha(3) = 0.6552168638155900d+00
2296	ros_Alpha(4) = ros_Alpha(3)
2297	!~~~> Gamma_i = \sum_j gamma_{i,j}
2298	ros_Gamma(1) = 0.5728200000000000d+00
2299	ros_Gamma(2) =-0.1769193891319233d+01
2300	ros_Gamma(3) = 0.7592633437920482d+00
2301	ros_Gamma(4) =-0.1049021087100450d+00
2302
2303	END SUBROUTINE Ros4
2304
2305	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2306	SUBROUTINE Rodas3
2307	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2308	! --- A STIFFLY-STABLE METHOD, 4 stages, order 3
2309	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2310
2311	IMPLICIT NONE
2312
2313	rosMethod = RD3
2314	!~~~> Name of the method
2315	ros_Name = 'RODAS-3'
2316	!~~~> Number of stages
2317	ros_S = 4
2318
2319	!~~~> The coefficient matrices A and C are strictly lower triangular.
2320	! The lower triangular (subdiagonal) elements are stored in row-wise order:
2321	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
2322	! The general mapping formula is:
2323	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
2324	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
2325
2326	ros_A(1) = 0.0d+00
2327	ros_A(2) = 2.0d+00
2328	ros_A(3) = 0.0d+00
2329	ros_A(4) = 2.0d+00
2330	ros_A(5) = 0.0d+00
2331	ros_A(6) = 1.0d+00
2332
2333	ros_C(1) = 4.0d+00
2334	ros_C(2) = 1.0d+00
2335	ros_C(3) =-1.0d+00
2336	ros_C(4) = 1.0d+00
2337	ros_C(5) =-1.0d+00
2338	ros_C(6) =-(8.0d+00/3.0d+00)
2339
2340	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
2341	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
2342	ros_NewF(1) = .TRUE.
2343	ros_NewF(2) = .FALSE.
2344	ros_NewF(3) = .TRUE.
2345	ros_NewF(4) = .TRUE.
2346	!~~~> M_i = Coefficients for new step solution
2347	ros_M(1) = 2.0d+00
2348	ros_M(2) = 0.0d+00
2349	ros_M(3) = 1.0d+00
2350	ros_M(4) = 1.0d+00
2351	!~~~> E_i = Coefficients for error estimator
2352	ros_E(1) = 0.0d+00
2353	ros_E(2) = 0.0d+00
2354	ros_E(3) = 0.0d+00
2355	ros_E(4) = 1.0d+00
2356	!~~~> ros_ELO = estimator of local order - the minimum between the
2357	! main and the embedded scheme orders plus 1
2358	ros_ELO = 3.0d+00
2359	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
2360	ros_Alpha(1) = 0.0d+00
2361	ros_Alpha(2) = 0.0d+00
2362	ros_Alpha(3) = 1.0d+00
2363	ros_Alpha(4) = 1.0d+00
2364	!~~~> Gamma_i = \sum_j gamma_{i,j}
2365	ros_Gamma(1) = 0.5d+00
2366	ros_Gamma(2) = 1.5d+00
2367	ros_Gamma(3) = 0.0d+00
2368	ros_Gamma(4) = 0.0d+00
2369
2370	END SUBROUTINE Rodas3
2371
2372	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2373	SUBROUTINE Rodas4
2374	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2375	! STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES
2376	!
2377	! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
2378	! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
2379	! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
2380	! SPRINGER-VERLAG (1996)
2381	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2382
2383	IMPLICIT NONE
2384
2385	rosMethod = RD4
2386	!~~~> Name of the method
2387	ros_Name = 'RODAS-4'
2388	!~~~> Number of stages
2389	ros_S = 6
2390
2391	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
2392	ros_Alpha(1) = 0.000d0
2393	ros_Alpha(2) = 0.386d0
2394	ros_Alpha(3) = 0.210d0
2395	ros_Alpha(4) = 0.630d0
2396	ros_Alpha(5) = 1.000d0
2397	ros_Alpha(6) = 1.000d0
2398
2399	!~~~> Gamma_i = \sum_j gamma_{i,j}
2400	ros_Gamma(1) = 0.2500000000000000d+00
2401	ros_Gamma(2) =-0.1043000000000000d+00
2402	ros_Gamma(3) = 0.1035000000000000d+00
2403	ros_Gamma(4) =-0.3620000000000023d-01
2404	ros_Gamma(5) = 0.0d0
2405	ros_Gamma(6) = 0.0d0
2406
2407	!~~~> The coefficient matrices A and C are strictly lower triangular.
2408	! The lower triangular (subdiagonal) elements are stored in row-wise order:
2409	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
2410	! The general mapping formula is: A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
2411	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
2412
2413	ros_A(1) = 0.1544000000000000d+01
2414	ros_A(2) = 0.9466785280815826d+00
2415	ros_A(3) = 0.2557011698983284d+00
2416	ros_A(4) = 0.3314825187068521d+01
2417	ros_A(5) = 0.2896124015972201d+01
2418	ros_A(6) = 0.9986419139977817d+00
2419	ros_A(7) = 0.1221224509226641d+01
2420	ros_A(8) = 0.6019134481288629d+01
2421	ros_A(9) = 0.1253708332932087d+02
2422	ros_A(10) =-0.6878860361058950d+00
2423	ros_A(11) = ros_A(7)
2424	ros_A(12) = ros_A(8)
2425	ros_A(13) = ros_A(9)
2426	ros_A(14) = ros_A(10)
2427	ros_A(15) = 1.0d+00
2428
2429	ros_C(1) =-0.5668800000000000d+01
2430	ros_C(2) =-0.2430093356833875d+01
2431	ros_C(3) =-0.2063599157091915d+00
2432	ros_C(4) =-0.1073529058151375d+00
2433	ros_C(5) =-0.9594562251023355d+01
2434	ros_C(6) =-0.2047028614809616d+02
2435	ros_C(7) = 0.7496443313967647d+01
2436	ros_C(8) =-0.1024680431464352d+02
2437	ros_C(9) =-0.3399990352819905d+02
2438	ros_C(10) = 0.1170890893206160d+02
2439	ros_C(11) = 0.8083246795921522d+01
2440	ros_C(12) =-0.7981132988064893d+01
2441	ros_C(13) =-0.3152159432874371d+02
2442	ros_C(14) = 0.1631930543123136d+02
2443	ros_C(15) =-0.6058818238834054d+01
2444
2445	!~~~> M_i = Coefficients for new step solution
2446	ros_M(1) = ros_A(7)
2447	ros_M(2) = ros_A(8)
2448	ros_M(3) = ros_A(9)
2449	ros_M(4) = ros_A(10)
2450	ros_M(5) = 1.0d+00
2451	ros_M(6) = 1.0d+00
2452
2453	!~~~> E_i = Coefficients for error estimator
2454	ros_E(1) = 0.0d+00
2455	ros_E(2) = 0.0d+00
2456	ros_E(3) = 0.0d+00
2457	ros_E(4) = 0.0d+00
2458	ros_E(5) = 0.0d+00
2459	ros_E(6) = 1.0d+00
2460
2461	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
2462	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
2463	ros_NewF(1) = .TRUE.
2464	ros_NewF(2) = .TRUE.
2465	ros_NewF(3) = .TRUE.
2466	ros_NewF(4) = .TRUE.
2467	ros_NewF(5) = .TRUE.
2468	ros_NewF(6) = .TRUE.
2469
2470	!~~~> ros_ELO = estimator of local order - the minimum between the
2471	! main and the embedded scheme orders plus 1
2472	ros_ELO = 4.0d0
2473
2474	END SUBROUTINE Rodas4
2475
2476
2477	END SUBROUTINE RosenbrockADJ ! and its internal procedures
2478
2479
2480	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2481	SUBROUTINE Fun_Template( T, Y, Ydot )
2482	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2483	! Template for the ODE function call.
2484	! Updates the rate coefficients (and possibly the fixed species) at each call
2485	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2486
2487	!~~~> Input variables
2488	KPP_REAL, INTENT(IN) :: T, Y(NVAR)
2489	!~~~> Output variables
2490	KPP_REAL, INTENT(OUT) :: Ydot(NVAR)
2491	!~~~> Local variables
2492	KPP_REAL :: Told
2493
2494	Told = TIME
2495	TIME = T
2496	CALL Update_SUN()
2497	CALL Update_RCONST()
2498	CALL Fun( Y, FIX, RCONST, Ydot )
2499	TIME = Told
2500
2501	END SUBROUTINE Fun_Template
2502
2503
2504	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2505	SUBROUTINE Jac_Template( T, Y, Jcb )
2506	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2507	! Template for the ODE Jacobian call.
2508	! Updates the rate coefficients (and possibly the fixed species) at each call
2509	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2510
2511	!~~~> Input variables
2512	KPP_REAL :: T, Y(NVAR)
2513	!~~~> Output variables
2514	#ifdef FULL_ALGEBRA
2515	KPP_REAL :: JV(LU_NONZERO), Jcb(NVAR,NVAR)
2516	#else
2517	KPP_REAL :: Jcb(LU_NONZERO)
2518	#endif
2519	!~~~> Local variables
2520	KPP_REAL :: Told
2521	#ifdef FULL_ALGEBRA
2522	INTEGER :: i, j
2523	#endif
2524
2525	Told = TIME
2526	TIME = T
2527	CALL Update_SUN()
2528	CALL Update_RCONST()
2529	#ifdef FULL_ALGEBRA
2530	CALL Jac_SP(Y, FIX, RCONST, JV)
2531	DO j=1,NVAR
2532	DO i=1,NVAR
2533	Jcb(i,j) = 0.0_dp
2534	END DO
2535	END DO
2536	DO i=1,LU_NONZERO
2537	Jcb(LU_IROW(i),LU_ICOL(i)) = JV(i)
2538	END DO
2539	#else
2540	CALL Jac_SP( Y, FIX, RCONST, Jcb )
2541	#endif
2542	TIME = Told
2543
2544	END SUBROUTINE Jac_Template
2545
2546
2547
2548	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2549	SUBROUTINE Hess_Template( T, Y, Hes )
2550	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2551	! Template for the ODE Hessian call.
2552	! Updates the rate coefficients (and possibly the fixed species) at each call
2553	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2554
2555	!~~~> Input variables
2556	KPP_REAL, INTENT(IN) :: T, Y(NVAR)
2557	!~~~> Output variables
2558	KPP_REAL, INTENT(OUT) :: Hes(NHESS)
2559	!~~~> Local variables
2560	KPP_REAL :: Told
2561
2562	Told = TIME
2563	TIME = T
2564	CALL Update_SUN()
2565	CALL Update_RCONST()
2566	CALL Hessian( Y, FIX, RCONST, Hes )
2567	TIME = Told
2568
2569	END SUBROUTINE Hess_Template
2570
2571	END MODULE KPP_ROOT_Integrator
2572
2573
2574
2575

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

Context Navigation

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

Download in other formats: