Home

Context Navigation

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

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

Rev	Line
[2696]	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 |