Home

Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

rosenbrock_soa.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: 70.0 KB

Line
1	MODULE KPP_ROOT_Integrator
2
3	USE KPP_ROOT_Precision
4	USE KPP_ROOT_Parameters
5	USE KPP_ROOT_Global
6	USE KPP_ROOT_Function
7	USE KPP_ROOT_Jacobian
8	USE KPP_ROOT_Hessian
9	USE KPP_ROOT_LinearAlgebra
10	USE KPP_ROOT_Rates
11
12	IMPLICIT NONE
13	PUBLIC
14	SAVE
15	!~~~> Statistics on the work performed by the Rosenbrock method
16	INTEGER :: Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng
17	INTEGER, PARAMETER :: ifun=11, ijac=12, istp=13, &
18	iacc=14, irej=15, idec=16, isol=17, &
19	isng=18, itexit=11,ihexit=12
20	!~~~> Checkpoints in memory
21	INTEGER, PARAMETER :: bufsize = 1500
22	INTEGER :: stack_ptr = 0 ! last written entry
23	KPP_REAL, DIMENSION(:), POINTER :: buf_H, buf_T
24	KPP_REAL, DIMENSION(:,:), POINTER :: buf_Y, buf_K
25	KPP_REAL, DIMENSION(:,:), POINTER :: buf_Y_tlm, buf_K_tlm
26
27	CONTAINS
28
29
30
31	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
32	SUBROUTINE INTEGRATE_SOA( NSOA, Y, Y_tlm, Y_adj, Y_soa, TIN, TOUT, &
33	AtolAdj, RtolAdj, ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U )
34	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
35
36	IMPLICIT NONE
37
38	!~~~> NSOA - No. of vectors to multiply SOA with
39	INTEGER :: NSOA
40	!~~~> Y - Forward model variables
41	KPP_REAL, INTENT(INOUT) :: Y(NVAR)
42	!~~~> Y_adj - Tangent linear variables
43	KPP_REAL, INTENT(INOUT) :: Y_tlm(NVAR,NSOA)
44	!~~~> Y_adj - First order adjoint
45	KPP_REAL, INTENT(INOUT) :: Y_adj(NVAR)
46	!~~~> Y_soa - Second order adjoint
47	KPP_REAL, INTENT(INOUT) :: Y_soa(NVAR,NSOA)
48	!~~~>
49	KPP_REAL, INTENT(IN) :: TIN ! TIN - Start Time
50	KPP_REAL, INTENT(IN) :: TOUT ! TOUT - End Time
51	!~~~> Optional input parameters and statistics
52	INTEGER, INTENT(IN), OPTIONAL :: ICNTRL_U(20)
53	KPP_REAL, INTENT(IN), OPTIONAL :: RCNTRL_U(20)
54	INTEGER, INTENT(OUT), OPTIONAL :: ISTATUS_U(20)
55	KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20)
56
57	INTEGER N_stp, N_acc, N_rej, N_sng, IERR
58	SAVE N_stp, N_acc, N_rej, N_sng
59	INTEGER i
60	KPP_REAL :: RCNTRL(20), RSTATUS(20)
61	KPP_REAL :: AtolAdj(NVAR), RtolAdj(NVAR)
62	INTEGER :: ICNTRL(20), ISTATUS(20)
63
64	ICNTRL(1:20) = 0
65	RCNTRL(1:20) = 0.0_dp
66	ISTATUS(1:20) = 0
67	RSTATUS(1:20) = 0.0_dp
68
69	ICNTRL(1) = 0 ! 0 = non-autonomous, 1 = autonomous
70	ICNTRL(2) = 1 ! 0 = scalar, 1 = vector tolerances
71	RCNTRL(3) = STEPMIN ! starting step
72	ICNTRL(4) = 5 ! choice of the method for forward and adjoint integration
73
74	! Tighter tolerances, especially atol, are needed for the full continuous adjoint
75	! (Atol on sensitivities is different than on concentrations)
76	! CADJ_ATOL(1:NVAR) = 1.0d-5
77	! CADJ_RTOL(1:NVAR) = 1.0d-4
78
79	! if optional parameters are given, and if they are >=0,
80	! then they overwrite default settings
81	IF (PRESENT(ICNTRL_U)) THEN
82	WHERE(ICNTRL_U(:) >= 0) ICNTRL(:) = ICNTRL_U(:)
83	ENDIF
84	IF (PRESENT(RCNTRL_U)) THEN
85	WHERE(RCNTRL_U(:) >= 0) RCNTRL(:) = RCNTRL_U(:)
86	ENDIF
87
88
89	CALL RosenbrockSOA(NSOA, &
90	Y, Y_tlm, Y_adj, Y_soa, &
91	TIN,TOUT, &
92	ATOL,RTOL, &
93	Fun_Template,Jac_Template,Hess_Template, &
94	RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR)
95
96
97	! N_stp = N_stp + ICNTRL(istp)
98	! N_acc = N_acc + ICNTRL(iacc)
99	! N_rej = N_rej + ICNTRL(irej)
100	! N_sng = N_sng + ICNTRL(isng)
101	! PRINT*,'Step=',N_stp,' Acc=',N_acc,' Rej=',N_rej, &
102	! ' Singular=',N_sng
103
104	IF (IERR < 0) THEN
105	print *,'RosenbrockSOA: Unsucessful step at T=', &
106	TIN,' (IERR=',IERR,')'
107	ENDIF
108
109	STEPMIN = RCNTRL(ihexit)
110	! if optional parameters are given for output they to return information
111	IF (PRESENT(ISTATUS_U)) ISTATUS_U(:) = ISTATUS(:)
112	IF (PRESENT(RSTATUS_U)) RSTATUS_U(:) = RSTATUS(:)
113
114	END SUBROUTINE INTEGRATE_SOA
115
116
117
118	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
119	SUBROUTINE RosenbrockSOA( NSOA, &
120	Y, Y_tlm, Y_adj, Y_soa, &
121	Tstart,Tend, &
122	AbsTol,RelTol, &
123	ode_Fun,ode_Jac , ode_Hess, &
124	RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR)
125	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
126	!
127	! ADJ = Adjoint of the Tangent Linear Model of a RosenbrockSOA Method
128	!
129	! Solves the system y'=F(t,y) using a RosenbrockSOA method defined by:
130	!
131	! G = 1/(H*gamma(1)) - ode_Jac(t0,Y0)
132	! T_i = t0 + Alpha(i)*H
133	! Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
134	! G * K_i = ode_Fun( T_i, Y_i ) + \sum_{j=1}^S C(i,j)/H * K_j +
135	! gamma(i)*dF/dT(t0, Y0)
136	! Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
137	!
138	! For details on RosenbrockSOA methods and their implementation consult:
139	! E. Hairer and G. Wanner
140	! "Solving ODEs II. Stiff and differential-algebraic problems".
141	! Springer series in computational mathematics, Springer-Verlag, 1996.
142	! The codes contained in the book inspired this implementation.
143	!
144	! (C) Adrian Sandu, August 2004
145	! Virginia Polytechnic Institute and State University
146	! Contact: sandu@cs.vt.edu
147	! This implementation is part of KPP - the Kinetic PreProcessor
148	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
149	!
150	!~~~> INPUT ARGUMENTS:
151	!
152	!- Y(NVAR) -> vector of initial conditions (at T=Tstart)
153	! NSOA -> dimension of linearized system,
154	! i.e. the number of sensitivity coefficients
155	!- Y_adj(NVAR) -> vector of initial sensitivity conditions (at T=Tstart)
156	!- [Tstart,Tend] = time range of integration
157	! (if Tstart>Tend the integration is performed backwards in time)
158	!- RelTol, AbsTol = user precribed accuracy
159	!- SUBROUTINE ode_Fun( T, Y, Ydot ) = ODE function,
160	! returns Ydot = Y' = F(T,Y)
161	!- SUBROUTINE ode_Fun( T, Y, Ydot ) = Jacobian of the ODE function,
162	! returns Jcb = dF/dY
163	!- ICNTRL(1:10) = integer inputs parameters
164	!- RCNTRL(1:10) = real inputs parameters
165	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
166	!
167	!~~~> OUTPUT ARGUMENTS:
168	!
169	!- Y(NVAR) -> vector of final states (at T->Tend)
170	!- Y_adj(NVAR) -> vector of final sensitivities (at T=Tend)
171	!- ICNTRL(11:20) -> integer output parameters
172	!- RCNTRL(11:20) -> real output parameters
173	!- IERR -> job status upon return
174	! - succes (positive value) or failure (negative value) -
175	! = 1 : Success
176	! = -1 : Improper value for maximal no of steps
177	! = -2 : Selected RosenbrockSOA method not implemented
178	! = -3 : Hmin/Hmax/Hstart must be positive
179	! = -4 : FacMin/FacMax/FacRej must be positive
180	! = -5 : Improper tolerance values
181	! = -6 : No of steps exceeds maximum bound
182	! = -7 : Step size too small
183	! = -8 : Matrix is repeatedly singular
184	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
185	!
186	!~~~> INPUT PARAMETERS:
187	!
188	! Note: For input parameters equal to zero the default values of the
189	! corresponding variables are used.
190	!
191	! ICNTRL(1) = 1: F = F(y) Independent of T (AUTONOMOUS)
192	! = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
193	! ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors
194	! = 1: AbsTol, RelTol are scalars
195	! ICNTRL(3) -> maximum number of integration steps
196	! For ICNTRL(3)=0) the default value of 100000 is used
197	!
198	! ICNTRL(4) -> selection of a particular Rosenbrock method
199	! = 0 : default method is Rodas3
200	! = 1 : method is Ros2
201	! = 2 : method is Ros3
202	! = 3 : method is Ros4
203	! = 4 : method is Rodas3
204	! = 5: method is Rodas4
205	!
206	! ICNTRL(5) -> Type of adjoint algorithm
207	! = 0 : default is discrete adjoint ( of method ICNTRL(4) )
208	! = 1 : no adjoint
209	! = 2 : discrete adjoint ( of method ICNTRL(4) )
210	! = 3 : fully adaptive continuous adjoint ( with method ICNTRL(6) )
211	! = 4 : simplified continuous adjoint ( with method ICNTRL(6) )
212	!
213	! ICNTRL(6) -> selection of a particular Rosenbrock method for the
214	! continuous adjoint integration - for cts adjoint it
215	! can be different than the forward method ICNTRL(4)
216	! Note 1: to avoid interpolation errors (which can be huge!)
217	! it is recommended to use only ICNTRL(6) = 1 or 4
218	! Note 2: the performance of the full continuous adjoint
219	! strongly depends on the forward solution accuracy Abs/RelTol
220	!
221	! RCNTRL(1) -> Hmin, lower bound for the integration step size
222	! It is strongly recommended to keep Hmin = ZERO
223	! RCNTRL(2) -> Hmax, upper bound for the integration step size
224	! RCNTRL(3) -> Hstart, starting value for the integration step size
225	!
226	! RCNTRL(4) -> FacMin, lower bound on step decrease factor (default=0.2)
227	! RCNTRL(5) -> FacMin,upper bound on step increase factor (default=6)
228	! RCNTRL(6) -> FacRej, step decrease factor after multiple rejections
229	! (default=0.1)
230	! RCNTRL(7) -> FacSafe, by which the new step is slightly smaller
231	! than the predicted value (default=0.9)
232	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
233	!
234	!~~~> OUTPUT PARAMETERS:
235	!
236	! Note: each call to RosenbrockSOA adds the corrent no. of fcn calls
237	! to previous value of ISTATUS(1), and similar for the other params.
238	! Set ISTATUS(1:10) = 0 before call to avoid this accumulation.
239	!
240	! ISTATUS(1) = No. of function calls
241	! ISTATUS(2) = No. of jacobian calls
242	! ISTATUS(3) = No. of steps
243	! ISTATUS(4) = No. of accepted steps
244	! ISTATUS(5) = No. of rejected steps (except at the beginning)
245	! ISTATUS(6) = No. of LU decompositions
246	! ISTATUS(7) = No. of forward/backward substitutions
247	! ISTATUS(8) = No. of singular matrix decompositions
248	!
249	! RSTATUS(1) -> Texit, the time corresponding to the
250	! computed Y upon return
251	! RSTATUS(2) -> Hexit, last accepted step before exit
252	! For multiple restarts, use Hexit as Hstart in the following run
253	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
254
255	IMPLICIT NONE
256
257	!~~~> Arguments
258	INTEGER, INTENT(IN) :: NSOA
259	KPP_REAL, INTENT(INOUT) :: Y(NVAR)
260	KPP_REAL, INTENT(INOUT) :: Y_tlm(NVAR,NSOA)
261	KPP_REAL, INTENT(INOUT) :: Y_adj(NVAR)
262	KPP_REAL, INTENT(INOUT) :: Y_soa(NVAR,NSOA)
263	KPP_REAL, INTENT(IN) :: Tstart, Tend
264	KPP_REAL, INTENT(IN) :: AbsTol(NVAR),RelTol(NVAR)
265	INTEGER, INTENT(IN) :: ICNTRL(20)
266	KPP_REAL, INTENT(IN) :: RCNTRL(20)
267	INTEGER, INTENT(INOUT) :: ISTATUS(20)
268	KPP_REAL, INTENT(INOUT) :: RSTATUS(20)
269	INTEGER, INTENT(OUT) :: IERR
270	!~~~> The method parameters
271	INTEGER, PARAMETER :: Smax = 6
272	INTEGER :: Method, ros_S
273	KPP_REAL, DIMENSION(Smax) :: ros_M, ros_E, ros_Alpha, ros_Gamma
274	KPP_REAL, DIMENSION(Smax*(Smax-1)/2) :: ros_A, ros_C
275	KPP_REAL :: ros_ELO
276	LOGICAL, DIMENSION(Smax) :: ros_NewF
277	CHARACTER(LEN=12) :: ros_Name
278	!~~~> Local variables
279	KPP_REAL :: Roundoff, FacMin, FacMax, FacRej, FacSafe
280	KPP_REAL :: Hmin, Hmax, Hstart, Hexit
281	KPP_REAL :: T
282	INTEGER :: i, UplimTol, Max_no_steps
283	LOGICAL :: Autonomous, VectorTol
284	!~~~> Parameters
285	KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
286	KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5
287	!~~~> Functions
288	EXTERNAL ode_Fun, ode_Jac, ode_Hess
289
290	!~~~> Initialize statistics
291	Nfun = ISTATUS(ifun)
292	Njac = ISTATUS(ijac)
293	Nstp = ISTATUS(istp)
294	Nacc = ISTATUS(iacc)
295	Nrej = ISTATUS(irej)
296	Ndec = ISTATUS(idec)
297	Nsol = ISTATUS(isol)
298	Nsng = ISTATUS(isng)
299
300	!~~~> Autonomous or time dependent ODE. Default is time dependent.
301	Autonomous = .NOT.(ICNTRL(1) == 0)
302
303	!~~~> For Scalar tolerances (ICNTRL(2).NE.0)
304	! the code uses AbsTol(1) and RelTol(1)
305	! For Vector tolerances (ICNTRL(2) == 0)
306	! the code uses AbsTol(1:NVAR) and RelTol(1:NVAR)
307	IF (ICNTRL(2) == 0) THEN
308	VectorTol = .TRUE.
309	UplimTol = NVAR
310	ELSE
311	VectorTol = .FALSE.
312	UplimTol = 1
313	END IF
314
315	!~~~> The maximum number of steps admitted
316	IF (ICNTRL(3) == 0) THEN
317	Max_no_steps = bufsize - 1
318	ELSEIF (Max_no_steps > 0) THEN
319	Max_no_steps=ICNTRL(3)
320	ELSE
321	PRINT * ,'User-selected max no. of steps: ICNTRL(3)=',ICNTRL(3)
322	CALL ros_ErrorMsg(-1,Tstart,ZERO,IERR)
323	RETURN
324	END IF
325
326	!~~~> The particular Rosenbrock method chosen
327	IF (ICNTRL(4) == 0) THEN
328	Method = 5
329	ELSEIF ( (ICNTRL(4) >= 1).AND.(ICNTRL(4) <= 5) ) THEN
330	Method = ICNTRL(4)
331	ELSE
332	PRINT * , 'User-selected Rosenbrock method: ICNTRL(4)=', Method
333	CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
334	RETURN
335	END IF
336
337
338	!~~~> Unit roundoff (1+Roundoff>1)
339	Roundoff = WLAMCH('E')
340
341	!~~~> Lower bound on the step size: (positive value)
342	IF (RCNTRL(1) == ZERO) THEN
343	Hmin = ZERO
344	ELSEIF (RCNTRL(1) > ZERO) THEN
345	Hmin = RCNTRL(1)
346	ELSE
347	PRINT * , 'User-selected Hmin: RCNTRL(1)=', RCNTRL(1)
348	CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
349	RETURN
350	END IF
351	!~~~> Upper bound on the step size: (positive value)
352	IF (RCNTRL(2) == ZERO) THEN
353	Hmax = ABS(Tend-Tstart)
354	ELSEIF (RCNTRL(2) > ZERO) THEN
355	Hmax = MIN(ABS(RCNTRL(2)),ABS(Tend-Tstart))
356	ELSE
357	PRINT * , 'User-selected Hmax: RCNTRL(2)=', RCNTRL(2)
358	CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
359	RETURN
360	END IF
361	!~~~> Starting step size: (positive value)
362	IF (RCNTRL(3) == ZERO) THEN
363	Hstart = MAX(Hmin,DeltaMin)
364	ELSEIF (RCNTRL(3) > ZERO) THEN
365	Hstart = MIN(ABS(RCNTRL(3)),ABS(Tend-Tstart))
366	ELSE
367	PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3)
368	CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
369	RETURN
370	END IF
371	!~~~> Step size can be changed s.t. FacMin < Hnew/Hexit < FacMax
372	IF (RCNTRL(4) == ZERO) THEN
373	FacMin = 0.2d0
374	ELSEIF (RCNTRL(4) > ZERO) THEN
375	FacMin = RCNTRL(4)
376	ELSE
377	PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4)
378	CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
379	RETURN
380	END IF
381	IF (RCNTRL(5) == ZERO) THEN
382	FacMax = 6.0d0
383	ELSEIF (RCNTRL(5) > ZERO) THEN
384	FacMax = RCNTRL(5)
385	ELSE
386	PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5)
387	CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
388	RETURN
389	END IF
390	!~~~> FacRej: Factor to decrease step after 2 succesive rejections
391	IF (RCNTRL(6) == ZERO) THEN
392	FacRej = 0.1d0
393	ELSEIF (RCNTRL(6) > ZERO) THEN
394	FacRej = RCNTRL(6)
395	ELSE
396	PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6)
397	CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
398	RETURN
399	END IF
400	!~~~> FacSafe: Safety Factor in the computation of new step size
401	IF (RCNTRL(7) == ZERO) THEN
402	FacSafe = 0.9d0
403	ELSEIF (RCNTRL(7) > ZERO) THEN
404	FacSafe = RCNTRL(7)
405	ELSE
406	PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7)
407	CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
408	RETURN
409	END IF
410	!~~~> Check if tolerances are reasonable
411	DO i=1,UplimTol
412	IF ( (AbsTol(i) <= ZERO) .OR. (RelTol(i) <= 10.d0*Roundoff) &
413	.OR. (RelTol(i) >= 1.0d0) ) THEN
414	PRINT * , ' AbsTol(',i,') = ',AbsTol(i)
415	PRINT * , ' RelTol(',i,') = ',RelTol(i)
416	CALL ros_ErrorMsg(-5,Tstart,ZERO,IERR)
417	RETURN
418	END IF
419	END DO
420
421
422	!~~~> Initialize the particular RosenbrockSOA method
423	SELECT CASE (Method)
424	CASE (1)
425	CALL Ros2(ros_S, ros_A, ros_C, ros_M, ros_E, &
426	ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
427	CASE (2)
428	CALL Ros3(ros_S, ros_A, ros_C, ros_M, ros_E, &
429	ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
430	CASE (3)
431	CALL Ros4(ros_S, ros_A, ros_C, ros_M, ros_E, &
432	ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
433	CASE (4)
434	CALL Rodas3(ros_S, ros_A, ros_C, ros_M, ros_E, &
435	ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
436	CASE (5)
437	CALL Rodas4(ros_S, ros_A, ros_C, ros_M, ros_E, &
438	ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
439	CASE DEFAULT
440	PRINT * , 'Unknown Rosenbrock method: ICNTRL(4)=', Method
441	CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
442	RETURN
443	END SELECT
444
445	!~~~> Allocate checkpoint space or open checkpoint files
446	CALL ros_AllocateDBuffers( ros_S )
447
448	!~~~> Forward Rosenbrock and TLM integration
449	CALL ros_TlmInt (NSOA, Y, Y_tlm, &
450	Tstart, Tend, T, &
451	AbsTol, RelTol, &
452	ode_Fun, ode_Jac, ode_Hess, &
453	!~~~> Rosenbrock method coefficients
454	ros_S, ros_M, ros_E, ros_A, ros_C, &
455	ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, &
456	!~~~> Integration parameters
457	Autonomous, VectorTol, Max_no_steps, &
458	Roundoff, Hmin, Hmax, Hstart, Hexit, &
459	FacMin, FacMax, FacRej, FacSafe, &
460	!~~~> Error indicator
461	IERR )
462
463	PRINT*,'FORWARD STATISTICS'
464	PRINT*,'Step=',Nstp,' Acc=',Nacc, &
465	' Rej=',Nrej, ' Singular=',Nsng
466	Nstp = 0
467	Nacc = 0
468	Nrej = 0
469	Nsng = 0
470
471	!~~~> If Forward integration failed return
472	IF (IERR<0) RETURN
473
474	!~~~> Backward ADJ and SOADJ Rosenbrock integration
475	CALL ros_SoaInt ( &
476	NSOA, Y_adj, Y_soa, &
477	Tstart, Tend, T, &
478	AbsTol, RelTol, &
479	ode_Fun, ode_Jac, ode_Hess, &
480	!~~~> RosenbrockSOA method coefficients
481	ros_S, ros_M, ros_E, ros_A, ros_C, &
482	ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, &
483	!~~~> Integration parameters
484	Autonomous, VectorTol, Max_no_steps, &
485	Roundoff, Hmin, Hmax, Hstart, &
486	FacMin, FacMax, FacRej, FacSafe, &
487	!~~~> Error indicator
488	IERR )
489
490
491	PRINT*,'ADJOINT STATISTICS'
492	PRINT*,'Step=',Nstp,' Acc=',Nacc, &
493	' Rej=',Nrej, ' Singular=',Nsng
494
495	!~~~> Free checkpoint space or close checkpoint files
496	CALL ros_FreeDBuffers
497
498	!~~~> Collect run statistics
499	ISTATUS(ifun) = Nfun
500	ISTATUS(ijac) = Njac
501	ISTATUS(istp) = Nstp
502	ISTATUS(iacc) = Nacc
503	ISTATUS(irej) = Nrej
504	ISTATUS(idec) = Ndec
505	ISTATUS(isol) = Nsol
506	ISTATUS(isng) = Nsng
507	!~~~> Last T and H
508	RSTATUS(itexit) = T
509	RSTATUS(ihexit) = Hexit
510
511	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
512	END SUBROUTINE RosenbrockSOA
513	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
514
515
516	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
517	SUBROUTINE ros_ErrorMsg(Code,T,H,IERR)
518	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
519	! Handles all error messages
520	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
521
522	KPP_REAL, INTENT(IN) :: T, H
523	INTEGER, INTENT(IN) :: Code
524	INTEGER, INTENT(OUT) :: IERR
525
526	IERR = Code
527	PRINT * , &
528	'Forced exit from RosenbrockSOA due to the following error:'
529
530	SELECT CASE (Code)
531	CASE (-1)
532	PRINT * , '--> Improper value for maximal no of steps'
533	CASE (-2)
534	PRINT * , '--> Selected RosenbrockSOA method not implemented'
535	CASE (-3)
536	PRINT * , '--> Hmin/Hmax/Hstart must be positive'
537	CASE (-4)
538	PRINT * , '--> FacMin/FacMax/FacRej must be positive'
539	CASE (-5)
540	PRINT * , '--> Improper tolerance values'
541	CASE (-6)
542	PRINT * , '--> No of steps exceeds maximum bound'
543	CASE (-7)
544	PRINT * , '--> Step size too small: T + 10*H = T', &
545	' or H < Roundoff'
546	CASE (-8)
547	PRINT * , '--> Matrix is repeatedly singular'
548	CASE (-9)
549	PRINT * , '--> Improper type of adjoint selected'
550	CASE DEFAULT
551	PRINT *, 'Unknown Error code: ', Code
552	END SELECT
553
554	PRINT *, "T=", T, "and H=", H
555
556	END SUBROUTINE ros_ErrorMsg
557
558
559
560
561	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
562	SUBROUTINE ros_TlmInt (NSOA, Y, Y_tlm, &
563	Tstart, Tend, T, &
564	AbsTol, RelTol, &
565	ode_Fun, ode_Jac, ode_Hess, &
566	!~~~> Rosenbrock method coefficients
567	ros_S, ros_M, ros_E, ros_A, ros_C, &
568	ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, &
569	!~~~> Integration parameters
570	Autonomous, VectorTol, Max_no_steps, &
571	Roundoff, Hmin, Hmax, Hstart, Hexit, &
572	FacMin, FacMax, FacRej, FacSafe, &
573	!~~~> Error indicator
574	IERR )
575	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
576	! Template for the implementation of a generic Rosenbrock method
577	! defined by ros_S (no of stages)
578	! and its coefficients ros_{A,C,M,E,Alpha,Gamma}
579	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
580
581	IMPLICIT NONE
582
583	!~~~> Input: the initial condition at Tstart; Output: the solution at T
584	KPP_REAL, INTENT(INOUT) :: Y(NVAR)
585	!~~~> Input: Number of sensitivity coefficients
586	INTEGER, INTENT(IN) :: NSOA
587	!~~~> Input: the initial sensitivites at Tstart; Output: the sensitivities at T
588	KPP_REAL, INTENT(INOUT) :: Y_tlm(NVAR,NSOA)
589	!~~~> Input: integration interval
590	KPP_REAL, INTENT(IN) :: Tstart,Tend
591	!~~~> Output: time at which the solution is returned (T=Tend if success)
592	KPP_REAL, INTENT(OUT) :: T
593	!~~~> Input: tolerances
594	KPP_REAL, INTENT(IN) :: AbsTol(NVAR), RelTol(NVAR)
595	!~~~> Input: ode function and its Jacobian
596	EXTERNAL ode_Fun, ode_Jac, ode_Hess
597	!~~~> Input: The Rosenbrock method parameters
598	INTEGER, INTENT(IN) :: ros_S
599	KPP_REAL, INTENT(IN) :: ros_M(ros_S), ros_E(ros_S), &
600	ros_Alpha(ros_S), ros_A(ros_S*(ros_S-1)/2), &
601	ros_Gamma(ros_S), ros_C(ros_S*(ros_S-1)/2), ros_ELO
602	LOGICAL, INTENT(IN) :: ros_NewF(ros_S)
603	!~~~> Input: integration parameters
604	LOGICAL, INTENT(IN) :: Autonomous, VectorTol
605	KPP_REAL, INTENT(IN) :: Hstart, Hmin, Hmax
606	INTEGER, INTENT(IN) :: Max_no_steps
607	KPP_REAL, INTENT(IN) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
608	!~~~> Output: last accepted step
609	KPP_REAL, INTENT(OUT) :: Hexit
610	!~~~> Output: Error indicator
611	INTEGER, INTENT(OUT) :: IERR
612	! ~~~~ Local variables
613	KPP_REAL :: Ystage(NVAR*ros_S), Fcn0(NVAR), Fcn(NVAR)
614	KPP_REAL :: K(NVAR*ros_S), Tmp(NVAR)
615	KPP_REAL :: Ystage_tlm(NVAR*ros_S,NSOA), Fcn0_tlm(NVAR,NSOA), Fcn_tlm(NVAR,NSOA)
616	KPP_REAL :: K_tlm(NVAR*ros_S,NSOA)
617	KPP_REAL :: Hes0(NHESS)
618	KPP_REAL :: dFdT(NVAR), dJdT(LU_NONZERO)
619	KPP_REAL :: Jac0(LU_NONZERO), Jac(LU_NONZERO), Ghimj(LU_NONZERO)
620	KPP_REAL :: H, Hnew, HC, HG, Fac, Tau
621	KPP_REAL :: Err, Yerr(NVAR), Ynew(NVAR), Ynew_tlm(NVAR,NSOA)
622	INTEGER :: Pivot(NVAR), Direction, ioffset, joffset, j, istage, mtlm
623	LOGICAL :: RejectLastH, RejectMoreH, Singular
624	!~~~> Local parameters
625	KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
626	KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5
627	!~~~> Locally called functions
628	! KPP_REAL WLAMCH
629	! EXTERNAL WLAMCH
630	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
631
632
633	!~~~> Initial preparations
634	T = Tstart
635	Hexit = 0.0_dp
636	H = MIN(Hstart,Hmax)
637	IF (ABS(H) <= 10.D0*Roundoff) H = DeltaMin
638
639	IF (Tend >= Tstart) THEN
640	Direction = +1
641	ELSE
642	Direction = -1
643	END IF
644
645	RejectLastH=.FALSE.
646	RejectMoreH=.FALSE.
647
648	!~~~> Time loop begins below
649
650	TimeLoop: DO WHILE ( (Direction > 0).AND.((T-Tend)+Roundoff <= ZERO) &
651	.OR. (Direction < 0).AND.((Tend-T)+Roundoff <= ZERO) )
652
653	IF ( Nstp > Max_no_steps ) THEN ! Too many steps
654	CALL ros_ErrorMsg(-6,T,H,IERR)
655	RETURN
656	END IF
657	IF ( ((T+0.1d0*H) == T).OR.(H <= Roundoff) ) THEN ! Step size too small
658	CALL ros_ErrorMsg(-7,T,H,IERR)
659	RETURN
660	END IF
661
662	!~~~> Limit H if necessary to avoid going beyond Tend
663	Hexit = H
664	H = MIN(H,ABS(Tend-T))
665
666	!~~~> Compute the function at current time
667	CALL ode_Fun(T,Y,Fcn0)
668
669	!~~~> Compute the Jacobian at current time
670	CALL ode_Jac(T,Y,Jac0)
671
672	!~~~> Compute the Hessian at current time
673	CALL ode_Hess(T,Y,Hes0)
674
675	!~~~> Compute the TLM function at current time
676	DO mtlm = 1, NSOA
677	CALL Jac_SP_Vec ( Jac0, Y_tlm(1:NVAR,mtlm), Fcn0_tlm(1:NVAR,mtlm) )
678	END DO
679
680	!~~~> Compute the function and Jacobian derivatives with respect to T
681	IF (.NOT.Autonomous) THEN
682	CALL ros_FunTimeDerivative ( T, Roundoff, Y, &
683	Fcn0, ode_Fun, dFdT )
684	CALL ros_JacTimeDerivative ( T, Roundoff, Y, &
685	Jac0, ode_Jac, dJdT )
686	END IF
687
688	!~~~> Repeat step calculation until current step accepted
689	UntilAccepted: DO
690
691	CALL ros_PrepareMatrix(H,Direction,ros_Gamma(1),&
692	Jac0,Ghimj,Pivot,Singular)
693	IF (Singular) THEN ! More than 5 consecutive failed decompositions
694	CALL ros_ErrorMsg(-8,T,H,IERR)
695	RETURN
696	END IF
697
698	!~~~> Compute the stages
699	Stage: DO istage = 1, ros_S
700
701	! Current istage vector is K(ioffset+1:ioffset+NVAR:ioffset+1+NVAR-1)
702	ioffset = NVAR*(istage-1)
703
704	! For the 1st istage the function has been computed previously
705	IF ( istage == 1 ) THEN
706	CALL WCOPY(NVAR,Y,1,Ystage(ioffset+1:ioffset+NVAR),1)
707	DO mtlm=1,NSOA
708	CALL WCOPY(NVAR,Y_tlm(1:NVAR,mtlm),1,Ystage_tlm(ioffset+1:ioffset+NVAR,mtlm),1)
709	END DO
710	CALL WCOPY(NVAR,Fcn0,1,Fcn,1)
711	CALL WCOPY(NVAR*NSOA,Fcn0_tlm,1,Fcn_tlm,1)
712	! istage>1 and a new function evaluation is needed at the current istage
713	ELSEIF ( ros_NewF(istage) ) THEN
714	CALL WCOPY(NVAR,Y,1,Ystage(ioffset+1:ioffset+NVAR),1)
715	DO mtlm=1,NSOA
716	CALL WCOPY(NVAR,Y_tlm(1:NVAR,mtlm),1,Ystage_tlm(ioffset+1:ioffset+NVAR,mtlm),1)
717	END DO
718	DO j = 1, istage-1
719	joffset = NVAR*(j-1)
720	CALL WAXPY(NVAR,ros_A((istage-1)*(istage-2)/2+j), &
721	K(joffset+1:joffset+NVAR),1,Ystage(ioffset+1:ioffset+NVAR),1)
722	DO mtlm=1,NSOA
723	CALL WAXPY(NVAR,ros_A((istage-1)*(istage-2)/2+j), &
724	K_tlm(joffset+1:joffset+NVAR,mtlm),1,Ystage_tlm(ioffset+1:ioffset+NVAR,mtlm),1)
725	END DO
726	END DO
727	Tau = T + ros_Alpha(istage)DirectionH
728	CALL ode_Fun(Tau,Ystage(ioffset+1:ioffset+NVAR),Fcn)
729	CALL ode_Jac(Tau,Ystage(ioffset+1:ioffset+NVAR),Jac)
730	DO mtlm=1,NSOA
731	CALL Jac_SP_Vec ( Jac, Ystage_tlm(ioffset+1:ioffset+NVAR,mtlm), Fcn_tlm(1:NVAR,mtlm) )
732	END DO
733	END IF ! if istage == 1 elseif ros_NewF(istage)
734
735	CALL WCOPY(NVAR,Fcn,1,K(ioffset+1:ioffset+NVAR),1)
736	DO mtlm=1,NSOA
737	CALL WCOPY(NVAR,Fcn_tlm(1:NVAR,mtlm),1,K_tlm(ioffset+1:ioffset+NVAR,mtlm),1)
738	END DO
739	DO j = 1, istage-1
740	HC = ros_C((istage-1)(istage-2)/2+j)/(DirectionH)
741	CALL WAXPY(NVAR,HC,K(NVAR(j-1)+1:NVARj),1,K(ioffset+1:ioffset+NVAR),1)
742	DO mtlm=1,NSOA
743	CALL WAXPY(NVAR,HC,K_tlm(NVAR(j-1)+1:NVARj,mtlm),1,K_tlm(ioffset+1:ioffset+NVAR,mtlm),1)
744	END DO
745	END DO
746	IF ((.NOT. Autonomous).AND.(ros_Gamma(istage).NE.ZERO)) THEN
747	HG = DirectionHros_Gamma(istage)
748	CALL WAXPY(NVAR,HG,dFdT,1,K(ioffset+1:ioffset+NVAR),1)
749	DO mtlm=1,NSOA
750	CALL Jac_SP_Vec ( dJdT, Ystage_tlm(ioffset+1:ioffset+NVAR,mtlm), Fcn_tlm(1:NVAR,mtlm) )
751	CALL WAXPY(NVAR,HG,Fcn_tlm(1:NVAR,mtlm),1,K_tlm(ioffset+1:ioffset+NVAR,mtlm),1)
752	END DO
753	END IF
754	CALL ros_Solve('N', Ghimj, Pivot, K(ioffset+1:ioffset+NVAR))
755	DO mtlm=1,NSOA
756	CALL Hess_Vec ( Hes0, K(ioffset+1:ioffset+NVAR), Y_tlm(1:NVAR,mtlm), Tmp )
757	CALL WAXPY(NVAR,ONE,Tmp,1,K_tlm(ioffset+1:ioffset+NVAR,mtlm),1)
758	CALL ros_Solve('N', Ghimj, Pivot, K_tlm(ioffset+1:ioffset+NVAR,mtlm))
759	END DO
760
761	END DO Stage
762
763
764	!~~~> Compute the new solution
765	CALL WCOPY(NVAR,Y,1,Ynew,1)
766	DO j=1,ros_S
767	CALL WAXPY(NVAR,ros_M(j),K(NVAR(j-1)+1:NVARj),1,Ynew,1)
768	END DO
769	DO mtlm=1,NSOA
770	CALL WCOPY(NVAR,Y_tlm(1:NVAR,mtlm),1,Ynew_tlm(1:NVAR,mtlm),1)
771	DO j=1,ros_S
772	joffset = NVAR*(j-1)
773	CALL WAXPY(NVAR,ros_M(j),K_tlm(joffset+1:joffset+NVAR,mtlm),1,Ynew_tlm(1:NVAR,mtlm),1)
774	END DO
775	END DO
776
777	!~~~> Compute the error estimation
778	CALL WSCAL(NVAR,ZERO,Yerr,1)
779	DO j=1,ros_S
780	CALL WAXPY(NVAR,ros_E(j),K(NVAR(j-1)+1:NVARj),1,Yerr,1)
781	END DO
782	Err = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol )
783
784	!~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax
785	Fac = MIN(FacMax,MAX(FacMin,FacSafe/Err**(ONE/ros_ELO)))
786	Hnew = H*Fac
787
788	!~~~> Check the error magnitude and adjust step size
789	Nstp = Nstp+1
790	IF ( (Err <= ONE).OR.(H <= Hmin) ) THEN !~~~> Accept step
791	Nacc = Nacc+1
792	!~~~> Checkpoints for stage values and vectors
793	CALL ros_DPush( ros_S, NSOA, T, H, Ystage, K, Ystage_tlm, K_tlm )
794	!~~~> Accept new solution, etc.
795	CALL WCOPY(NVAR,Ynew,1,Y,1)
796	CALL WCOPY(NVAR*NSOA,Ynew_tlm,1,Y_tlm,1)
797	T = T + Direction*H
798	Hnew = MAX(Hmin,MIN(Hnew,Hmax))
799	IF (RejectLastH) THEN ! No step size increase after a rejected step
800	Hnew = MIN(Hnew,H)
801	END IF
802	RejectLastH = .FALSE.
803	RejectMoreH = .FALSE.
804	H = Hnew
805	EXIT UntilAccepted ! EXIT THE LOOP: WHILE STEP NOT ACCEPTED
806	ELSE !~~~> Reject step
807	IF (RejectMoreH) THEN
808	Hnew = H*FacRej
809	END IF
810	RejectMoreH = RejectLastH
811	RejectLastH = .TRUE.
812	H = Hnew
813	IF (Nacc >= 1) THEN
814	Nrej = Nrej+1
815	END IF
816	END IF ! Err <= 1
817
818	END DO UntilAccepted
819
820	END DO TimeLoop
821
822	!~~~> Succesful exit
823	IERR = 1 !~~~> The integration was successful
824
825	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
826	END SUBROUTINE ros_TlmInt
827	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
828
829
830
831	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
832	SUBROUTINE ros_SoaInt ( &
833	NSOA, Lambda, Sigma, &
834	Tstart, Tend, T, &
835	AbsTol, RelTol, &
836	ode_Fun, ode_Jac, ode_Hess, &
837	!~~~> RosenbrockSOA method coefficients
838	ros_S, ros_M, ros_E, ros_A, ros_C, &
839	ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, &
840	!~~~> Integration parameters
841	Autonomous, VectorTol, Max_no_steps, &
842	Roundoff, Hmin, Hmax, Hstart, &
843	FacMin, FacMax, FacRej, FacSafe, &
844	!~~~> Error indicator
845	IERR )
846	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
847	! Template for the implementation of a generic RosenbrockSOA method
848	! defined by ros_S (no of stages)
849	! and its coefficients ros_{A,C,M,E,Alpha,Gamma}
850	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
851
852	IMPLICIT NONE
853
854	!~~~> Input: the initial condition at Tstart; Output: the solution at T
855	INTEGER, INTENT(IN) :: NSOA
856	!~~~> First order adjoint
857	KPP_REAL, INTENT(INOUT) :: Lambda(NVAR)
858	!~~~> Second order adjoint
859	KPP_REAL, INTENT(INOUT) :: Sigma(NVAR,NSOA)
860	!~~~> Input: integration interval
861	KPP_REAL, INTENT(IN) :: Tstart,Tend
862	!~~~> Output: time at which the solution is returned (T=Tend if success)
863	KPP_REAL, INTENT(OUT) :: T
864	!~~~> Input: tolerances
865	KPP_REAL, INTENT(IN) :: AbsTol(NVAR), RelTol(NVAR)
866	!~~~> Input: ode function and its Jacobian
867	EXTERNAL ode_Fun, ode_Jac, ode_Hess
868	!~~~> Input: The RosenbrockSOA method parameters
869	INTEGER, INTENT(IN) :: ros_S
870	KPP_REAL, INTENT(IN) :: ros_M(ros_S), ros_E(ros_S), &
871	ros_Alpha(ros_S), ros_A(ros_S*(ros_S-1)/2), &
872	ros_Gamma(ros_S), ros_C(ros_S*(ros_S-1)/2), ros_ELO
873	LOGICAL, INTENT(IN) :: ros_NewF(ros_S)
874	!~~~> Input: integration parameters
875	LOGICAL, INTENT(IN) :: Autonomous, VectorTol
876	KPP_REAL, INTENT(IN) :: Hstart, Hmin, Hmax
877	INTEGER, INTENT(IN) :: Max_no_steps
878	KPP_REAL, INTENT(IN) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
879	!~~~> Output: Error indicator
880	INTEGER, INTENT(OUT) :: IERR
881	! ~~~~ Local variables
882	!KPP_REAL :: Ystage_adj(NVAR,NSOA)
883	!KPP_REAL :: dFdT(NVAR)
884	KPP_REAL :: Ystage(NVARros_S), K(NVARros_S)
885	KPP_REAL :: Ystage_tlm(NVARros_S,NSOA), K_tlm(NVARros_S,NSOA)
886	KPP_REAL :: U(NVARros_S), V(NVARros_S)
887	KPP_REAL :: W(NVARros_S,NSOA), Z(NVARros_S,NSOA)
888	KPP_REAL :: Jac(LU_NONZERO), dJdT(LU_NONZERO), Ghimj(LU_NONZERO)
889	KPP_REAL :: Hes0(NHESS), Hes1(NHESS), dHdT(NHESS)
890	KPP_REAL :: Tmp(NVAR), Tmp2(NVAR)
891	KPP_REAL :: H, HC, HA, Tau
892	INTEGER :: Pivot(NVAR), Direction, i, ioffset, joffset
893	INTEGER :: msoa, j, istage
894	!~~~> Local parameters
895	KPP_REAL, PARAMETER :: ZERO = 0.0d0, ONE = 1.0d0
896	KPP_REAL, PARAMETER :: DeltaMin = 1.0d-5
897	!~~~> Locally called functions
898	! KPP_REAL WLAMCH
899	! EXTERNAL WLAMCH
900	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
901
902
903
904	IF (Tend >= Tstart) THEN
905	Direction = +1
906	ELSE
907	Direction = -1
908	END IF
909
910	!~~~> Time loop begins below
911	TimeLoop: DO WHILE ( stack_ptr > 0 )
912
913	!~~~> Recover checkpoints for stage values and vectors
914	CALL ros_DPop( ros_S, NSOA, T, H, Ystage, K, Ystage_tlm, K_tlm )
915
916	Nstp = Nstp+1
917
918	!~~~> Compute LU decomposition
919	CALL ode_Jac(T,Ystage(1:NVAR),Ghimj)
920	CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1)
921	Tau = ONE/(DirectionHros_Gamma(1))
922	DO i=1,NVAR
923	Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+Tau
924	END DO
925	CALL ros_Decomp( Ghimj, Pivot, j )
926
927	!~~~> Compute Hessian at the beginning of the interval
928	CALL ode_Hess(T,Ystage(1),Hes0)
929
930	!~~~> Compute the stages
931	Stage: DO istage = ros_S, 1, -1
932
933	!~~~> Current istage offset.
934	ioffset = NVAR*(istage-1)
935
936	!~~~> Compute U
937	CALL WCOPY(NVAR,Lambda,1,U(ioffset+1:ioffset+NVAR),1)
938	CALL WSCAL(NVAR,ros_M(istage),U(ioffset+1:ioffset+NVAR),1)
939	DO j = istage+1, ros_S
940	joffset = NVAR*(j-1)
941	HA = ros_A((j-1)*(j-2)/2+istage)
942	HC = ros_C((j-1)(j-2)/2+istage)/(DirectionH)
943	CALL WAXPY(NVAR,HA,V(joffset+1:joffset+NVAR),1,U(ioffset+1:ioffset+NVAR),1)
944	CALL WAXPY(NVAR,HC,U(joffset+1:joffset+NVAR),1,U(ioffset+1:ioffset+NVAR),1)
945	END DO
946	CALL ros_Solve('T', Ghimj, Pivot, U(ioffset+1:ioffset+NVAR))
947	!~~~> Compute W
948	DO msoa = 1, NSOA
949	CALL WCOPY(NVAR,Sigma(1:NVAR,msoa),1,W(ioffset+1:ioffset+NVAR,msoa),1)
950	CALL WSCAL(NVAR,ros_M(istage),W(ioffset+1:ioffset+NVAR,msoa),1)
951	END DO
952	DO j = istage+1, ros_S
953	joffset = NVAR*(j-1)
954	HA = ros_A((j-1)*(j-2)/2+istage)
955	HC = ros_C((j-1)(j-2)/2+istage)/(DirectionH)
956	DO msoa = 1, NSOA
957	CALL WAXPY(NVAR,HA, &
958	Z(joffset+1:joffset+NVAR,msoa),1,W(ioffset+1:ioffset+NVAR,msoa),1)
959	CALL WAXPY(NVAR,HC, &
960	W(joffset+1:joffset+NVAR,msoa),1,W(ioffset+1:ioffset+NVAR,msoa),1)
961	END DO
962	END DO
963	DO msoa = 1, NSOA
964	CALL HessTR_Vec( Hes0, U(ioffset+1:ioffset+NVAR), Ystage_tlm(ioffset+1:ioffset+NVAR,msoa), Tmp )
965	CALL WAXPY(NVAR,ONE,Tmp,1,W(ioffset+1:ioffset+NVAR,msoa),1)
966	CALL ros_Solve('T', Ghimj, Pivot, W(ioffset+1:ioffset+NVAR,msoa))
967	END DO
968	!~~~> Compute V
969	Tau = T + ros_Alpha(istage)DirectionH
970	CALL ode_Jac(Tau,Ystage(ioffset+1:ioffset+NVAR),Jac)
971	CALL JacTR_SP_Vec(Jac,U(ioffset+1:ioffset+NVAR),V(ioffset+1:ioffset+NVAR))
972	!~~~> Compute Z
973	CALL ode_Hess(T,Ystage(ioffset+1:ioffset+NVAR),Hes1)
974	DO msoa = 1, NSOA
975	CALL JacTR_SP_Vec(Jac,W(ioffset+1:ioffset+NVAR,msoa),Z(ioffset+1:ioffset+NVAR,msoa))
976	CALL HessTR_Vec( Hes1, U(ioffset+1:ioffset+NVAR), Ystage_tlm(ioffset+1:ioffset+NVAR,msoa), Tmp )
977	CALL WAXPY(NVAR,ONE,Tmp,1,Z(ioffset+1:ioffset+NVAR,msoa),1)
978	END DO
979
980	END DO Stage
981
982	IF (.NOT.Autonomous) THEN
983	!~~~> Compute the Jacobian derivative with respect to T.
984	! Last "Jac" computed for stage 1
985	CALL ros_JacTimeDerivative ( T, Roundoff, Ystage(1), &
986	Jac, ode_Jac, dJdT )
987	!~~~> Compute the Hessian derivative with respect to T.
988	! Last "Jac" computed for stage 1
989	CALL ros_HesTimeDerivative ( T, Roundoff, Ystage(1), &
990	Hes0, ode_Hess, dHdT )
991	END IF
992
993	!~~~> Compute the new solution
994
995	!~~~> Compute Lambda
996	DO istage=1,ros_S
997	ioffset = NVAR*(istage-1)
998	! Add V_i
999	CALL WAXPY(NVAR,ONE,V(ioffset+1:ioffset+NVAR),1,Lambda,1)
1000	! Add (H0xK_i)^T * U_i
1001	CALL HessTR_Vec ( Hes0, U(ioffset+1:ioffset+NVAR), K(ioffset+1:ioffset+NVAR), Tmp )
1002	CALL WAXPY(NVAR,ONE,Tmp,1,Lambda,1)
1003	END DO
1004	! Add H * dJac_dT_0^T * \sum(gamma_i U_i)
1005	! Tmp holds sum gamma_i U_i
1006	IF (.NOT.Autonomous) THEN
1007	Tmp(1:NVAR) = ZERO
1008	DO istage = 1, ros_S
1009	ioffset = NVAR*(istage-1)
1010	CALL WAXPY(NVAR,ros_Gamma(istage),U(ioffset+1:ioffset+NVAR),1,Tmp,1)
1011	END DO
1012	CALL JacTR_SP_Vec(dJdT,Tmp,Tmp2)
1013	CALL WAXPY(NVAR,H,Tmp2,1,Lambda,1)
1014	END IF ! .NOT.Autonomous
1015
1016	!~~~> Compute Sigma
1017	DO msoa = 1, NSOA
1018
1019	DO istage=1,ros_S
1020	ioffset = NVAR*(istage-1)
1021	! Add Z_i
1022	CALL WAXPY(NVAR,ONE,Z(ioffset+1:ioffset+NVAR,msoa),1,Sigma(1:NVAR,msoa),1)
1023	! Add (Hess_0 x K_i)^T * W_i
1024	CALL HessTR_Vec ( Hes0, W(ioffset+1:ioffset+NVAR,msoa), K(ioffset+1:ioffset+NVAR), Tmp )
1025	CALL WAXPY(NVAR,ONE,Tmp,1,Sigma(1:NVAR,msoa),1)
1026	! Add (Hess_0 x K_tlm_i)^T * U_i
1027	CALL HessTR_Vec ( Hes0, U(ioffset+1:ioffset+NVAR), K_tlm(ioffset+1:ioffset+NVAR,msoa), Tmp )
1028	CALL WAXPY(NVAR,ONE,Tmp,1,Sigma(1:NVAR,msoa),1)
1029	END DO
1030
1031	!~~~> Add high derivative terms
1032	DO istage=1,ros_S
1033	ioffset = NVAR*(istage-1)
1034	CALL ros_HighDerivative ( T, Roundoff, Ystage(1), Hes0, K(ioffset+1:ioffset+NVAR), &
1035	U(ioffset+1:ioffset+NVAR), Ystage_tlm(1:NVAR,msoa), ode_Hess, Tmp)
1036	CALL WAXPY(NVAR,ONE,Tmp,1,Sigma(1:NVAR,msoa),1)
1037	END DO
1038
1039	IF (.NOT.Autonomous) THEN
1040	! Add H * dJac_dT_0^T * \sum(gamma_i W_i)
1041	! Tmp holds sum gamma_i W_i
1042	Tmp(1:NVAR) = ZERO
1043	DO istage = 1, ros_S
1044	ioffset = NVAR*(istage-1)
1045	CALL WAXPY(NVAR,ros_Gamma(istage),W(ioffset+1:ioffset+NVAR,msoa),1,Tmp,1)
1046	END DO
1047	CALL JacTR_SP_Vec(dJdT,Tmp,Tmp2)
1048	CALL WAXPY(NVAR,H,Tmp2,1,Sigma(1:NVAR,msoa),1)
1049	! Add H * ( dHess_dT_0 x Y_tlm_0)^T * \sum(gamma_i U_i)
1050	! Tmp holds sum gamma_i U_i
1051	Tmp(1:NVAR) = ZERO
1052	DO istage = 1, ros_S
1053	ioffset = NVAR*(istage-1)
1054	CALL WAXPY(NVAR,ros_Gamma(istage),U(ioffset+1:ioffset+NVAR),1,Tmp,1)
1055	END DO
1056	CALL HessTR_Vec ( dHdT, Tmp, Ystage_tlm(ioffset+1:ioffset+NVAR,msoa), Tmp2 )
1057	CALL WAXPY(NVAR,H,Tmp2,1,Sigma(1:NVAR,msoa),1)
1058	END IF ! .NOT.Autonomous
1059
1060	END DO ! msoa
1061
1062
1063	END DO TimeLoop
1064
1065	!~~~> Save last state
1066
1067	!~~~> Succesful exit
1068	IERR = 1 !~~~> The integration was successful
1069
1070	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1071	END SUBROUTINE ros_SoaInt
1072	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1073
1074
1075
1076	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1077	KPP_REAL FUNCTION ros_ErrorNorm ( Y, Ynew, Yerr, &
1078	AbsTol, RelTol, VectorTol )
1079	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1080	!~~~> Computes the "scaled norm" of the error vector Yerr
1081	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1082	IMPLICIT NONE
1083
1084	! Input arguments
1085	KPP_REAL, INTENT(IN) :: Y(NVAR), Ynew(NVAR), &
1086	Yerr(NVAR), AbsTol(NVAR), RelTol(NVAR)
1087	LOGICAL, INTENT(IN) :: VectorTol
1088	! Local variables
1089	KPP_REAL :: Err, Scale, Ymax
1090	INTEGER :: i
1091	KPP_REAL, PARAMETER :: ZERO = 0.0d0
1092
1093	Err = ZERO
1094	DO i=1,NVAR
1095	Ymax = MAX(ABS(Y(i)),ABS(Ynew(i)))
1096	IF (VectorTol) THEN
1097	Scale = AbsTol(i)+RelTol(i)*Ymax
1098	ELSE
1099	Scale = AbsTol(1)+RelTol(1)*Ymax
1100	END IF
1101	Err = Err+(Yerr(i)/Scale)**2
1102	END DO
1103	Err = SQRT(Err/NVAR)
1104
1105	ros_ErrorNorm = MAX(Err,1.0d-10)
1106
1107	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1108	END FUNCTION ros_ErrorNorm
1109	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1110
1111
1112	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1113	SUBROUTINE ros_FunTimeDerivative ( T, Roundoff, Y, &
1114	Fcn0, ode_Fun, dFdT )
1115	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1116	!~~~> The time partial derivative of the function by finite differences
1117	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1118	IMPLICIT NONE
1119
1120	!~~~> Input arguments
1121	KPP_REAL, INTENT(IN) :: T, Roundoff, Y(NVAR), Fcn0(NVAR)
1122	EXTERNAL ode_Fun
1123	!~~~> Output arguments
1124	KPP_REAL, INTENT(OUT) :: dFdT(NVAR)
1125	!~~~> Local variables
1126	KPP_REAL :: Delta
1127	KPP_REAL, PARAMETER :: ONE = 1.0d0, DeltaMin = 1.0d-6
1128
1129	Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T))
1130	CALL ode_Fun(T+Delta,Y,dFdT)
1131	CALL WAXPY(NVAR,(-ONE),Fcn0,1,dFdT,1)
1132	CALL WSCAL(NVAR,(ONE/Delta),dFdT,1)
1133
1134	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1135	END SUBROUTINE ros_FunTimeDerivative
1136	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1137
1138
1139	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1140	SUBROUTINE ros_JacTimeDerivative ( T, Roundoff, Y, &
1141	Jac0, ode_Jac, dJdT )
1142	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1143	!~~~> The time partial derivative of the Jacobian by finite differences
1144	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1145	IMPLICIT NONE
1146
1147	!~~~> Input arguments
1148	KPP_REAL, INTENT(IN) :: T, Roundoff, Y(NVAR), Jac0(LU_NONZERO)
1149	EXTERNAL ode_Jac
1150	!~~~> Output arguments
1151	KPP_REAL, INTENT(OUT) :: dJdT(LU_NONZERO)
1152	!~~~> Local variables
1153	KPP_REAL Delta
1154	KPP_REAL, PARAMETER :: ONE = 1.0d0, DeltaMin = 1.0d-6
1155
1156	Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T))
1157	CALL ode_Jac( T+Delta, Y, dJdT )
1158	CALL WAXPY(LU_NONZERO,(-ONE),Jac0,1,dJdT,1)
1159	CALL WSCAL(LU_NONZERO,(ONE/Delta),dJdT,1)
1160
1161	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1162	END SUBROUTINE ros_JacTimeDerivative
1163	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1164
1165
1166	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1167	SUBROUTINE ros_HesTimeDerivative ( T, Roundoff, Y, Hes0, ode_Hess, dHdT )
1168	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1169	!~~~> The time partial derivative of the Hessian by finite differences
1170	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1171	IMPLICIT NONE
1172
1173	!~~~> Input arguments
1174	KPP_REAL, INTENT(IN) :: T, Roundoff, Y(NVAR), Hes0(NHESS)
1175	EXTERNAL ode_Hess
1176	!~~~> Output arguments
1177	KPP_REAL, INTENT(OUT) :: dHdT(NHESS)
1178	!~~~> Local variables
1179	KPP_REAL Delta
1180	KPP_REAL, PARAMETER :: ONE = 1.0d0, DeltaMin = 1.0d-6
1181
1182	Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T))
1183	CALL ode_Hess( T+Delta, Y, dHdT )
1184	CALL WAXPY(NHESS,(-ONE),Hes0,1,dHdT,1)
1185	CALL WSCAL(NHESS,(ONE/Delta),dHdT,1)
1186
1187	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1188	END SUBROUTINE ros_HesTimeDerivative
1189	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1190
1191
1192	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1193	SUBROUTINE ros_HighDerivative ( T, Roundoff, Y, Hes0, K, U, Y_tlm, &
1194	ode_Hess, Term)
1195	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1196	!~~~> High order derivative by finite differences:
1197	! d/dy { (Hes0 x K_i)^T * U_i } * Y_tlm
1198	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1199	IMPLICIT NONE
1200
1201	!~~~> Input arguments
1202	KPP_REAL, INTENT(IN) :: T, Roundoff, Y(NVAR), Hes0(NHESS)
1203	KPP_REAL, INTENT(IN) :: K(NVAR), U(NVAR), Y_tlm(NVAR)
1204	EXTERNAL ode_Hess
1205	!~~~> Output arguments
1206	KPP_REAL, INTENT(OUT) :: Term(NVAR)
1207	!~~~> Local variables
1208	KPP_REAL :: Delta, Y1(NVAR), Hes1(NHESS), Tmp(NVAR)
1209	KPP_REAL, PARAMETER :: ONE = 1.0d0, DeltaMin = 1.0d-6
1210
1211	CALL HessTR_Vec ( Hes0, U, K, Tmp )
1212
1213	Delta = SQRT(Roundoff)*MAX(DeltaMin,ABS(T))
1214	Y1(1:NVAR) = Y(1:NVAR) + Delta*Y_tlm(1:NVAR)
1215	CALL ode_Hess( T, Y1, Hes1 )
1216	! Add (Hess_0 x K_i)^T * U_i
1217	CALL HessTR_Vec ( Hes1, U, K, Term )
1218
1219	CALL WAXPY(NVAR,(-ONE),Tmp,1,Term,1)
1220	CALL WSCAL(NVAR,(ONE/Delta),Term,1)
1221
1222	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1223	END SUBROUTINE ros_HighDerivative
1224	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1225
1226
1227	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1228	SUBROUTINE ros_PrepareMatrix ( H, Direction, gam, &
1229	Jac0, Ghimj, Pivot, Singular )
1230	! --- --- --- --- --- --- --- --- --- --- --- --- ---
1231	! Prepares the LHS matrix for stage calculations
1232	! 1. Construct Ghimj = 1/(H*ham) - Jac0
1233	! "(Gamma H) Inverse Minus Jacobian"
1234	! 2. Repeat LU decomposition of Ghimj until successful.
1235	! -half the step size if LU decomposition fails and retry
1236	! -exit after 5 consecutive fails
1237	! --- --- --- --- --- --- --- --- --- --- --- --- ---
1238	IMPLICIT NONE
1239
1240	!~~~> Input arguments
1241	KPP_REAL, INTENT(IN) :: gam, Jac0(LU_NONZERO)
1242	INTEGER, INTENT(IN) :: Direction
1243	!~~~> Output arguments
1244	KPP_REAL, INTENT(OUT) :: Ghimj(LU_NONZERO)
1245	LOGICAL, INTENT(OUT) :: Singular
1246	INTEGER, INTENT(OUT) :: Pivot(NVAR)
1247	!~~~> Inout arguments
1248	KPP_REAL, INTENT(INOUT) :: H ! step size is decreased when LU fails
1249	!~~~> Local variables
1250	INTEGER :: i, ising, Nconsecutive
1251	KPP_REAL :: ghinv
1252	KPP_REAL, PARAMETER :: ONE = 1.0d0, HALF = 0.5d0
1253
1254	Nconsecutive = 0
1255	Singular = .TRUE.
1256
1257	DO WHILE (Singular)
1258
1259	!~~~> Construct Ghimj = 1/(H*ham) - Jac0
1260	CALL WCOPY(LU_NONZERO,Jac0,1,Ghimj,1)
1261	CALL WSCAL(LU_NONZERO,(-ONE),Ghimj,1)
1262	ghinv = ONE/(DirectionHgam)
1263	DO i=1,NVAR
1264	Ghimj(LU_DIAG(i)) = Ghimj(LU_DIAG(i))+ghinv
1265	END DO
1266	!~~~> Compute LU decomposition
1267	CALL ros_Decomp( Ghimj, Pivot, ising )
1268	IF (ising == 0) THEN
1269	!~~~> If successful done
1270	Singular = .FALSE.
1271	ELSE ! ising .ne. 0
1272	!~~~> If unsuccessful half the step size; if 5 consecutive fails then return
1273	Nsng = Nsng+1
1274	Nconsecutive = Nconsecutive+1
1275	Singular = .TRUE.
1276	PRINT*,'Warning: LU Decomposition returned ising = ',ising
1277	IF (Nconsecutive <= 5) THEN ! Less than 5 consecutive failed decomps
1278	H = H*HALF
1279	ELSE ! More than 5 consecutive failed decompositions
1280	RETURN
1281	END IF ! Nconsecutive
1282	END IF ! ising
1283
1284	END DO ! WHILE Singular
1285
1286	END SUBROUTINE ros_PrepareMatrix
1287
1288
1289	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1290	SUBROUTINE ros_Decomp( A, Pivot, ising )
1291	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1292	! Template for the LU decomposition
1293	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1294	IMPLICIT NONE
1295	!~~~> Inout variables
1296	KPP_REAL, INTENT(INOUT) :: A(LU_NONZERO)
1297	!~~~> Output variables
1298	INTEGER, INTENT(OUT) :: Pivot(NVAR), ising
1299
1300	CALL KppDecomp ( A, ising )
1301	!~~~> Note: for a full matrix use Lapack:
1302	! CALL DGETRF( NVAR, NVAR, A, NVAR, Pivot, ising )
1303	Pivot(1) = 1
1304
1305	Ndec = Ndec + 1
1306
1307	END SUBROUTINE ros_Decomp
1308
1309
1310	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1311	SUBROUTINE ros_Solve( C, A, Pivot, b )
1312	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1313	! Template for the forward/backward substitution (using pre-computed LU decomp)
1314	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1315	IMPLICIT NONE
1316	!~~~> Input variables
1317	CHARACTER, INTENT(IN) :: C
1318	KPP_REAL, INTENT(IN) :: A(LU_NONZERO)
1319	INTEGER, INTENT(IN) :: Pivot(NVAR)
1320	!~~~> InOut variables
1321	KPP_REAL, INTENT(INOUT) :: b(NVAR)
1322
1323	SELECT CASE (C)
1324	CASE ('N')
1325	CALL KppSolve( A, b )
1326	CASE ('T')
1327	CALL KppSolveTR( A, b, b )
1328	CASE DEFAULT
1329	PRINT*,'Unknown C = (',C,') in ros_Solve'
1330	STOP
1331	END SELECT
1332	!~~~> Note: for a full matrix use Lapack:
1333	! NRHS = 1
1334	! CALL DGETRS( C, NVAR , NRHS, A, NVAR, Pivot, b, NVAR, INFO )
1335
1336	Nsol = Nsol+1
1337
1338	END SUBROUTINE ros_Solve
1339
1340
1341	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1342	SUBROUTINE Ros2 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
1343	ros_Gamma,ros_NewF,ros_ELO,ros_Name)
1344	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1345	! --- AN L-STABLE METHOD, 2 stages, order 2
1346	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1347
1348	IMPLICIT NONE
1349
1350	INTEGER, PARAMETER :: S = 2
1351	INTEGER, INTENT(OUT) :: ros_S
1352	KPP_REAL, DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
1353	KPP_REAL, DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
1354	KPP_REAL, INTENT(OUT) :: ros_ELO
1355	LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
1356	CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
1357	KPP_REAL :: g
1358
1359	g = 1.0d0 + 1.0d0/SQRT(2.0d0)
1360
1361	!~~~> Name of the method
1362	ros_Name = 'ROS-2'
1363	!~~~> Number of stages
1364	ros_S = S
1365
1366	!~~~> The coefficient matrices A and C are strictly lower triangular.
1367	! The lower triangular (subdiagonal) elements are stored in row-wise order:
1368	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1369	! The general mapping formula is:
1370	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
1371	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
1372
1373	ros_A(1) = (1.d0)/g
1374	ros_C(1) = (-2.d0)/g
1375	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1376	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1377	ros_NewF(1) = .TRUE.
1378	ros_NewF(2) = .TRUE.
1379	!~~~> M_i = Coefficients for new step solution
1380	ros_M(1)= (3.d0)/(2.d0*g)
1381	ros_M(2)= (1.d0)/(2.d0*g)
1382	! E_i = Coefficients for error estimator
1383	ros_E(1) = 1.d0/(2.d0*g)
1384	ros_E(2) = 1.d0/(2.d0*g)
1385	!~~~> ros_ELO = estimator of local order - the minimum between the
1386	! main and the embedded scheme orders plus one
1387	ros_ELO = 2.0d0
1388	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1389	ros_Alpha(1) = 0.0d0
1390	ros_Alpha(2) = 1.0d0
1391	!~~~> Gamma_i = \sum_j gamma_{i,j}
1392	ros_Gamma(1) = g
1393	ros_Gamma(2) =-g
1394
1395	END SUBROUTINE Ros2
1396
1397
1398	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1399	SUBROUTINE Ros3 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
1400	ros_Gamma,ros_NewF,ros_ELO,ros_Name)
1401	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1402	! --- AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations
1403	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1404
1405	IMPLICIT NONE
1406
1407	INTEGER, PARAMETER :: S = 3
1408	INTEGER, INTENT(OUT) :: ros_S
1409	KPP_REAL, DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
1410	KPP_REAL, DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
1411	KPP_REAL, INTENT(OUT) :: ros_ELO
1412	LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
1413	CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
1414
1415	!~~~> Name of the method
1416	ros_Name = 'ROS-3'
1417	!~~~> Number of stages
1418	ros_S = S
1419
1420	!~~~> The coefficient matrices A and C are strictly lower triangular.
1421	! The lower triangular (subdiagonal) elements are stored in row-wise order:
1422	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1423	! The general mapping formula is:
1424	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
1425	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
1426
1427	ros_A(1)= 1.d0
1428	ros_A(2)= 1.d0
1429	ros_A(3)= 0.d0
1430
1431	ros_C(1) = -0.10156171083877702091975600115545d+01
1432	ros_C(2) = 0.40759956452537699824805835358067d+01
1433	ros_C(3) = 0.92076794298330791242156818474003d+01
1434	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1435	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1436	ros_NewF(1) = .TRUE.
1437	ros_NewF(2) = .TRUE.
1438	ros_NewF(3) = .FALSE.
1439	!~~~> M_i = Coefficients for new step solution
1440	ros_M(1) = 0.1d+01
1441	ros_M(2) = 0.61697947043828245592553615689730d+01
1442	ros_M(3) = -0.42772256543218573326238373806514d+00
1443	! E_i = Coefficients for error estimator
1444	ros_E(1) = 0.5d+00
1445	ros_E(2) = -0.29079558716805469821718236208017d+01
1446	ros_E(3) = 0.22354069897811569627360909276199d+00
1447	!~~~> ros_ELO = estimator of local order - the minimum between the
1448	! main and the embedded scheme orders plus 1
1449	ros_ELO = 3.0d0
1450	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1451	ros_Alpha(1)= 0.0d+00
1452	ros_Alpha(2)= 0.43586652150845899941601945119356d+00
1453	ros_Alpha(3)= 0.43586652150845899941601945119356d+00
1454	!~~~> Gamma_i = \sum_j gamma_{i,j}
1455	ros_Gamma(1)= 0.43586652150845899941601945119356d+00
1456	ros_Gamma(2)= 0.24291996454816804366592249683314d+00
1457	ros_Gamma(3)= 0.21851380027664058511513169485832d+01
1458
1459	END SUBROUTINE Ros3
1460
1461	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1462
1463
1464	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1465	SUBROUTINE Ros4 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
1466	ros_Gamma,ros_NewF,ros_ELO,ros_Name)
1467	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1468	! L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES
1469	! L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
1470	!
1471	! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
1472	! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
1473	! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
1474	! SPRINGER-VERLAG (1990)
1475	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1476
1477	IMPLICIT NONE
1478
1479	INTEGER, PARAMETER :: S=4
1480	INTEGER, INTENT(OUT) :: ros_S
1481	KPP_REAL, DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
1482	KPP_REAL, DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
1483	KPP_REAL, INTENT(OUT) :: ros_ELO
1484	LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
1485	CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
1486
1487
1488	!~~~> Name of the method
1489	ros_Name = 'ROS-4'
1490	!~~~> Number of stages
1491	ros_S = S
1492
1493	!~~~> The coefficient matrices A and C are strictly lower triangular.
1494	! The lower triangular (subdiagonal) elements are stored in row-wise order:
1495	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1496	! The general mapping formula is:
1497	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
1498	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
1499
1500	ros_A(1) = 0.2000000000000000d+01
1501	ros_A(2) = 0.1867943637803922d+01
1502	ros_A(3) = 0.2344449711399156d+00
1503	ros_A(4) = ros_A(2)
1504	ros_A(5) = ros_A(3)
1505	ros_A(6) = 0.0D0
1506
1507	ros_C(1) =-0.7137615036412310d+01
1508	ros_C(2) = 0.2580708087951457d+01
1509	ros_C(3) = 0.6515950076447975d+00
1510	ros_C(4) =-0.2137148994382534d+01
1511	ros_C(5) =-0.3214669691237626d+00
1512	ros_C(6) =-0.6949742501781779d+00
1513	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1514	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1515	ros_NewF(1) = .TRUE.
1516	ros_NewF(2) = .TRUE.
1517	ros_NewF(3) = .TRUE.
1518	ros_NewF(4) = .FALSE.
1519	!~~~> M_i = Coefficients for new step solution
1520	ros_M(1) = 0.2255570073418735d+01
1521	ros_M(2) = 0.2870493262186792d+00
1522	ros_M(3) = 0.4353179431840180d+00
1523	ros_M(4) = 0.1093502252409163d+01
1524	!~~~> E_i = Coefficients for error estimator
1525	ros_E(1) =-0.2815431932141155d+00
1526	ros_E(2) =-0.7276199124938920d-01
1527	ros_E(3) =-0.1082196201495311d+00
1528	ros_E(4) =-0.1093502252409163d+01
1529	!~~~> ros_ELO = estimator of local order - the minimum between the
1530	! main and the embedded scheme orders plus 1
1531	ros_ELO = 4.0d0
1532	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1533	ros_Alpha(1) = 0.D0
1534	ros_Alpha(2) = 0.1145640000000000d+01
1535	ros_Alpha(3) = 0.6552168638155900d+00
1536	ros_Alpha(4) = ros_Alpha(3)
1537	!~~~> Gamma_i = \sum_j gamma_{i,j}
1538	ros_Gamma(1) = 0.5728200000000000d+00
1539	ros_Gamma(2) =-0.1769193891319233d+01
1540	ros_Gamma(3) = 0.7592633437920482d+00
1541	ros_Gamma(4) =-0.1049021087100450d+00
1542
1543	END SUBROUTINE Ros4
1544
1545	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1546	SUBROUTINE Rodas3 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
1547	ros_Gamma,ros_NewF,ros_ELO,ros_Name)
1548	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1549	! --- A STIFFLY-STABLE METHOD, 4 stages, order 3
1550	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1551
1552	IMPLICIT NONE
1553
1554	INTEGER, PARAMETER :: S=4
1555	INTEGER, INTENT(OUT) :: ros_S
1556	KPP_REAL, DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
1557	KPP_REAL, DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
1558	KPP_REAL, INTENT(OUT) :: ros_ELO
1559	LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
1560	CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
1561
1562
1563	!~~~> Name of the method
1564	ros_Name = 'RODAS-3'
1565	!~~~> Number of stages
1566	ros_S = S
1567
1568	!~~~> The coefficient matrices A and C are strictly lower triangular.
1569	! The lower triangular (subdiagonal) elements are stored in row-wise order:
1570	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1571	! The general mapping formula is:
1572	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
1573	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
1574
1575	ros_A(1) = 0.0d+00
1576	ros_A(2) = 2.0d+00
1577	ros_A(3) = 0.0d+00
1578	ros_A(4) = 2.0d+00
1579	ros_A(5) = 0.0d+00
1580	ros_A(6) = 1.0d+00
1581
1582	ros_C(1) = 4.0d+00
1583	ros_C(2) = 1.0d+00
1584	ros_C(3) =-1.0d+00
1585	ros_C(4) = 1.0d+00
1586	ros_C(5) =-1.0d+00
1587	ros_C(6) =-(8.0d+00/3.0d+00)
1588
1589	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1590	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1591	ros_NewF(1) = .TRUE.
1592	ros_NewF(2) = .FALSE.
1593	ros_NewF(3) = .TRUE.
1594	ros_NewF(4) = .TRUE.
1595	!~~~> M_i = Coefficients for new step solution
1596	ros_M(1) = 2.0d+00
1597	ros_M(2) = 0.0d+00
1598	ros_M(3) = 1.0d+00
1599	ros_M(4) = 1.0d+00
1600	!~~~> E_i = Coefficients for error estimator
1601	ros_E(1) = 0.0d+00
1602	ros_E(2) = 0.0d+00
1603	ros_E(3) = 0.0d+00
1604	ros_E(4) = 1.0d+00
1605	!~~~> ros_ELO = estimator of local order - the minimum between the
1606	! main and the embedded scheme orders plus 1
1607	ros_ELO = 3.0d+00
1608	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1609	ros_Alpha(1) = 0.0d+00
1610	ros_Alpha(2) = 0.0d+00
1611	ros_Alpha(3) = 1.0d+00
1612	ros_Alpha(4) = 1.0d+00
1613	!~~~> Gamma_i = \sum_j gamma_{i,j}
1614	ros_Gamma(1) = 0.5d+00
1615	ros_Gamma(2) = 1.5d+00
1616	ros_Gamma(3) = 0.0d+00
1617	ros_Gamma(4) = 0.0d+00
1618
1619	END SUBROUTINE Rodas3
1620
1621	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1622	SUBROUTINE Rodas4 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
1623	ros_Gamma,ros_NewF,ros_ELO,ros_Name)
1624	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1625	! STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES
1626	!
1627	! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
1628	! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
1629	! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
1630	! SPRINGER-VERLAG (1996)
1631	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1632
1633	IMPLICIT NONE
1634
1635	INTEGER, PARAMETER :: S=6
1636	INTEGER, INTENT(OUT) :: ros_S
1637	KPP_REAL, DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
1638	KPP_REAL, DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
1639	KPP_REAL, INTENT(OUT) :: ros_ELO
1640	LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
1641	CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
1642
1643
1644	!~~~> Name of the method
1645	ros_Name = 'RODAS-4'
1646	!~~~> Number of stages
1647	ros_S = S
1648
1649	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1650	ros_Alpha(1) = 0.000d0
1651	ros_Alpha(2) = 0.386d0
1652	ros_Alpha(3) = 0.210d0
1653	ros_Alpha(4) = 0.630d0
1654	ros_Alpha(5) = 1.000d0
1655	ros_Alpha(6) = 1.000d0
1656
1657	!~~~> Gamma_i = \sum_j gamma_{i,j}
1658	ros_Gamma(1) = 0.2500000000000000d+00
1659	ros_Gamma(2) =-0.1043000000000000d+00
1660	ros_Gamma(3) = 0.1035000000000000d+00
1661	ros_Gamma(4) =-0.3620000000000023d-01
1662	ros_Gamma(5) = 0.0d0
1663	ros_Gamma(6) = 0.0d0
1664
1665	!~~~> The coefficient matrices A and C are strictly lower triangular.
1666	! The lower triangular (subdiagonal) elements are stored in row-wise order:
1667	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1668	! The general mapping formula is: A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
1669	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
1670
1671	ros_A(1) = 0.1544000000000000d+01
1672	ros_A(2) = 0.9466785280815826d+00
1673	ros_A(3) = 0.2557011698983284d+00
1674	ros_A(4) = 0.3314825187068521d+01
1675	ros_A(5) = 0.2896124015972201d+01
1676	ros_A(6) = 0.9986419139977817d+00
1677	ros_A(7) = 0.1221224509226641d+01
1678	ros_A(8) = 0.6019134481288629d+01
1679	ros_A(9) = 0.1253708332932087d+02
1680	ros_A(10) =-0.6878860361058950d+00
1681	ros_A(11) = ros_A(7)
1682	ros_A(12) = ros_A(8)
1683	ros_A(13) = ros_A(9)
1684	ros_A(14) = ros_A(10)
1685	ros_A(15) = 1.0d+00
1686
1687	ros_C(1) =-0.5668800000000000d+01
1688	ros_C(2) =-0.2430093356833875d+01
1689	ros_C(3) =-0.2063599157091915d+00
1690	ros_C(4) =-0.1073529058151375d+00
1691	ros_C(5) =-0.9594562251023355d+01
1692	ros_C(6) =-0.2047028614809616d+02
1693	ros_C(7) = 0.7496443313967647d+01
1694	ros_C(8) =-0.1024680431464352d+02
1695	ros_C(9) =-0.3399990352819905d+02
1696	ros_C(10) = 0.1170890893206160d+02
1697	ros_C(11) = 0.8083246795921522d+01
1698	ros_C(12) =-0.7981132988064893d+01
1699	ros_C(13) =-0.3152159432874371d+02
1700	ros_C(14) = 0.1631930543123136d+02
1701	ros_C(15) =-0.6058818238834054d+01
1702
1703	!~~~> M_i = Coefficients for new step solution
1704	ros_M(1) = ros_A(7)
1705	ros_M(2) = ros_A(8)
1706	ros_M(3) = ros_A(9)
1707	ros_M(4) = ros_A(10)
1708	ros_M(5) = 1.0d+00
1709	ros_M(6) = 1.0d+00
1710
1711	!~~~> E_i = Coefficients for error estimator
1712	ros_E(1) = 0.0d+00
1713	ros_E(2) = 0.0d+00
1714	ros_E(3) = 0.0d+00
1715	ros_E(4) = 0.0d+00
1716	ros_E(5) = 0.0d+00
1717	ros_E(6) = 1.0d+00
1718
1719	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1720	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1721	ros_NewF(1) = .TRUE.
1722	ros_NewF(2) = .TRUE.
1723	ros_NewF(3) = .TRUE.
1724	ros_NewF(4) = .TRUE.
1725	ros_NewF(5) = .TRUE.
1726	ros_NewF(6) = .TRUE.
1727
1728	!~~~> ros_ELO = estimator of local order - the minimum between the
1729	! main and the embedded scheme orders plus 1
1730	ros_ELO = 4.0d0
1731
1732	END SUBROUTINE Rodas4
1733
1734
1735	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1736	SUBROUTINE Fun_Template( T, Y, Ydot )
1737	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1738	! Template for the ODE function call.
1739	! Updates the rate coefficients (and possibly the fixed species) at each call
1740	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1741	!~~~> Input variables
1742	KPP_REAL T, Y(NVAR)
1743	!~~~> Output variables
1744	KPP_REAL Ydot(NVAR)
1745	!~~~> Local variables
1746	KPP_REAL Told
1747
1748	Told = TIME
1749	TIME = T
1750	CALL Update_SUN()
1751	CALL Update_RCONST()
1752	CALL Fun( Y, FIX, RCONST, Ydot )
1753	TIME = Told
1754
1755	Nfun = Nfun+1
1756
1757	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1758	END SUBROUTINE Fun_Template
1759	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1760
1761
1762	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1763	SUBROUTINE Jac_Template( T, Y, Jcb )
1764	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1765	! Template for the ODE Jacobian call.
1766	! Updates the rate coefficients (and possibly the fixed species) at each call
1767	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1768
1769	!~~~> Input variables
1770	KPP_REAL T, Y(NVAR)
1771	!~~~> Output variables
1772	KPP_REAL Jcb(LU_NONZERO)
1773	!~~~> Local variables
1774	KPP_REAL Told
1775
1776	Told = TIME
1777	TIME = T
1778	CALL Update_SUN()
1779	CALL Update_RCONST()
1780	CALL Jac_SP( Y, FIX, RCONST, Jcb )
1781	TIME = Told
1782
1783	Njac = Njac+1
1784
1785	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1786	END SUBROUTINE Jac_Template
1787	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1788
1789	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1790	SUBROUTINE Hess_Template( T, Y, Hes )
1791	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1792	! Template for the ODE Hessian call.
1793	! Updates the rate coefficients (and possibly the fixed species) at each call
1794	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1795	!~~~> Input variables
1796	KPP_REAL T, Y(NVAR)
1797	!~~~> Output variables
1798	KPP_REAL Hes(NHESS)
1799	!~~~> Local variables
1800	KPP_REAL Told
1801
1802	Told = TIME
1803	TIME = T
1804	CALL Update_SUN()
1805	CALL Update_RCONST()
1806	CALL Hessian( Y, FIX, RCONST, Hes )
1807	TIME = Told
1808
1809	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1810	END SUBROUTINE Hess_Template
1811	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1812
1813	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1814	SUBROUTINE ros_AllocateDBuffers( S )
1815	!~~~> Allocate buffer space for discrete adjoint
1816	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1817	INTEGER :: i, S
1818
1819	ALLOCATE( buf_H(bufsize), STAT=i )
1820	IF (i/=0) THEN
1821	PRINT*,'Failed allocation of buffer H'; STOP
1822	END IF
1823	ALLOCATE( buf_T(bufsize), STAT=i )
1824	IF (i/=0) THEN
1825	PRINT*,'Failed allocation of buffer T'; STOP
1826	END IF
1827	ALLOCATE( buf_Y(NVAR*S,bufsize), STAT=i )
1828	IF (i/=0) THEN
1829	PRINT*,'Failed allocation of buffer Y'; STOP
1830	END IF
1831	ALLOCATE( buf_K(NVAR*S,bufsize), STAT=i )
1832	IF (i/=0) THEN
1833	PRINT*,'Failed allocation of buffer K'; STOP
1834	END IF
1835	ALLOCATE( buf_Y_tlm(NVAR*S,bufsize), STAT=i )
1836	IF (i/=0) THEN
1837	PRINT*,'Failed allocation of buffer Y_tlm'; STOP
1838	END IF
1839	ALLOCATE( buf_K_tlm(NVAR*S,bufsize), STAT=i )
1840	IF (i/=0) THEN
1841	PRINT*,'Failed allocation of buffer K_tlm'; STOP
1842	END IF
1843
1844	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1845	END SUBROUTINE ros_AllocateDBuffers
1846	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1847
1848
1849	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1850	SUBROUTINE ros_FreeDBuffers
1851	!~~~> Dallocate buffer space for discrete adjoint
1852	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1853	INTEGER :: i
1854
1855	DEALLOCATE( buf_H, STAT=i )
1856	IF (i/=0) THEN
1857	PRINT*,'Failed deallocation of buffer H'; STOP
1858	END IF
1859	DEALLOCATE( buf_T, STAT=i )
1860	IF (i/=0) THEN
1861	PRINT*,'Failed deallocation of buffer T'; STOP
1862	END IF
1863	DEALLOCATE( buf_Y, STAT=i )
1864	IF (i/=0) THEN
1865	PRINT*,'Failed deallocation of buffer Y'; STOP
1866	END IF
1867	DEALLOCATE( buf_K, STAT=i )
1868	IF (i/=0) THEN
1869	PRINT*,'Failed deallocation of buffer K'; STOP
1870	END IF
1871	DEALLOCATE( buf_Y_tlm, STAT=i )
1872	IF (i/=0) THEN
1873	PRINT*,'Failed deallocation of buffer Y_tlm'; STOP
1874	END IF
1875	DEALLOCATE( buf_K_tlm, STAT=i )
1876	IF (i/=0) THEN
1877	PRINT*,'Failed deallocation of buffer K_tlm'; STOP
1878	END IF
1879
1880	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1881	END SUBROUTINE ros_FreeDBuffers
1882	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1883
1884
1885	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1886	SUBROUTINE ros_DPush( S, NSOA, T, H, Ystage, K, Ystage_tlm, K_tlm )
1887	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1888	!~~~> Saves the next trajectory snapshot for discrete adjoints
1889	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1890
1891	INTEGER, INTENT(IN) :: S ! no of stages
1892	INTEGER, INTENT(IN) :: NSOA ! no of second order adjoints
1893	KPP_REAL :: T, H, Ystage(NVARS), K(NVARS)
1894	KPP_REAL :: Ystage_tlm(NVARS,NSOA), K_tlm(NVARS,NSOA)
1895
1896	stack_ptr = stack_ptr + 1
1897	IF ( stack_ptr > bufsize ) THEN
1898	PRINT*,'Push failed: buffer overflow'
1899	STOP
1900	END IF
1901	buf_H( stack_ptr ) = H
1902	buf_T( stack_ptr ) = T
1903	CALL WCOPY(NVARS,Ystage,1,buf_Y(1:NVARS,stack_ptr),1)
1904	CALL WCOPY(NVARS,K,1,buf_K(1:NVARS,stack_ptr),1)
1905	CALL WCOPY(NVARSNSOA,Ystage_tlm,1,buf_Y_tlm(1:NVARSNSOA,stack_ptr),1)
1906	CALL WCOPY(NVARSNSOA,K_tlm,1,buf_K_tlm(1:NVARSNSOA,stack_ptr),1)
1907
1908	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1909	END SUBROUTINE ros_DPush
1910	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1911
1912
1913	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1914	SUBROUTINE ros_DPop( S, NSOA, T, H, Ystage, K, Ystage_tlm, K_tlm )
1915	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1916	!~~~> Retrieves the next trajectory snapshot for discrete adjoints
1917	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1918
1919	INTEGER, INTENT(IN) :: S ! no of stages
1920	INTEGER, INTENT(IN) :: NSOA ! no of second order adjoints
1921	KPP_REAL, INTENT(OUT) :: T, H, Ystage(NVARS), K(NVARS)
1922	KPP_REAL, INTENT(OUT) :: Ystage_tlm(NVARS,NSOA), K_tlm(NVARS,NSOA)
1923
1924	IF ( stack_ptr <= 0 ) THEN
1925	PRINT*,'Pop failed: empty buffer'
1926	STOP
1927	END IF
1928	H = buf_H( stack_ptr )
1929	T = buf_T( stack_ptr )
1930	CALL WCOPY(NVARS,buf_Y(1:NVARS,stack_ptr),1,Ystage,1)
1931	CALL WCOPY(NVARS,buf_K(1:NVARS,stack_ptr),1,K,1)
1932	CALL WCOPY(NVARSNSOA,buf_Y_tlm(1:NVARSNSOA,stack_ptr),1,Ystage_tlm,1)
1933	CALL WCOPY(NVARSNSOA,buf_K_tlm(1:NVARSNSOA,stack_ptr),1,K_tlm,1)
1934
1935	stack_ptr = stack_ptr - 1
1936
1937	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1938	END SUBROUTINE ros_DPop
1939	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1940
1941	END MODULE KPP_ROOT_Integrator
1942

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

Context Navigation

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

Download in other formats: