Home

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp4palm/templates/rosenbrock_vec.f90 @ 3201

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

Line
1	! NOTE: This routine is a vectorised version of the routine in
2	! rosenbrock_mz.f90 and can only be used with KP4.
3	!
4	SUBROUTINE Rosenbrock(Y,Tstart,Tend, &
5	AbsTol,RelTol, &
6	RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR)
7	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
8	!
9	! Solves the system y'=F(t,y) using a Rosenbrock method defined by:
10	!
11	! G = 1/(H*gamma(1)) - Jac(t0,Y0)
12	! T_i = t0 + Alpha(i)*H
13	! Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
14	! G * K_i = Fun( T_i, Y_i ) + \sum_{j=1}^S C(i,j)/H * K_j +
15	! gamma(i)*dF/dT(t0, Y0)
16	! Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
17	!
18	! For details on Rosenbrock methods and their implementation consult:
19	! E. Hairer and G. Wanner
20	! "Solving ODEs II. Stiff and differential-algebraic problems".
21	! Springer series in computational mathematics, Springer-Verlag, 1996.
22	! The codes contained in the book inspired this implementation.
23	!
24	! (C) Adrian Sandu, August 2004
25	! Virginia Polytechnic Institute and State University
26	! Contact: sandu@cs.vt.edu
27	! This implementation is part of KPP - the Kinetic PreProcessor
28	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
29	!
30	!~~~> INPUT ARGUMENTS:
31	!
32	!- Y(NVAR) = vector of initial conditions (at T=Tstart)
33	!- [Tstart,Tend] = time range of integration
34	! (if Tstart>Tend the integration is performed backwards in time)
35	!- RelTol, AbsTol = user precribed accuracy
36	!- SUBROUTINE Fun( T, Y, Ydot ) = ODE function,
37	! returns Ydot = Y' = F(T,Y)
38	!- SUBROUTINE Jac( T, Y, Jcb ) = Jacobian of the ODE function,
39	! returns Jcb = dFun/dY
40	!- ICNTRL(1:20) = integer inputs parameters
41	!- RCNTRL(1:20) = real inputs parameters
42	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
43	!
44	!~~~> OUTPUT ARGUMENTS:
45	!
46	!- Y(NVAR) -> vector of final states (at T->Tend)
47	!- ISTATUS(1:20) -> integer output parameters
48	!- RSTATUS(1:20) -> real output parameters
49	!- IERR -> job status upon return
50	! success (positive value) or
51	! failure (negative value)
52	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
53	!
54	!~~~> INPUT PARAMETERS:
55	!
56	! Note: For input parameters equal to zero the default values of the
57	! corresponding variables are used.
58	!
59	! ICNTRL(1) = 1: F = F(y) Independent of T (AUTONOMOUS)
60	! = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
61	!
62	! ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors
63	! = 1: AbsTol, RelTol are scalars
64	!
65	! ICNTRL(3) -> selection of a particular Rosenbrock method
66	! = 0 : default method is Rodas3
67	! = 1 : method is Ros2
68	! = 2 : method is Ros3
69	! = 3 : method is Ros4
70	! = 4 : method is Rodas3
71	! = 5: method is Rodas4
72	!
73	! ICNTRL(4) -> maximum number of integration steps
74	! For ICNTRL(4)=0) the default value of 100000 is used
75	!
76	! RCNTRL(1) -> Hmin, lower bound for the integration step size
77	! It is strongly recommended to keep Hmin = ZERO
78	! RCNTRL(2) -> Hmax, upper bound for the integration step size
79	! RCNTRL(3) -> Hstart, starting value for the integration step size
80	!
81	! RCNTRL(4) -> FacMin, lower bound on step decrease factor (default=0.2)
82	! RCNTRL(5) -> FacMax, upper bound on step increase factor (default=6)
83	! RCNTRL(6) -> FacRej, step decrease factor after multiple rejections
84	! (default=0.1)
85	! RCNTRL(7) -> FacSafe, by which the new step is slightly smaller
86	! than the predicted value (default=0.9)
87	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
88	!
89	!~~~> OUTPUT PARAMETERS:
90	!
91	! Note: each call to Rosenbrock adds the current no. of fcn calls
92	! to previous value of ISTATUS(1), and similar for the other params.
93	! Set ISTATUS(1:20) = 0 before call to avoid this accumulation.
94	!
95	! ISTATUS(1) = No. of function calls
96	! ISTATUS(2) = No. of jacobian calls
97	! ISTATUS(3) = No. of steps
98	! ISTATUS(4) = No. of accepted steps
99	! ISTATUS(5) = No. of rejected steps (except at the beginning)
100	! ISTATUS(6) = No. of LU decompositions
101	! ISTATUS(7) = No. of forward/backward substitutions
102	! ISTATUS(8) = No. of singular matrix decompositions
103	!
104	! RSTATUS(1) -> Texit, the time corresponding to the
105	! computed Y upon return
106	! RSTATUS(2) -> Hexit, last accepted step before exit
107	! For multiple restarts, use Hexit as Hstart in the following run
108	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
109
110	IMPLICIT NONE
111
112	!~~~> Arguments
113	REAL(kind=dp), INTENT(INOUT) :: Y(:,:)
114	REAL(kind=dp), INTENT(IN) :: Tstart,Tend
115	REAL(kind=dp), INTENT(IN) :: AbsTol(NVAR),RelTol(NVAR)
116	INTEGER, INTENT(IN) :: ICNTRL(20)
117	REAL(kind=dp), INTENT(IN) :: RCNTRL(20)
118	INTEGER, INTENT(INOUT) :: ISTATUS(20)
119	REAL(kind=dp), INTENT(INOUT) :: RSTATUS(20)
120	INTEGER, INTENT(OUT) :: IERR
121	!~~~> The method parameters
122	INTEGER, PARAMETER :: Smax = 6
123	INTEGER :: Method, ros_S
124	REAL(kind=dp), DIMENSION(Smax) :: ros_M, ros_E, ros_Alpha, ros_Gamma
125	REAL(kind=dp), DIMENSION(Smax*(Smax-1)/2) :: ros_A, ros_C
126	REAL(kind=dp) :: ros_ELO
127	LOGICAL, DIMENSION(Smax) :: ros_NewF
128	CHARACTER(LEN=12) :: ros_Name
129	!~~~> Local variables
130	REAL(kind=dp) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
131	REAL(kind=dp) :: Hmin, Hmax, Hstart, Hexit
132	REAL(kind=dp) :: Texit
133	INTEGER :: i, UplimTol, Max_no_steps
134	LOGICAL :: Autonomous, VectorTol
135	!~~~> Parameters
136	REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp, ONE = 1.0_dp
137	REAL(kind=dp), PARAMETER :: DeltaMin = 1.0E-5_dp
138
139	!~~~> Initialize statistics
140	Nfun = ISTATUS(ifun)
141	Njac = ISTATUS(ijac)
142	Nstp = ISTATUS(istp)
143	Nacc = ISTATUS(iacc)
144	Nrej = ISTATUS(irej)
145	Ndec = ISTATUS(idec)
146	Nsol = ISTATUS(isol)
147	Nsng = ISTATUS(isng)
148
149	!~~~> Autonomous or time dependent ODE. Default is time dependent.
150	Autonomous = .NOT.(ICNTRL(1) == 0)
151
152	!~~~> For Scalar tolerances (ICNTRL(2).NE.0) the code uses AbsTol(1) and RelTol(1)
153	! For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:NVAR) and RelTol(1:NVAR)
154	IF (ICNTRL(2) == 0) THEN
155	VectorTol = .TRUE.
156	UplimTol = NVAR
157	ELSE
158	VectorTol = .FALSE.
159	UplimTol = 1
160	END IF
161
162	!~~~> The particular Rosenbrock method chosen
163	IF (ICNTRL(3) == 0) THEN
164	Method = 4
165	ELSEIF ( (ICNTRL(3) >= 1).AND.(ICNTRL(3) <= 5) ) THEN
166	Method = ICNTRL(3)
167	ELSE
168	PRINT * , 'User-selected Rosenbrock method: ICNTRL(3)=', Method
169	!CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
170	IERRV(:) = -2
171	RETURN
172	END IF
173
174	!~~~> The maximum number of steps admitted
175	IF (ICNTRL(4) == 0) THEN
176	Max_no_steps = 100000 ! 200
177	ELSEIF (ICNTRL(4) > 0) THEN
178	Max_no_steps=ICNTRL(4)
179	ELSE
180	PRINT * ,'User-selected max no. of steps: ICNTRL(4)=',ICNTRL(4)
181	!CALL ros_ErrorMsg(-1,Tstart,ZERO,IERR)
182	IERRV(:) = -1
183	RETURN
184	END IF
185
186	!~~~> Unit roundoff (1+Roundoff>1)
187	Roundoff = WLAMCH('E')
188
189	!~~~> Lower bound on the step size: (positive value)
190	IF (RCNTRL(1) == ZERO) THEN
191	Hmin = ZERO
192	ELSEIF (RCNTRL(1) > ZERO) THEN
193	Hmin = RCNTRL(1)
194	ELSE
195	PRINT * , 'User-selected Hmin: RCNTRL(1)=', RCNTRL(1)
196	!CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
197	IERRV(:) = -3
198	RETURN
199	END IF
200	!~~~> Upper bound on the step size: (positive value)
201	IF (RCNTRL(2) == ZERO) THEN
202	Hmax = ABS(Tend-Tstart)
203	ELSEIF (RCNTRL(2) > ZERO) THEN
204	Hmax = MIN(ABS(RCNTRL(2)),ABS(Tend-Tstart))
205	ELSE
206	PRINT * , 'User-selected Hmax: RCNTRL(2)=', RCNTRL(2)
207	!CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
208	IERRV(:) = -3
209	RETURN
210	END IF
211	!~~~> Starting step size: (positive value)
212	IF (RCNTRL(3) == ZERO) THEN
213	Hstart = MAX(Hmin,DeltaMin)
214	ELSEIF (RCNTRL(3) > ZERO) THEN
215	Hstart = MIN(ABS(RCNTRL(3)),ABS(Tend-Tstart))
216	ELSE
217	PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3)
218	!CALL ros_ErrorMsg(-3,Tstart,ZERO,IERR)
219	IERRV(:) = -3
220	RETURN
221	END IF
222	!~~~> Step size can be changed s.t. FacMin < Hnew/Hexit < FacMax
223	IF (RCNTRL(4) == ZERO) THEN
224	FacMin = 0.2_dp
225	ELSEIF (RCNTRL(4) > ZERO) THEN
226	FacMin = RCNTRL(4)
227	ELSE
228	PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4)
229	!CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
230	IERRV(:) = -4
231	RETURN
232	END IF
233	IF (RCNTRL(5) == ZERO) THEN
234	FacMax = 6.0_dp
235	ELSEIF (RCNTRL(5) > ZERO) THEN
236	FacMax = RCNTRL(5)
237	ELSE
238	PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5)
239	!CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
240	IERRV(:) = -4
241	RETURN
242	END IF
243	!~~~> FacRej: Factor to decrease step after 2 succesive rejections
244	IF (RCNTRL(6) == ZERO) THEN
245	FacRej = 0.1_dp
246	ELSEIF (RCNTRL(6) > ZERO) THEN
247	FacRej = RCNTRL(6)
248	ELSE
249	PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6)
250	!CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
251	IERRV(:) = -4
252	RETURN
253	END IF
254	!~~~> FacSafe: Safety Factor in the computation of new step size
255	IF (RCNTRL(7) == ZERO) THEN
256	FacSafe = 0.9_dp
257	ELSEIF (RCNTRL(7) > ZERO) THEN
258	FacSafe = RCNTRL(7)
259	ELSE
260	PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7)
261	!CALL ros_ErrorMsg(-4,Tstart,ZERO,IERR)
262	IERRV(:) = -4
263	RETURN
264	END IF
265	!~~~> Check if tolerances are reasonable
266	DO i=1,UplimTol
267	IF ( (AbsTol(i) <= ZERO) .OR. (RelTol(i) <= 10.0_dp*Roundoff) &
268	.OR. (RelTol(i) >= 1.0_dp) ) THEN
269	PRINT * , ' AbsTol(',i,') = ',AbsTol(i)
270	PRINT * , ' RelTol(',i,') = ',RelTol(i)
271	!CALL ros_ErrorMsg(-5,Tstart,ZERO,IERR)
272	IERRV(:) = -5
273	RETURN
274	END IF
275	END DO
276
277
278	!~~~> Initialize the particular Rosenbrock method
279	SELECT CASE (Method)
280	CASE (1)
281	CALL Ros2(ros_S, ros_A, ros_C, ros_M, ros_E, &
282	ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
283	CASE (2)
284	CALL Ros3(ros_S, ros_A, ros_C, ros_M, ros_E, &
285	ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
286	CASE (3)
287	CALL Ros4(ros_S, ros_A, ros_C, ros_M, ros_E, &
288	ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
289	CASE (4)
290	CALL Rodas3(ros_S, ros_A, ros_C, ros_M, ros_E, &
291	ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
292	CASE (5)
293	CALL Rodas4(ros_S, ros_A, ros_C, ros_M, ros_E, &
294	ros_Alpha, ros_Gamma, ros_NewF, ros_ELO, ros_Name)
295	CASE DEFAULT
296	PRINT * , 'Unknown Rosenbrock method: ICNTRL(3)=', Method
297	!CALL ros_ErrorMsg(-2,Tstart,ZERO,IERR)
298	IERRV(:) = -2
299	RETURN
300	END SELECT
301
302	!~~~> CALL Rosenbrock method
303	CALL ros_Integrator(Y,Tstart,Tend,Texit, &
304	AbsTol, RelTol, &
305	! Rosenbrock method coefficients
306	ros_S, ros_M, ros_E, ros_A, ros_C, &
307	ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, &
308	! Integration parameters
309	Autonomous, VectorTol, Max_no_steps, &
310	Roundoff, Hmin, Hmax, Hstart, Hexit, &
311	FacMin, FacMax, FacRej, FacSafe, &
312	! Error indicator
313	IERR)
314
315
316	!~~~> Collect run statistics
317	ISTATUS(ifun) = Nfun
318	ISTATUS(ijac) = Njac
319	ISTATUS(istp) = Nstp
320	ISTATUS(iacc) = Nacc
321	ISTATUS(irej) = Nrej
322	ISTATUS(idec) = Ndec
323	ISTATUS(isol) = Nsol
324	ISTATUS(isng) = Nsng
325	!~~~> Last T and H
326	RSTATUS(itexit) = Texit
327	RSTATUS(ihexit) = Hexit
328
329	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
330	CONTAINS ! SUBROUTINES internal to Rosenbrock
331	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
332
333	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
334	SUBROUTINE ros_ErrorMsg(Code,T,H,IERR)
335	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
336	! Handles all error messages
337	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
338
339	REAL(kind=dp), INTENT(IN) :: T, H
340	INTEGER, INTENT(IN) :: Code
341	INTEGER, INTENT(INOUT) :: IERR
342
343	IERR = Code
344	PRINT * , &
345	'Forced exit from Rosenbrock due to the following error:'
346	IF ((Code>=-8).AND.(Code<=-1)) THEN
347	PRINT *, IERR_NAMES(Code)
348	ELSE
349	PRINT *, 'Unknown Error code: ', Code
350	ENDIF
351
352	PRINT *, "T=", T, "and H=", H
353
354	END SUBROUTINE ros_ErrorMsg
355
356	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
357	SUBROUTINE ros_Integrator (Y, Tstart, Tend, Tout, &
358	AbsTol, RelTol, &
359	!~~~> Rosenbrock method coefficients
360	ros_S, ros_M, ros_E, ros_A, ros_C, &
361	ros_Alpha, ros_Gamma, ros_ELO, ros_NewF, &
362	!~~~> Integration parameters
363	Autonomous, VectorTol, Max_no_steps, &
364	Roundoff, Hmin, Hmax, Hstart, Hexit, &
365	FacMin, FacMax, FacRej, FacSafe, &
366	!~~~> Error indicator
367	IERR )
368	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
369	! Template for the implementation of a generic Rosenbrock method
370	! defined by ros_S (no of stages)
371	! and its coefficients ros_{A,C,M,E,Alpha,Gamma}
372	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
373
374	USE messy_main_tools_kp4_compress, ONLY: kco_initialize, kco_compress, &
375	kco_finalize, cell_done
376	IMPLICIT NONE
377
378	!~~~> Input: the initial condition at Tstart; Output: the solution at T
379	REAL(kind=dp), INTENT(INOUT) :: Y(:,:)
380	!~~~> Input: integration interval
381	REAL(kind=dp), INTENT(IN) :: Tstart,Tend
382	!~~~> Output: time at which the solution is returned (T=Tend if success)
383	REAL(kind=dp), INTENT(OUT) :: Tout
384	!~~~> Input: tolerances
385	REAL(kind=dp), INTENT(IN) :: AbsTol(NVAR), RelTol(NVAR)
386	!~~~> Input: The Rosenbrock method parameters
387	INTEGER, INTENT(IN) :: ros_S
388	REAL(kind=dp), INTENT(IN) :: ros_M(ros_S), ros_E(ros_S), &
389	ros_Alpha(ros_S), ros_A(ros_S*(ros_S-1)/2), &
390	ros_Gamma(ros_S), ros_C(ros_S*(ros_S-1)/2), ros_ELO
391	LOGICAL, INTENT(IN) :: ros_NewF(ros_S)
392	!~~~> Input: integration parameters
393	LOGICAL, INTENT(IN) :: Autonomous, VectorTol
394	REAL(kind=dp), INTENT(IN) :: Hstart, Hmin, Hmax
395	INTEGER, INTENT(IN) :: Max_no_steps
396	REAL(kind=dp), INTENT(IN) :: Roundoff, FacMin, FacMax, FacRej, FacSafe
397	!~~~> Output: last accepted step
398	REAL(kind=dp), INTENT(OUT) :: Hexit
399	!~~~> Output: Error indicator
400	INTEGER, INTENT(OUT) :: IERR
401	! ~~~~ Local variables
402	REAL(kind=dp) :: Ynew(VL_DIM,NVAR), Fcn0(VL_DIM,NVAR), Fcn(VL_DIM,NVAR)
403	REAL(kind=dp) :: K(VL_DIM,NVAR*ros_S)
404	REAL(kind=dp) :: Jac0(VL_DIM,LU_NONZERO), Ghimj(VL_DIM,LU_NONZERO)
405	REAL(kind=dp) :: Yerr(VL_DIM,NVAR)
406
407	!~~~> Local parameters
408	INTEGER :: ioffset, i, j, istage
409	REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp, ONE = 1.0_dp
410	REAL(kind=dp), PARAMETER :: DeltaMin = 1.0E-5_dp
411	!~~~> Locally called functions
412	! REAL(kind=dp) WLAMCH
413	! EXTERNAL WLAMCH
414	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
415	real(kind=dp) :: deltat
416	!kk logical,parameter :: fixed_step=.true.
417
418	integer,save :: icount=0
419
420	! vector compression algorithm
421
422	REAL(kind=dp), dimension(VL_dim) :: T
423	REAL(kind=dp), dimension(VL_dim) :: H, Hnew, HC, HC2, HG, Fac, Tau, Err
424	LOGICAL, dimension(VL_dim) :: RejectLastH, RejectMoreH
425	LOGICAL, dimension(VL_dim) :: Accept_step
426	integer :: vl_save
427
428	REAL(kind=dp), dimension(size(FIX,1),size(FIX,2)) :: FIX_lo
429
430	!~~~> Initial preparations
431	T(1:VL) = Tstart
432	Hexit = 0.0_dp
433	H(1:VL) = MIN(Hstart,Hmax)
434	Tout = Tend
435	icount = icount+1
436	where (ABS(H(1:VL)) <= 10.0_dp*Roundoff) H(1:VL) = DeltaMin
437
438	deltat=Tend-Tstart
439
440	RejectLastH(1:VL)=.FALSE.
441	RejectMoreH(1:VL)=.FALSE.
442
443	!~~~> Time loop begins below
444
445	Nstp = 0
446	IERRV(1:VL) = 1 ! successful integration
447	vl_save = VL
448	FIX_lo(1:VL,:) = FIX(1:VL,:)
449	Kacc(1:VL) = 0
450	Krej(1:VL) = 0
451
452	! intialize compression algorithm
453
454	IF (.not. l_fixed_step) then
455	call kco_initialize (VL, y, IERRV, Kacc, Krej)
456	end if
457
458	TimeLoop: DO
459
460	IF ( Nstp > Max_no_steps ) THEN ! Too many steps
461	!CALL ros_ErrorMsg(-6,T(1),H(1),IERR)
462	IERRV(1:VL) = -6
463	! call abort() !Vorlaeufig
464	RETURN
465	END IF
466
467	Where ( (T(1:VL)+0.1_dp*H(1:VL)) == T(1:VL)) ! Step size too small
468	IERRV(1:VL) = -7
469	endwhere
470
471	!!$ if( any(error_flag)) then
472	!!$ CALL ros_ErrorMsg(-7,T(1),H(1),IERR)
473	!!$ call abort() !Vorlaeufig
474	!!$ RETURN
475	!!$ end if
476
477	!~~~> Limit H if necessary to avoid going beyond Tend
478	!kk Hexit = H
479	H(1:VL) = MIN(H(1:VL),ABS(Tend-T(1:VL)))
480
481	!~~~> Compute the function at current time
482	CALL Fun( Y, FIX_lo, RCONST, Fcn0)
483
484
485	!~~~> Compute the Jacobian at current time
486	CALL Jac_SP( Y, FIX_lo, RCONST, Jac0)
487
488
489	!~~~> Repeat step calculation until current step accepted
490
491	CALL ros_PrepareMatrix(H,ros_Gamma(1), Jac0,Ghimj)
492
493	!~~~> Compute the stages
494	Stage: DO istage = 1, ros_S
495
496	! Current istage offset. Current istage vector is K(ioffset+1:ioffset+NVAR)
497	ioffset = NVAR*(istage-1)
498
499	! For the 1st istage the function has been computed previously
500	IF ( istage == 1 ) THEN
501	Fcn(1:VL,1:NVAR) = Fcn0(1:VL,1:NVAR)
502	! istage>1 and a new function evaluation is needed at the current istage
503	ELSEIF ( ros_NewF(istage) ) THEN
504	select case (istage-1)
505	case (1)
506	Ynew(1:VL,1:NVAR) = Y(1:VL,1:NVAR)
507	Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + K(1:VL,NVAR(1- 1)+ 1:NVAR1) &
508	ros_A((istage-1)(istage-2)/2+ 1)
509	case (2)
510	Ynew(1:VL,1:NVAR) = Y(1:VL,1:NVAR)
511	Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + K(1:VL,NVAR(2- 1)+ 1:NVAR2) &
512	ros_A((istage-1)(istage-2)/2+ 2)
513	Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + K(1:VL,NVAR(2- 1)+ 1:NVAR2) &
514	ros_A((istage-1)(istage-2)/2+ 2)
515	case DEFAULT
516	Ynew(1:VL,1:NVAR) = Y(1:VL,1:NVAR)
517	DO j = 1, istage-1
518	Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + K(1:VL,NVAR(j- 1)+ 1:NVARj) &
519	ros_A((istage-1)(istage-2)/2+ j)
520	END DO
521	end select
522
523	CALL Fun( Ynew, FIX_lo, RCONST, Fcn)
524	END IF ! if istage == 1 elseif ros_NewF(istage)
525
526	select case (istage-1)
527	case (1)
528	K(1:VL,ioffset+ 1:ioffset+ NVAR) = Fcn(1:VL,1:NVAR)
529	j = 1
530	HC(1:VL) = ros_C((istage-1)*(istage-2)/2+ j)/(H(1:VL))
531	do i=1,NVAR
532	K(1:VL,ioffset+ i) = K(1:VL,ioffset+ i) + K(1:VL,NVAR(j- 1)+ i) HC(1:VL)
533	end do
534	case (2)
535	K(1:VL,ioffset+ 1:ioffset+ NVAR) = Fcn(1:VL,1:NVAR)
536	HC(1:VL) = ros_C((istage-1)*(istage-2)/2+ 1)/(H(1:VL))
537	HC2(1:VL) = ros_C((istage-1)*(istage-2)/2+ 2)/(H(1:VL))
538	do i=1,NVAR
539	K(1:VL,ioffset+ i) = K(1:VL,ioffset+ i) + K(1:VL,NVAR(1- 1)+ i) HC(1:VL)
540	K(1:VL,ioffset+ i) = K(1:VL,ioffset+ i) + K(1:VL,NVAR(2- 1)+ i) HC2(1:VL)
541	end do
542	case DEFAULT
543	K(1:VL,ioffset+ 1:ioffset+ NVAR) = Fcn(1:VL,1:NVAR)
544	DO j = 1, istage-1
545	HC(1:VL) = ros_C((istage-1)*(istage-2)/2+ j)/(H(1:VL))
546	do i=1,NVAR
547	K(1:VL,ioffset+ i) = K(1:VL,ioffset+ i) + K(1:VL,NVAR(j- 1)+ i) HC(1:VL)
548	end do
549	END DO
550	end select
551
552	CALL ros_Solve(Ghimj, K(:,ioffset+1:ioffset+NVAR))
553
554	END DO Stage
555
556
557	if(ros_S == 3) then
558	YNEW(1:VL,1:NVAR) = Y(1:VL,1:NVAR)
559	Yerr(1:VL,1:NVAR) = ZERO
560
561	!~~~> Compute the new solution
562	Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + K(1:VL,NVAR(1-1)+1:NVAR1) * ros_m(1)
563	Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + K(1:VL,NVAR(2-1)+1:NVAR2) * ros_m(2)
564	Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + K(1:VL,NVAR(3-1)+1:NVAR3) * ros_m(3)
565
566	!~~~> Compute the error estimation
567
568	Yerr(1:VL,1:NVAR) = Yerr(1:VL,1:NVAR) + K(1:VL,NVAR(1-1)+1:NVAR1) * ros_E(1)
569	Yerr(1:VL,1:NVAR) = Yerr(1:VL,1:NVAR) + K(1:VL,NVAR(2-1)+1:NVAR2) * ros_E(2)
570	Yerr(1:VL,1:NVAR) = Yerr(1:VL,1:NVAR) + K(1:VL,NVAR(3-1)+1:NVAR3) * ros_E(3)
571	else
572	YNEW(1:VL,1:NVAR) = Y(1:VL,1:NVAR)
573	Yerr(1:VL,1:NVAR) = ZERO
574	DO j=1,ros_S
575	!~~~> Compute the new solution
576	Ynew(1:VL,1:NVAR) = Ynew(1:VL,1:NVAR) + &
577	K(1:VL,NVAR(j-1)+1:NVARj) * ros_m(j)
578
579	!~~~> Compute the error estimation
580
581	Yerr(1:VL,1:NVAR) = Yerr(1:VL,1:NVAR) + &
582	K(1:VL,NVAR(j-1)+1:NVARj) * ros_E(j)
583	END DO
584	end if
585
586
587	!~~~> Check the error magnitude and adjust step size
588	Nstp = Nstp+1
589
590	if(l_fixed_step) then
591	Err = 0.0
592	H = t_steps(Nstp)*deltat
593	Hnew = t_steps(Nstp)*deltat
594	else
595	Err(1:VL) = ros_ErrorNorm ( Y, Ynew, Yerr, AbsTol, RelTol, VectorTol )
596	!~~~> New step size is bounded by FacMin <= Hnew/H <= FacMax
597	where(Err(1:VL) >= 10.0_dp*Roundoff)
598	Fac(1:VL) = MIN(FacMax,MAX(FacMin,FacSafe/Err(1:VL)**(ONE/ros_ELO)))
599	elsewhere
600	Fac(1:VL) = 1000. !Force large time step
601	endwhere
602	Hnew(1:VL) = H(1:VL)*Fac(1:VL)
603	end if
604
605	Accept_step(1:VL) = (Err(1:VL) <= ONE).OR.(H(1:VL) <= Hmin)
606	where (Accept_step(1:VL))
607	Kacc(1:VL) = Kacc(1:VL)+1
608	T(1:VL) = T(1:VL) + H(1:VL)
609	Hnew(1:VL) = MAX(Hmin,MIN(Hnew(1:VL),Hmax))
610	endwhere
611	do i=1,NVAR
612	where (Accept_step(1:VL))
613	Y(1:VL,i) = MAX(Ynew(1:VL,i),ZERO)
614	endwhere
615	end do
616	where (Accept_step(1:VL) .and. RejectLastH(1:VL))
617	Hnew(1:VL) = MIN(Hnew(1:VL),H(1:VL))
618	endwhere
619	where (Accept_step(1:VL))
620	RejectLastH(1:VL) = .FALSE.
621	RejectMoreH(1:VL) = .FALSE.
622	H(1:VL) = Hnew(1:VL)
623	endwhere
624	!kk ELSE !~~~> Reject step
625	where (.NOT. Accept_step(1:VL) .and. RejectMoreH(1:VL))
626	Hnew(1:VL) = H(1:VL)*FacRej
627	endwhere
628	where (.NOT. Accept_step(1:VL))
629	RejectMoreH(1:VL) = RejectLastH(1:VL)
630	RejectLastH(1:VL) = .TRUE.
631	H(1:VL) = Hnew(1:VL)
632	endwhere
633	where ( .NOT. Accept_step(1:VL) .and. Kacc(1:VL) >= 1)
634	Krej(1:VL) = Krej(1:VL)+1
635	endwhere
636
637
638	if(l_fixed_step) then
639	if(Nstp >= Nfsteps) EXIT
640	else
641	cell_done(1:VL) = ( .NOT. (T(1:VL)-Tend)+Roundoff <= ZERO ) .OR. &
642	(IERRV(1:VL) < 0) ! mz_pj_20070516
643	call kco_compress (VL, T, H, Hnew, IERRV, y, RCONST, FIX_lo, &
644	RejectLastH, RejectMoreH, Kacc, Krej)
645	if(VL == 0) EXIT
646	end if
647
648	!kk write(0,*) ' aaa ',icount,Nstp,count(cell_done),count(Accept_step),VL,minval(T)
649	END DO TimeLoop
650
651	!~~~> Succesful exit
652	IERR = 1 !~~~> The integration was successful
653
654	IF (.not. l_fixed_step) then
655	call kco_finalize (y, IERRV, Kacc, Krej)
656	end if
657
658	VL = VL_save
659
660	!kk write(0,*) 'bbb ',icount,Nstp,maxval(Kacc),maxval(Krej)
661	!kk if(icount >= 10) call abort()
662
663	END SUBROUTINE ros_Integrator
664
665
666	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
667	FUNCTION ros_ErrorNorm ( Y, Ynew, Yerr, &
668	AbsTol, RelTol, VectorTol )
669	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
670	!~~~> Computes the "scaled norm" of the error vector Yerr
671	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
672	IMPLICIT NONE
673
674	! Input arguments
675	REAL(kind=dp), INTENT(IN) :: Y(:,:), Ynew(:,:), &
676	Yerr(:,:), AbsTol(NVAR), RelTol(NVAR)
677	LOGICAL, INTENT(IN) :: VectorTol
678	! Local variables
679	REAL(kind=dp), dimension(VL) :: Err, Scale, Ymax
680	INTEGER :: i
681	REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp
682	REAL(kind=dp), dimension(VL) :: ros_ErrorNorm
683
684	Err = ZERO
685	DO i=1,NVAR
686	Ymax = MAX( ABS(Y(1:VL,i)),ABS(Ynew(1:VL,i)))
687	IF (VectorTol) THEN
688	Scale = AbsTol(i)+ RelTol(i)*Ymax
689	ELSE
690	Scale = AbsTol(1)+ RelTol(1)*Ymax
691	END IF
692	Err = Err+ (Yerr(1:VL,i)/Scale)**2
693	END DO
694	Err = SQRT(Err/NVAR)
695
696	ros_ErrorNorm = Err
697
698	END FUNCTION ros_ErrorNorm
699
700
701	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
702
703	SUBROUTINE ros_PrepareMatrix ( H, gam, Jac0, Ghimj )
704
705	! --- --- --- --- --- --- --- --- --- --- --- --- ---
706	! Prepares the LHS matrix for stage calculations
707	! 1. Construct Ghimj = 1/(H*ham) - Jac0
708	! "(Gamma H) Inverse Minus Jacobian"
709	! 2. Repeat LU decomposition of Ghimj until successful.
710	! -half the step size if LU decomposition fails and retry
711	! -exit after 5 consecutive fails
712	! --- --- --- --- --- --- --- --- --- --- --- --- ---
713	IMPLICIT NONE
714
715	!~~~> Input arguments
716	REAL(kind=dp), INTENT(IN) :: Jac0(:,:)
717	REAL(kind=dp), INTENT(IN) :: gam
718	!~~~> Output arguments
719	REAL(kind=dp), INTENT(OUT) :: Ghimj(:,:)
720	!~~~> Inout arguments
721	REAL(kind=dp), INTENT(INOUT) :: H(:) ! step size is decreased when LU fails
722	!~~~> Local variables
723	INTEGER :: i, ising, Nconsecutive
724	REAL(kind=dp), dimension(size(H)) :: ghinv
725	REAL(kind=dp), PARAMETER :: ONE = 1.0_dp, HALF = 0.5_dp
726
727	Nconsecutive = 0
728
729
730	!~~~> Construct Ghimj = 1/(H*gam) - Jac0
731	Ghimj(1:VL,1:LU_NONZERO) = -ONE * JAC0(1:VL,1:LU_NONZERO)
732	ghinv(1:VL) = ONE/(H(1:VL)*gam)
733	DO i=1,NVAR
734	Ghimj(1:VL,LU_DIAG(i)) = Ghimj(1:VL,LU_DIAG(i))+ghinv(1:VL)
735	END DO
736
737	!~~~> Compute LU decomposition
738	!CDIR noiexpand
739	CALL KppDecomp( Ghimj, ising )
740
741	Ndec = Ndec + VL
742
743	END SUBROUTINE ros_PrepareMatrix
744
745
746	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
747	SUBROUTINE ros_Solve( A, b )
748	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
749	! Template for the forward/backward substitution (using pre-computed LU decomposition)
750	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
751	IMPLICIT NONE
752	!~~~> Input variables
753	REAL(kind=dp), INTENT(IN) :: A(:,:)
754	!~~~> InOut variables
755	REAL(kind=dp), INTENT(INOUT) :: b(:,:)
756
757	CALL KppSolve( A, b )
758
759	Nsol = Nsol+VL
760
761	END SUBROUTINE ros_Solve
762
763
764
765	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
766	SUBROUTINE Ros2 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
767	ros_Gamma,ros_NewF,ros_ELO,ros_Name)
768	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
769	! --- AN L-STABLE METHOD, 2 stages, order 2
770	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
771
772	IMPLICIT NONE
773
774	INTEGER, PARAMETER :: S=2
775	INTEGER, INTENT(OUT) :: ros_S
776	REAL(kind=dp), DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
777	REAL(kind=dp), DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
778	REAL(kind=dp), INTENT(OUT) :: ros_ELO
779	LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
780	CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
781	DOUBLE PRECISION g
782
783	g = 1.0_dp + 1.0_dp/SQRT(2.0_dp)
784
785	!~~~> Name of the method
786	ros_Name = 'ROS-2'
787	!~~~> Number of stages
788	ros_S = S
789
790	!~~~> The coefficient matrices A and C are strictly lower triangular.
791	! The lower triangular (subdiagonal) elements are stored in row-wise order:
792	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
793	! The general mapping formula is:
794	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
795	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
796
797	ros_A(1) = (1.0_dp)/g
798	ros_C(1) = (-2.0_dp)/g
799	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
800	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
801	ros_NewF(1) = .TRUE.
802	ros_NewF(2) = .TRUE.
803	!~~~> M_i = Coefficients for new step solution
804	ros_M(1)= (3.0_dp)/(2.0_dp*g)
805	ros_M(2)= (1.0_dp)/(2.0_dp*g)
806	! E_i = Coefficients for error estimator
807	ros_E(1) = 1.0_dp/(2.0_dp*g)
808	ros_E(2) = 1.0_dp/(2.0_dp*g)
809	!~~~> ros_ELO = estimator of local order - the minimum between the
810	! main and the embedded scheme orders plus one
811	ros_ELO = 2.0_dp
812	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
813	ros_Alpha(1) = 0.0_dp
814	ros_Alpha(2) = 1.0_dp
815	!~~~> Gamma_i = \sum_j gamma_{i,j}
816	ros_Gamma(1) = g
817	ros_Gamma(2) =-g
818
819	END SUBROUTINE Ros2
820
821
822	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
823	SUBROUTINE Ros3 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
824	ros_Gamma,ros_NewF,ros_ELO,ros_Name)
825	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
826	! --- AN L-STABLE METHOD, 3 stages, order 3, 2 function evaluations
827	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
828
829	IMPLICIT NONE
830
831	INTEGER, PARAMETER :: S=3
832	INTEGER, INTENT(OUT) :: ros_S
833	REAL(kind=dp), DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
834	REAL(kind=dp), DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
835	REAL(kind=dp), INTENT(OUT) :: ros_ELO
836	LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
837	CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
838
839	!~~~> Name of the method
840	ros_Name = 'ROS-3'
841	!~~~> Number of stages
842	ros_S = S
843
844	!~~~> The coefficient matrices A and C are strictly lower triangular.
845	! The lower triangular (subdiagonal) elements are stored in row-wise order:
846	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
847	! The general mapping formula is:
848	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
849	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
850
851	ros_A(1)= 1.0_dp
852	ros_A(2)= 1.0_dp
853	ros_A(3)= 0.0_dp
854
855	ros_C(1) = -0.10156171083877702091975600115545E+01_dp
856	ros_C(2) = 0.40759956452537699824805835358067E+01_dp
857	ros_C(3) = 0.92076794298330791242156818474003E+01_dp
858	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
859	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
860	ros_NewF(1) = .TRUE.
861	ros_NewF(2) = .TRUE.
862	ros_NewF(3) = .FALSE.
863	!~~~> M_i = Coefficients for new step solution
864	ros_M(1) = 0.1E+01_dp
865	ros_M(2) = 0.61697947043828245592553615689730E+01_dp
866	ros_M(3) = -0.42772256543218573326238373806514E+00_dp
867	! E_i = Coefficients for error estimator
868	ros_E(1) = 0.5E+00_dp
869	ros_E(2) = -0.29079558716805469821718236208017E+01_dp
870	ros_E(3) = 0.22354069897811569627360909276199E+00_dp
871	!~~~> ros_ELO = estimator of local order - the minimum between the
872	! main and the embedded scheme orders plus 1
873	ros_ELO = 3.0_dp
874	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
875	ros_Alpha(1)= 0.0E+00_dp
876	ros_Alpha(2)= 0.43586652150845899941601945119356E+00_dp
877	ros_Alpha(3)= 0.43586652150845899941601945119356E+00_dp
878	!~~~> Gamma_i = \sum_j gamma_{i,j}
879	ros_Gamma(1)= 0.43586652150845899941601945119356E+00_dp
880	ros_Gamma(2)= 0.24291996454816804366592249683314E+00_dp
881	ros_Gamma(3)= 0.21851380027664058511513169485832E+01_dp
882
883	END SUBROUTINE Ros3
884
885	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
886
887
888	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
889	SUBROUTINE Ros4 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
890	ros_Gamma,ros_NewF,ros_ELO,ros_Name)
891	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
892	! L-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 4 STAGES
893	! L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
894	!
895	! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
896	! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
897	! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
898	! SPRINGER-VERLAG (1990)
899	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
900
901	IMPLICIT NONE
902
903	INTEGER, PARAMETER :: S=4
904	INTEGER, INTENT(OUT) :: ros_S
905	REAL(kind=dp), DIMENSION(4), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
906	REAL(kind=dp), DIMENSION(6), INTENT(OUT) :: ros_A, ros_C
907	REAL(kind=dp), INTENT(OUT) :: ros_ELO
908	LOGICAL, DIMENSION(4), INTENT(OUT) :: ros_NewF
909	CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
910	DOUBLE PRECISION g
911
912	!~~~> Name of the method
913	ros_Name = 'ROS-4'
914	!~~~> Number of stages
915	ros_S = S
916
917	!~~~> The coefficient matrices A and C are strictly lower triangular.
918	! The lower triangular (subdiagonal) elements are stored in row-wise order:
919	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
920	! The general mapping formula is:
921	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
922	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
923
924	ros_A(1) = 0.2000000000000000E+01_dp
925	ros_A(2) = 0.1867943637803922E+01_dp
926	ros_A(3) = 0.2344449711399156E+00_dp
927	ros_A(4) = ros_A(2)
928	ros_A(5) = ros_A(3)
929	ros_A(6) = 0.0_dp
930
931	ros_C(1) =-0.7137615036412310E+01_dp
932	ros_C(2) = 0.2580708087951457E+01_dp
933	ros_C(3) = 0.6515950076447975E+00_dp
934	ros_C(4) =-0.2137148994382534E+01_dp
935	ros_C(5) =-0.3214669691237626E+00_dp
936	ros_C(6) =-0.6949742501781779E+00_dp
937	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
938	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
939	ros_NewF(1) = .TRUE.
940	ros_NewF(2) = .TRUE.
941	ros_NewF(3) = .TRUE.
942	ros_NewF(4) = .FALSE.
943	!~~~> M_i = Coefficients for new step solution
944	ros_M(1) = 0.2255570073418735E+01_dp
945	ros_M(2) = 0.2870493262186792E+00_dp
946	ros_M(3) = 0.4353179431840180E+00_dp
947	ros_M(4) = 0.1093502252409163E+01_dp
948	!~~~> E_i = Coefficients for error estimator
949	ros_E(1) =-0.2815431932141155E+00_dp
950	ros_E(2) =-0.7276199124938920E-01_dp
951	ros_E(3) =-0.1082196201495311E+00_dp
952	ros_E(4) =-0.1093502252409163E+01_dp
953	!~~~> ros_ELO = estimator of local order - the minimum between the
954	! main and the embedded scheme orders plus 1
955	ros_ELO = 4.0_dp
956	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
957	ros_Alpha(1) = 0.0_dp
958	ros_Alpha(2) = 0.1145640000000000E+01_dp
959	ros_Alpha(3) = 0.6552168638155900E+00_dp
960	ros_Alpha(4) = ros_Alpha(3)
961	!~~~> Gamma_i = \sum_j gamma_{i,j}
962	ros_Gamma(1) = 0.5728200000000000E+00_dp
963	ros_Gamma(2) =-0.1769193891319233E+01_dp
964	ros_Gamma(3) = 0.7592633437920482E+00_dp
965	ros_Gamma(4) =-0.1049021087100450E+00_dp
966
967	END SUBROUTINE Ros4
968
969	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
970	SUBROUTINE Rodas3 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
971	ros_Gamma,ros_NewF,ros_ELO,ros_Name)
972	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
973	! --- A STIFFLY-STABLE METHOD, 4 stages, order 3
974	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
975
976	IMPLICIT NONE
977
978	INTEGER, PARAMETER :: S=4
979	INTEGER, INTENT(OUT) :: ros_S
980	REAL(kind=dp), DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
981	REAL(kind=dp), DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
982	REAL(kind=dp), INTENT(OUT) :: ros_ELO
983	LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
984	CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
985	DOUBLE PRECISION g
986
987	!~~~> Name of the method
988	ros_Name = 'RODAS-3'
989	!~~~> Number of stages
990	ros_S = S
991
992	!~~~> The coefficient matrices A and C are strictly lower triangular.
993	! The lower triangular (subdiagonal) elements are stored in row-wise order:
994	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
995	! The general mapping formula is:
996	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
997	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
998
999	ros_A(1) = 0.0E+00_dp
1000	ros_A(2) = 2.0E+00_dp
1001	ros_A(3) = 0.0E+00_dp
1002	ros_A(4) = 2.0E+00_dp
1003	ros_A(5) = 0.0E+00_dp
1004	ros_A(6) = 1.0E+00_dp
1005
1006	ros_C(1) = 4.0E+00_dp
1007	ros_C(2) = 1.0E+00_dp
1008	ros_C(3) =-1.0E+00_dp
1009	ros_C(4) = 1.0E+00_dp
1010	ros_C(5) =-1.0E+00_dp
1011	ros_C(6) =-(8.0E+00_dp/3.0E+00_dp)
1012
1013	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1014	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1015	ros_NewF(1) = .TRUE.
1016	ros_NewF(2) = .FALSE.
1017	ros_NewF(3) = .TRUE.
1018	ros_NewF(4) = .TRUE.
1019	!~~~> M_i = Coefficients for new step solution
1020	ros_M(1) = 2.0E+00_dp
1021	ros_M(2) = 0.0E+00_dp
1022	ros_M(3) = 1.0E+00_dp
1023	ros_M(4) = 1.0E+00_dp
1024	!~~~> E_i = Coefficients for error estimator
1025	ros_E(1) = 0.0E+00_dp
1026	ros_E(2) = 0.0E+00_dp
1027	ros_E(3) = 0.0E+00_dp
1028	ros_E(4) = 1.0E+00_dp
1029	!~~~> ros_ELO = estimator of local order - the minimum between the
1030	! main and the embedded scheme orders plus 1
1031	ros_ELO = 3.0E+00_dp
1032	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1033	ros_Alpha(1) = 0.0E+00_dp
1034	ros_Alpha(2) = 0.0E+00_dp
1035	ros_Alpha(3) = 1.0E+00_dp
1036	ros_Alpha(4) = 1.0E+00_dp
1037	!~~~> Gamma_i = \sum_j gamma_{i,j}
1038	ros_Gamma(1) = 0.5E+00_dp
1039	ros_Gamma(2) = 1.5E+00_dp
1040	ros_Gamma(3) = 0.0E+00_dp
1041	ros_Gamma(4) = 0.0E+00_dp
1042
1043	END SUBROUTINE Rodas3
1044
1045	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1046	SUBROUTINE Rodas4 (ros_S,ros_A,ros_C,ros_M,ros_E,ros_Alpha,&
1047	ros_Gamma,ros_NewF,ros_ELO,ros_Name)
1048	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1049	! STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4, WITH 6 STAGES
1050	!
1051	! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
1052	! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
1053	! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
1054	! SPRINGER-VERLAG (1996)
1055	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1056
1057	IMPLICIT NONE
1058
1059	INTEGER, PARAMETER :: S=6
1060	INTEGER, INTENT(OUT) :: ros_S
1061	REAL(kind=dp), DIMENSION(S), INTENT(OUT) :: ros_M,ros_E,ros_Alpha,ros_Gamma
1062	REAL(kind=dp), DIMENSION(S*(S-1)/2), INTENT(OUT) :: ros_A, ros_C
1063	REAL(kind=dp), INTENT(OUT) :: ros_ELO
1064	LOGICAL, DIMENSION(S), INTENT(OUT) :: ros_NewF
1065	CHARACTER(LEN=12), INTENT(OUT) :: ros_Name
1066	DOUBLE PRECISION g
1067
1068	!~~~> Name of the method
1069	ros_Name = 'RODAS-4'
1070	!~~~> Number of stages
1071	ros_S = 6
1072
1073	!~~~> Y_stage_i ~ Y( T + H*Alpha_i )
1074	ros_Alpha(1) = 0.000_dp
1075	ros_Alpha(2) = 0.386_dp
1076	ros_Alpha(3) = 0.210_dp
1077	ros_Alpha(4) = 0.630_dp
1078	ros_Alpha(5) = 1.000_dp
1079	ros_Alpha(6) = 1.000_dp
1080
1081	!~~~> Gamma_i = \sum_j gamma_{i,j}
1082	ros_Gamma(1) = 0.2500000000000000E+00_dp
1083	ros_Gamma(2) =-0.1043000000000000E+00_dp
1084	ros_Gamma(3) = 0.1035000000000000E+00_dp
1085	ros_Gamma(4) =-0.3620000000000023E-01_dp
1086	ros_Gamma(5) = 0.0_dp
1087	ros_Gamma(6) = 0.0_dp
1088
1089	!~~~> The coefficient matrices A and C are strictly lower triangular.
1090	! The lower triangular (subdiagonal) elements are stored in row-wise order:
1091	! A(2,1) = ros_A(1), A(3,1)=ros_A(2), A(3,2)=ros_A(3), etc.
1092	! The general mapping formula is: A(i,j) = ros_A( (i-1)*(i-2)/2 + j )
1093	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j )
1094
1095	ros_A(1) = 0.1544000000000000E+01_dp
1096	ros_A(2) = 0.9466785280815826E+00_dp
1097	ros_A(3) = 0.2557011698983284E+00_dp
1098	ros_A(4) = 0.3314825187068521E+01_dp
1099	ros_A(5) = 0.2896124015972201E+01_dp
1100	ros_A(6) = 0.9986419139977817E+00_dp
1101	ros_A(7) = 0.1221224509226641E+01_dp
1102	ros_A(8) = 0.6019134481288629E+01_dp
1103	ros_A(9) = 0.1253708332932087E+02_dp
1104	ros_A(10) =-0.6878860361058950E+00_dp
1105	ros_A(11) = ros_A(7)
1106	ros_A(12) = ros_A(8)
1107	ros_A(13) = ros_A(9)
1108	ros_A(14) = ros_A(10)
1109	ros_A(15) = 1.0E+00_dp
1110
1111	ros_C(1) =-0.5668800000000000E+01_dp
1112	ros_C(2) =-0.2430093356833875E+01_dp
1113	ros_C(3) =-0.2063599157091915E+00_dp
1114	ros_C(4) =-0.1073529058151375E+00_dp
1115	ros_C(5) =-0.9594562251023355E+01_dp
1116	ros_C(6) =-0.2047028614809616E+02_dp
1117	ros_C(7) = 0.7496443313967647E+01_dp
1118	ros_C(8) =-0.1024680431464352E+02_dp
1119	ros_C(9) =-0.3399990352819905E+02_dp
1120	ros_C(10) = 0.1170890893206160E+02_dp
1121	ros_C(11) = 0.8083246795921522E+01_dp
1122	ros_C(12) =-0.7981132988064893E+01_dp
1123	ros_C(13) =-0.3152159432874371E+02_dp
1124	ros_C(14) = 0.1631930543123136E+02_dp
1125	ros_C(15) =-0.6058818238834054E+01_dp
1126
1127	!~~~> M_i = Coefficients for new step solution
1128	ros_M(1) = ros_A(7)
1129	ros_M(2) = ros_A(8)
1130	ros_M(3) = ros_A(9)
1131	ros_M(4) = ros_A(10)
1132	ros_M(5) = 1.0E+00_dp
1133	ros_M(6) = 1.0E+00_dp
1134
1135	!~~~> E_i = Coefficients for error estimator
1136	ros_E(1) = 0.0E+00_dp
1137	ros_E(2) = 0.0E+00_dp
1138	ros_E(3) = 0.0E+00_dp
1139	ros_E(4) = 0.0E+00_dp
1140	ros_E(5) = 0.0E+00_dp
1141	ros_E(6) = 1.0E+00_dp
1142
1143	!~~~> Does the stage i require a new function evaluation (ros_NewF(i)=TRUE)
1144	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i)=FALSE)
1145	ros_NewF(1) = .TRUE.
1146	ros_NewF(2) = .TRUE.
1147	ros_NewF(3) = .TRUE.
1148	ros_NewF(4) = .TRUE.
1149	ros_NewF(5) = .TRUE.
1150	ros_NewF(6) = .TRUE.
1151
1152	!~~~> ros_ELO = estimator of local order - the minimum between the
1153	! main and the embedded scheme orders plus 1
1154	ros_ELO = 4.0_dp
1155
1156	END SUBROUTINE Rodas4
1157
1158	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1159	! End of the set of internal Rosenbrock subroutines
1160	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1161	END SUBROUTINE Rosenbrock

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |