Home

Context Navigation

source: palm/trunk/SOURCE/chem_gasphase_mod.f90 @ 2723

Last change on this file since 2723 was 2718, checked in by maronga, 7 years ago
deleting of deprecated files; headers updated where needed
File size: 79.2 KB

Rev	Line
[2657]	1	MODULE chem_gasphase_mod
[2668]	2	#if defined( __chem )
[2657]	3	!------------------------------------------------------------------------------!
	4	!
	5	! ****Module chem_gasphase_mod is automatically generated by kpp4palm ****
	6	!
	7	! *******Please do NOT change this Code *******
	8	!
	9	!------------------------------------------------------------------------------!
[2696]	10	! This file is part of the PALM model system.
[2657]	11	!
	12	! PALM is free software: you can redistribute it and/or modify it under the
	13	! terms of the GNU General Public License as published by the Free Software
[2696]	14	! Foundation, either version 3 of the License, or (at your option) any later
[2657]	15	! version.
	16	!
[2696]	17	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
[2657]	18	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	19	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	20	!
	21	! You should have received a copy of the GNU General Public License along with
[2696]	22	! PALM. If not, see <http://www.gnu.org/licenses/>.
[2657]	23	!
[2718]	24	! Copyright 1997-2018 Leibniz Universitaet Hannover
[2657]	25	!--------------------------------------------------------------------------------!
	26	!
	27	! Current revisions:
	28	! ------------------
	29	!
[2678]	30	!
[2657]	31	! Former revisions:
	32	! -----------------
	33	! $Id: module_header 2460 2017-09-13 14:47:48Z forkel $
	34	!
	35	!
	36	! Variables for photolyis added
	37	!
	38	!
	39	!
	40	!
	41	!
	42	! Nov. 2016: Intial version (Klaus Ketelsen)
	43	!
	44	!------------------------------------------------------------------------------!
	45	!
	46
	47
	48	! Set kpp Double Precision to PALM Default Precision
	49
	50	USE kinds, ONLY: dp=>wp
	51
	52	USE pegrid, ONLY: myid,threads_per_task
	53
	54	IMPLICIT NONE
	55	PRIVATE
	56	!SAVE ! NOTE: OCCURS AGAIN IN AUTOMATICALLY GENERATED CODE ...
	57
	58	! PUBLIC :: IERR_NAMES
	59
	60	! PUBLIC :: SPC_NAMES,EQN_NAMES,EQN_TAGS,REQ_HET,REQ_AEROSOL,REQ_PHOTRAT &
	61	! ,REQ_MCFCT,IP_MAX,jname
	62
	63	PUBLIC :: eqn_names, phot_names,spc_names
	64	PUBLIC :: nmaxfixsteps
	65	PUBLIC :: atol,rtol
	66	PUBLIC :: nspec,nreact
	67	PUBLIC :: temp
	68	PUBLIC :: phot
	69	PUBLIC :: rconst
	70	PUBLIC :: nvar
	71	PUBLIC :: nphot
	72
	73	PUBLIC :: initialize,integrate,update_rconst
	74	PUBLIC :: chem_gasphase_integrate
	75	PUBLIC :: initialize_kpp_ctrl
	76
	77	! END OF MODULE HEADER TEMPLATE
	78
	79	! Variables used for vector mode
	80
	81	LOGICAL, PARAMETER :: l_vector = .FALSE.
	82	INTEGER, PARAMETER :: i_lu_di = 0
	83	INTEGER, PARAMETER :: vl_dim = 1
	84	INTEGER :: vl
	85
	86	INTEGER :: vl_glo
	87	INTEGER :: is,ie
	88
	89	INTEGER, DIMENSION(vl_dim) :: kacc,krej
	90	INTEGER, DIMENSION(vl_dim) :: ierrv
	91	LOGICAL :: data_loaded = .false.
	92	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	93	!
	94	! Parameter Module File
	95	!
	96	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
	97	! (http://www.cs.vt.edu/~asandu/Software/KPP)
	98	! KPP is distributed under GPL,the general public licence
	99	! (http://www.gnu.org/copyleft/gpl.html)
	100	! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
	101	! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
	102	! With important contributions from:
	103	! M. Damian,Villanova University,USA
	104	! R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
	105	!
	106	! File : chem_gasphase_mod_Parameters.f90
[2678]	107	! Time : Fri Dec 8 11:54:15 2017
	108	! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
[2657]	109	! Equation file : chem_gasphase_mod.kpp
	110	! Output root filename : chem_gasphase_mod
	111	!
	112	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	113
	114
	115
	116
	117
	118
	119	! NSPEC - Number of chemical species
[2678]	120	INTEGER,PARAMETER :: nspec = 2
[2657]	121	! NVAR - Number of Variable species
[2678]	122	INTEGER,PARAMETER :: nvar = 2
[2657]	123	! NVARACT - Number of Active species
[2678]	124	INTEGER,PARAMETER :: nvaract = 2
[2657]	125	! NFIX - Number of Fixed species
	126	INTEGER,PARAMETER :: nfix = 1
	127	! NREACT - Number of reactions
[2678]	128	INTEGER,PARAMETER :: nreact = 2
[2657]	129	! NVARST - Starting of variables in conc. vect.
	130	INTEGER,PARAMETER :: nvarst = 1
	131	! NFIXST - Starting of fixed in conc. vect.
[2678]	132	INTEGER,PARAMETER :: nfixst = 3
[2657]	133	! NONZERO - Number of nonzero entries in Jacobian
[2678]	134	INTEGER,PARAMETER :: nonzero = 2
[2657]	135	! LU_NONZERO - Number of nonzero entries in LU factoriz. of Jacobian
[2678]	136	INTEGER,PARAMETER :: lu_nonzero = 2
[2657]	137	! CNVAR - (NVAR+1) Number of elements in compressed row format
[2678]	138	INTEGER,PARAMETER :: cnvar = 3
[2657]	139	! CNEQN - (NREACT+1) Number stoicm elements in compressed col format
[2678]	140	INTEGER,PARAMETER :: cneqn = 3
[2657]	141	! NHESS - Length of Sparse Hessian
	142	INTEGER,PARAMETER :: nhess = 1
	143	! NLOOKAT - Number of species to look at
	144	INTEGER,PARAMETER :: nlookat = 0
	145	! NMONITOR - Number of species to monitor
	146	INTEGER,PARAMETER :: nmonitor = 0
	147	! NMASS - Number of atoms to check mass balance
	148	INTEGER,PARAMETER :: nmass = 1
	149
	150	! Index declaration for variable species in C and VAR
	151	! VAR(ind_spc) = C(ind_spc)
	152
[2678]	153	INTEGER,PARAMETER,PUBLIC :: ind_pm10 = 1
	154	INTEGER,PARAMETER,PUBLIC :: ind_pm25 = 2
[2657]	155
	156	! Index declaration for fixed species in C
	157	! C(ind_spc)
	158
	159
	160	! Index declaration for fixed species in FIX
	161	! FIX(indf_spc) = C(ind_spc) = C(NVAR+indf_spc)
	162
	163
	164	! NJVRP - Length of sparse Jacobian JVRP
[2678]	165	INTEGER,PARAMETER :: njvrp = 2
[2657]	166
	167	! NSTOICM - Length of Sparse Stoichiometric Matrix
	168	INTEGER,PARAMETER :: nstoicm = 1
	169
	170
	171	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	172	!
	173	! Global Data Module File
	174	!
	175	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
	176	! (http://www.cs.vt.edu/~asandu/Software/KPP)
	177	! KPP is distributed under GPL,the general public licence
	178	! (http://www.gnu.org/copyleft/gpl.html)
	179	! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
	180	! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
	181	! With important contributions from:
	182	! M. Damian,Villanova University,USA
	183	! R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
	184	!
	185	! File : chem_gasphase_mod_Global.f90
[2678]	186	! Time : Fri Dec 8 11:54:15 2017
	187	! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
[2657]	188	! Equation file : chem_gasphase_mod.kpp
	189	! Output root filename : chem_gasphase_mod
	190	!
	191	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	192
	193
	194
	195
	196
	197
	198	! Declaration of global variables
	199
	200	! C - Concentration of all species
	201	REAL(kind=dp):: c(nspec)
	202	! VAR - Concentrations of variable species (global)
	203	REAL(kind=dp):: var(nvar)
	204	! FIX - Concentrations of fixed species (global)
	205	REAL(kind=dp):: fix(nfix)
	206	! VAR,FIX are chunks of array C
	207	equivalence( c(1),var(1))
	208	! RCONST - Rate constants (global)
	209	REAL(kind=dp):: rconst(nreact)
	210	! TIME - Current integration time
	211	REAL(kind=dp):: time
	212	! SUN - Sunlight intensity between [0,1]
	213	REAL(kind=dp):: sun
	214	! TEMP - Temperature
	215	REAL(dp),dimension(:),allocatable :: temp
	216	! RTOLS - (scalar) Relative tolerance
	217	REAL(kind=dp):: rtols
	218	! TSTART - Integration start time
	219	REAL(kind=dp):: tstart
	220	! TEND - Integration end time
	221	REAL(kind=dp):: tend
	222	! DT - Integration step
	223	REAL(kind=dp):: dt
	224	! ATOL - Absolute tolerance
	225	REAL(kind=dp):: atol(nvar)
	226	! RTOL - Relative tolerance
	227	REAL(kind=dp):: rtol(nvar)
	228	! STEPMIN - Lower bound for integration step
	229	REAL(kind=dp):: stepmin
	230	! STEPMAX - Upper bound for integration step
	231	REAL(kind=dp):: stepmax
	232	! CFACTOR - Conversion factor for concentration units
	233	REAL(kind=dp):: cfactor
	234	! DDMTYPE - DDM sensitivity w.r.t.: 0=init.val.,1=params
	235	INTEGER :: ddmtype
	236
	237	! INLINED global variable declarations
	238
	239	! declaration of global variable declarations for photolysis will come from
	240
	241	! INLINED global variable declarations
	242
	243
	244
	245	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	246	!
	247	! Sparse Jacobian Data Structures File
	248	!
	249	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
	250	! (http://www.cs.vt.edu/~asandu/Software/KPP)
	251	! KPP is distributed under GPL,the general public licence
	252	! (http://www.gnu.org/copyleft/gpl.html)
	253	! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
	254	! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
	255	! With important contributions from:
	256	! M. Damian,Villanova University,USA
	257	! R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
	258	!
	259	! File : chem_gasphase_mod_JacobianSP.f90
[2678]	260	! Time : Fri Dec 8 11:54:15 2017
	261	! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
[2657]	262	! Equation file : chem_gasphase_mod.kpp
	263	! Output root filename : chem_gasphase_mod
	264	!
	265	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	266
	267
	268
	269
	270
	271
	272	! Sparse Jacobian Data
	273
	274
[2678]	275	INTEGER,PARAMETER,DIMENSION(2):: lu_irow = (/ &
	276	1, 2 /)
	277
	278	INTEGER,PARAMETER,DIMENSION(2):: lu_icol = (/ &
	279	1, 2 /)
	280
	281	INTEGER,PARAMETER,DIMENSION(3):: lu_crow = (/ &
[2657]	282	1, 2, 3 /)
	283
[2678]	284	INTEGER,PARAMETER,DIMENSION(3):: lu_diag = (/ &
[2657]	285	1, 2, 3 /)
	286
	287
	288
	289	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	290	!
	291	! Utility Data Module File
	292	!
	293	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
	294	! (http://www.cs.vt.edu/~asandu/Software/KPP)
	295	! KPP is distributed under GPL,the general public licence
	296	! (http://www.gnu.org/copyleft/gpl.html)
	297	! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
	298	! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
	299	! With important contributions from:
	300	! M. Damian,Villanova University,USA
	301	! R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
	302	!
	303	! File : chem_gasphase_mod_Monitor.f90
[2678]	304	! Time : Fri Dec 8 11:54:15 2017
	305	! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
[2657]	306	! Equation file : chem_gasphase_mod.kpp
	307	! Output root filename : chem_gasphase_mod
	308	!
	309	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	310
	311
	312
	313
	314
[2678]	315	CHARACTER(len=15),PARAMETER,DIMENSION(2):: spc_names = (/ &
	316	'PM10 ','PM25 ' /)
[2657]	317
	318	INTEGER,DIMENSION(1):: lookat
	319	INTEGER,DIMENSION(1):: monitor
	320	CHARACTER(len=15),DIMENSION(1):: smass
[2678]	321	CHARACTER(len=100),PARAMETER,DIMENSION(2):: eqn_names = (/ &
[2657]	322	'PM10 --> PM10 ',&
	323	'PM25 --> PM25 ' /)
	324
	325	! INLINED global variables
	326
	327	! inline f90_data: declaration of global variables for photolysis
	328	! REAL(kind=dp):: phot(nphot)must eventually be moved to global later for
	329	INTEGER,PARAMETER :: nphot = 1
	330	! phot photolysis frequencies
	331	REAL(kind=dp):: phot(nphot)
	332
	333	INTEGER,PARAMETER,PUBLIC :: j_no2 = 1
	334
	335	CHARACTER(len=15),PARAMETER,DIMENSION(nphot):: phot_names = (/ &
	336	'J_NO2 '/)
	337
	338	! End INLINED global variables
	339
	340
	341	! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
	342
	343	! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
	344
	345	! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
	346
	347	! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
	348
	349
	350	! variable definations from individual module headers
	351
	352	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	353	!
	354	! Initialization File
	355	!
	356	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
	357	! (http://www.cs.vt.edu/~asandu/Software/KPP)
	358	! KPP is distributed under GPL,the general public licence
	359	! (http://www.gnu.org/copyleft/gpl.html)
	360	! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
	361	! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
	362	! With important contributions from:
	363	! M. Damian,Villanova University,USA
	364	! R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
	365	!
	366	! File : chem_gasphase_mod_Initialize.f90
[2678]	367	! Time : Fri Dec 8 11:54:15 2017
	368	! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
[2657]	369	! Equation file : chem_gasphase_mod.kpp
	370	! Output root filename : chem_gasphase_mod
	371	!
	372	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	373
	374
	375
	376
	377
	378	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	379	!
	380	! Numerical Integrator (Time-Stepping) File
	381	!
	382	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
	383	! (http://www.cs.vt.edu/~asandu/Software/KPP)
	384	! KPP is distributed under GPL,the general public licence
	385	! (http://www.gnu.org/copyleft/gpl.html)
	386	! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
	387	! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
	388	! With important contributions from:
	389	! M. Damian,Villanova University,USA
	390	! R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
	391	!
	392	! File : chem_gasphase_mod_Integrator.f90
[2678]	393	! Time : Fri Dec 8 11:54:15 2017
	394	! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
[2657]	395	! Equation file : chem_gasphase_mod.kpp
	396	! Output root filename : chem_gasphase_mod
	397	!
	398	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	399
	400
	401
	402	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	403	!
	404	! INTEGRATE - Integrator routine
	405	! Arguments :
	406	! TIN - Start Time for Integration
	407	! TOUT - End Time for Integration
	408	!
	409	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	410
	411	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
	412	! Rosenbrock - Implementation of several Rosenbrock methods: !
	413	! *Ros2 !
	414	! *Ros3 !
	415	! *Ros4 !
	416	! *Rodas3 !
	417	! *Rodas4 !
	418	! By default the code employs the KPP sparse linear algebra routines !
	419	! Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK) !
	420	! !
	421	! (C) Adrian Sandu,August 2004 !
	422	! Virginia Polytechnic Institute and State University !
	423	! Contact: sandu@cs.vt.edu !
	424	! Revised by Philipp Miehe and Adrian Sandu,May 2006 ! !
	425	! This implementation is part of KPP - the Kinetic PreProcessor !
	426	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
	427
	428
	429	SAVE
	430
	431	!~~~> statistics on the work performed by the rosenbrock method
	432	INTEGER,PARAMETER :: nfun=1,njac=2,nstp=3,nacc=4,&
	433	nrej=5,ndec=6,nsol=7,nsng=8,&
	434	ntexit=1,nhexit=2,nhnew = 3
	435
	436	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	437	!
	438	! Linear Algebra Data and Routines File
	439	!
	440	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
	441	! (http://www.cs.vt.edu/~asandu/Software/KPP)
	442	! KPP is distributed under GPL,the general public licence
	443	! (http://www.gnu.org/copyleft/gpl.html)
	444	! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
	445	! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
	446	! With important contributions from:
	447	! M. Damian,Villanova University,USA
	448	! R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
	449	!
	450	! File : chem_gasphase_mod_LinearAlgebra.f90
[2678]	451	! Time : Fri Dec 8 11:54:15 2017
	452	! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
[2657]	453	! Equation file : chem_gasphase_mod.kpp
	454	! Output root filename : chem_gasphase_mod
	455	!
	456	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	457
	458
	459
	460
	461
	462
	463	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	464	!
	465	! The ODE Jacobian of Chemical Model File
	466	!
	467	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
	468	! (http://www.cs.vt.edu/~asandu/Software/KPP)
	469	! KPP is distributed under GPL,the general public licence
	470	! (http://www.gnu.org/copyleft/gpl.html)
	471	! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
	472	! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
	473	! With important contributions from:
	474	! M. Damian,Villanova University,USA
	475	! R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
	476	!
	477	! File : chem_gasphase_mod_Jacobian.f90
[2678]	478	! Time : Fri Dec 8 11:54:15 2017
	479	! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
[2657]	480	! Equation file : chem_gasphase_mod.kpp
	481	! Output root filename : chem_gasphase_mod
	482	!
	483	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	484
	485
	486
	487
	488
	489
	490	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	491	!
	492	! The ODE Function of Chemical Model File
	493	!
	494	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
	495	! (http://www.cs.vt.edu/~asandu/Software/KPP)
	496	! KPP is distributed under GPL,the general public licence
	497	! (http://www.gnu.org/copyleft/gpl.html)
	498	! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
	499	! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
	500	! With important contributions from:
	501	! M. Damian,Villanova University,USA
	502	! R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
	503	!
	504	! File : chem_gasphase_mod_Function.f90
[2678]	505	! Time : Fri Dec 8 11:54:15 2017
	506	! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
[2657]	507	! Equation file : chem_gasphase_mod.kpp
	508	! Output root filename : chem_gasphase_mod
	509	!
	510	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	511
	512
	513
	514
	515
	516	! A - Rate for each equation
	517	REAL(kind=dp):: a(nreact)
	518
	519	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	520	!
	521	! The Reaction Rates File
	522	!
	523	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
	524	! (http://www.cs.vt.edu/~asandu/Software/KPP)
	525	! KPP is distributed under GPL,the general public licence
	526	! (http://www.gnu.org/copyleft/gpl.html)
	527	! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
	528	! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
	529	! With important contributions from:
	530	! M. Damian,Villanova University,USA
	531	! R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
	532	!
	533	! File : chem_gasphase_mod_Rates.f90
[2678]	534	! Time : Fri Dec 8 11:54:15 2017
	535	! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
[2657]	536	! Equation file : chem_gasphase_mod.kpp
	537	! Output root filename : chem_gasphase_mod
	538	!
	539	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	540
	541
	542
	543
	544
	545	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	546	!
	547	! Auxiliary Routines File
	548	!
	549	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
	550	! (http://www.cs.vt.edu/~asandu/Software/KPP)
	551	! KPP is distributed under GPL,the general public licence
	552	! (http://www.gnu.org/copyleft/gpl.html)
	553	! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
	554	! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
	555	! With important contributions from:
	556	! M. Damian,Villanova University,USA
	557	! R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
	558	!
	559	! File : chem_gasphase_mod_Util.f90
[2678]	560	! Time : Fri Dec 8 11:54:15 2017
	561	! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
[2657]	562	! Equation file : chem_gasphase_mod.kpp
	563	! Output root filename : chem_gasphase_mod
	564	!
	565	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	566
	567
	568
	569
	570
	571
	572	! header MODULE initialize_kpp_ctrl_template
	573
	574	! notes:
	575	! - l_vector is automatically defined by kp4
	576	! - vl_dim is automatically defined by kp4
	577	! - i_lu_di is automatically defined by kp4
	578	! - wanted is automatically defined by xmecca
	579	! - icntrl rcntrl are automatically defined by kpp
	580	! - "USE messy_main_tools" is in MODULE_header of messy_mecca_kpp.f90
	581	! - SAVE will be automatically added by kp4
	582
	583	!SAVE
	584
	585	! for fixed time step control
	586	! ... max. number of fixed time steps (sum must be 1)
	587	INTEGER,PARAMETER :: nmaxfixsteps = 50
	588	! ... switch for fixed time stepping
	589	LOGICAL,PUBLIC :: l_fixed_step = .false.
	590	INTEGER,PUBLIC :: nfsteps = 1
	591	! ... number of kpp control PARAMETERs
	592	INTEGER,PARAMETER,PUBLIC :: nkppctrl = 20
	593	!
	594	INTEGER, DIMENSION(nkppctrl),PUBLIC :: icntrl = 0
	595	REAL(dp),DIMENSION(nkppctrl),PUBLIC :: rcntrl = 0.0_dp
	596	REAL(dp),DIMENSION(nmaxfixsteps),PUBLIC :: t_steps = 0.0_dp
	597
	598	! END header MODULE initialize_kpp_ctrl_template
	599
	600
	601	! Interface Block
	602
	603	INTERFACE initialize
	604	MODULE PROCEDURE initialize
	605	END INTERFACE initialize
	606
	607	INTERFACE integrate
	608	MODULE PROCEDURE integrate
	609	END INTERFACE integrate
	610
	611	INTERFACE fun
	612	MODULE PROCEDURE fun
	613	END INTERFACE fun
	614
	615	INTERFACE kppsolve
	616	MODULE PROCEDURE kppsolve
	617	END INTERFACE kppsolve
	618
	619	INTERFACE kppdecomp
	620	MODULE PROCEDURE kppdecomp
	621	END INTERFACE kppdecomp
	622
	623	INTERFACE wlamch
	624	MODULE PROCEDURE wlamch
	625	END INTERFACE wlamch
	626
	627	INTERFACE wlamch_add
	628	MODULE PROCEDURE wlamch_add
	629	END INTERFACE wlamch_add
	630
	631	INTERFACE jac_sp
	632	MODULE PROCEDURE jac_sp
	633	END INTERFACE jac_sp
	634
	635	INTERFACE k_arr
	636	MODULE PROCEDURE k_arr
	637	END INTERFACE k_arr
	638
	639	INTERFACE update_rconst
	640	MODULE PROCEDURE update_rconst
	641	END INTERFACE update_rconst
	642
	643	INTERFACE arr2
	644	MODULE PROCEDURE arr2
	645	END INTERFACE arr2
	646
	647	INTERFACE initialize_kpp_ctrl
	648	MODULE PROCEDURE initialize_kpp_ctrl
	649	END INTERFACE initialize_kpp_ctrl
	650
	651	INTERFACE error_output
	652	MODULE PROCEDURE error_output
	653	END INTERFACE error_output
	654
	655	!interface not working INTERFACE wcopy
	656	!interface not working MODULE PROCEDURE wcopy
	657	!interface not working END INTERFACE wcopy
	658
	659	INTERFACE wscal
	660	MODULE PROCEDURE wscal
	661	END INTERFACE wscal
	662
	663	!interface not working INTERFACE waxpy
	664	!interface not working MODULE PROCEDURE waxpy
	665	!interface not working END INTERFACE waxpy
	666
	667	INTERFACE rosenbrock
	668	MODULE PROCEDURE rosenbrock
	669	END INTERFACE rosenbrock
	670
	671	INTERFACE funtemplate
	672	MODULE PROCEDURE funtemplate
	673	END INTERFACE funtemplate
	674
	675	INTERFACE jactemplate
	676	MODULE PROCEDURE jactemplate
	677	END INTERFACE jactemplate
	678
	679	INTERFACE chem_gasphase_integrate
	680	MODULE PROCEDURE chem_gasphase_integrate
	681	END INTERFACE chem_gasphase_integrate
	682
	683	INTERFACE fill_temp
	684	MODULE PROCEDURE fill_temp
	685	END INTERFACE fill_temp
	686	PUBLIC fill_temp
	687
	688
	689	CONTAINS
	690
	691	SUBROUTINE initialize()
	692
	693
	694
	695	INTEGER :: i
	696	REAL(kind=dp):: x
	697
	698	cfactor = 1.000000e+00_dp
	699
	700	x = (0.)*cfactor
	701	DO i = 1,nvar
	702	var(i) = x
	703	ENDDO
	704
	705	x = (0.)*cfactor
	706	DO i = 1,nfix
	707	fix(i) = x
	708	ENDDO
	709
	710	! constant rate coefficients
	711	! END constant rate coefficients
	712
	713	! INLINED initializations
	714
	715	! End INLINED initializations
	716
	717
	718	END SUBROUTINE initialize
	719
	720	SUBROUTINE integrate( tin,tout,&
	721	icntrl_u,rcntrl_u,istatus_u,rstatus_u,ierr_u)
	722
	723
	724	REAL(kind=dp),INTENT(in):: tin ! start time
	725	REAL(kind=dp),INTENT(in):: tout ! END time
	726	! OPTIONAL input PARAMETERs and statistics
	727	INTEGER, INTENT(in), OPTIONAL :: icntrl_u(20)
	728	REAL(kind=dp),INTENT(in), OPTIONAL :: rcntrl_u(20)
	729	INTEGER, INTENT(out),OPTIONAL :: istatus_u(20)
	730	REAL(kind=dp),INTENT(out),OPTIONAL :: rstatus_u(20)
	731	INTEGER, INTENT(out),OPTIONAL :: ierr_u
	732
	733	REAL(kind=dp):: rcntrl(20),rstatus(20)
	734	INTEGER :: icntrl(20),istatus(20),ierr
	735
	736	INTEGER,SAVE :: ntotal = 0
	737
	738	icntrl(:) = 0
	739	rcntrl(:) = 0.0_dp
	740	istatus(:) = 0
	741	rstatus(:) = 0.0_dp
	742
	743	!~~~> fine-tune the integrator:
	744	icntrl(1) = 0 ! 0 - non- autonomous,1 - autonomous
	745	icntrl(2) = 0 ! 0 - vector tolerances,1 - scalars
	746
	747	! IF OPTIONAL PARAMETERs are given,and IF they are >0,
	748	! THEN they overwrite default settings.
	749	IF (present(icntrl_u))THEN
	750	where(icntrl_u(:)> 0)icntrl(:) = icntrl_u(:)
	751	ENDIF
	752	IF (present(rcntrl_u))THEN
	753	where(rcntrl_u(:)> 0)rcntrl(:) = rcntrl_u(:)
	754	ENDIF
	755
	756
	757	CALL rosenbrock(nvar,var,tin,tout, &
	758	atol,rtol, &
	759	rcntrl,icntrl,rstatus,istatus,ierr)
	760
	761	!~~~> debug option: show no of steps
	762	! ntotal = ntotal + istatus(nstp)
	763	! PRINT*,'NSTEPS=',ISTATUS(Nstp),' (',Ntotal,')',' O3=',VAR(ind_O3)
	764
	765	stepmin = rstatus(nhexit)
	766	! IF OPTIONAL PARAMETERs are given for output they
	767	! are updated with the RETURN information
	768	IF (present(istatus_u))istatus_u(:) = istatus(:)
	769	IF (present(rstatus_u))rstatus_u(:) = rstatus(:)
	770	IF (present(ierr_u)) ierr_u = ierr
	771
	772	END SUBROUTINE integrate
	773
	774	SUBROUTINE fun(v,f,rct,vdot)
	775
	776	! V - Concentrations of variable species (local)
	777	REAL(kind=dp):: v(nvar)
	778	! F - Concentrations of fixed species (local)
	779	REAL(kind=dp):: f(nfix)
	780	! RCT - Rate constants (local)
	781	REAL(kind=dp):: rct(nreact)
	782	! Vdot - Time derivative of variable species concentrations
	783	REAL(kind=dp):: vdot(nvar)
	784
	785
	786	! Computation of equation rates
	787
	788	! Aggregate function
	789	vdot(1) = 0
	790	vdot(2) = 0
	791
	792	END SUBROUTINE fun
	793
	794	SUBROUTINE kppsolve(jvs,x)
	795
	796	! JVS - sparse Jacobian of variables
	797	REAL(kind=dp):: jvs(lu_nonzero)
	798	! X - Vector for variables
	799	REAL(kind=dp):: x(nvar)
	800
	801	x(2) = x(2)/ jvs(2)
	802	x(1) = x(1)/ jvs(1)
	803
	804	END SUBROUTINE kppsolve
	805
	806	SUBROUTINE kppdecomp( jvs,ier)
	807	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	808	! Sparse LU factorization
	809	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	810
	811
	812	INTEGER :: ier
	813	REAL(kind=dp):: jvs(lu_nonzero),w(nvar),a
	814	INTEGER :: k,kk,j,jj
	815
	816	a = 0. ! mz_rs_20050606
	817	ier = 0
	818	DO k=1,nvar
	819	! mz_rs_20050606: don't check if real value == 0
	820	! IF(jvs( lu_diag(k)).eq. 0.)THEN
	821	IF(abs(jvs(lu_diag(k)))< tiny(a))THEN
	822	ier = k
	823	RETURN
	824	ENDIF
	825	DO kk = lu_crow(k),lu_crow(k+ 1)- 1
	826	w( lu_icol(kk)) = jvs(kk)
	827	ENDDO
	828	DO kk = lu_crow(k),lu_diag(k)- 1
	829	j = lu_icol(kk)
	830	a = - w(j)/ jvs( lu_diag(j))
	831	w(j) = - a
	832	DO jj = lu_diag(j)+ 1,lu_crow(j+ 1)- 1
	833	w( lu_icol(jj)) = w( lu_icol(jj))+ a*jvs(jj)
	834	ENDDO
	835	ENDDO
	836	DO kk = lu_crow(k),lu_crow(k+ 1)- 1
	837	jvs(kk) = w( lu_icol(kk))
	838	ENDDO
	839	ENDDO
	840
	841	END SUBROUTINE kppdecomp
	842
	843	REAL(kind=dp)FUNCTION wlamch( c)
	844	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	845	! returns epsilon machine
	846	! after LAPACK
	847	! replace this by the function from the optimized LAPACK implementation:
	848	! CALL SLAMCH('E') or CALL DLAMCH('E')
	849	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	850	! USE chem_gasphase_mod_Precision
	851
	852	CHARACTER :: c
	853	INTEGER :: i
	854	REAL(kind=dp),SAVE :: eps
	855	REAL(kind=dp) :: suma
	856	REAL(kind=dp),PARAMETER :: one=1.0_dp, half=0.5_dp
	857	LOGICAL,SAVE :: first=.true.
	858
	859	IF (first)THEN
	860	first = .false.
	861	eps = half**(16)
	862	DO i = 17,80
	863	eps = eps*half
	864	CALL wlamch_add(one,eps,suma)
	865	IF (suma.le.one)goto 10
	866	ENDDO
	867	PRINT*,'ERROR IN WLAMCH. EPS < ',Eps
	868	RETURN
	869	10 eps = eps*2
	870	i = i- 1
	871	ENDIF
	872
	873	wlamch = eps
	874
	875	END FUNCTION wlamch
	876
	877	SUBROUTINE wlamch_add( a,b,suma)
	878	! USE chem_gasphase_mod_Precision
	879
	880	REAL(kind=dp)a,b,suma
	881	suma = a + b
	882
	883	END SUBROUTINE wlamch_add
	884
	885	SUBROUTINE jac_sp(v,f,rct,jvs)
	886
	887	! V - Concentrations of variable species (local)
	888	REAL(kind=dp):: v(nvar)
	889	! F - Concentrations of fixed species (local)
	890	REAL(kind=dp):: f(nfix)
	891	! RCT - Rate constants (local)
	892	REAL(kind=dp):: rct(nreact)
	893	! JVS - sparse Jacobian of variables
	894	REAL(kind=dp):: jvs(lu_nonzero)
	895
	896
	897	! Local variables
	898	! B - Temporary array
[2678]	899	REAL(kind=dp):: b(2)
[2657]	900
	901	! B(1) = dA(1)/dV(1)
	902	b(1) = rct(1)
	903	! B(2) = dA(2)/dV(2)
	904	b(2) = rct(2)
	905
	906	! Construct the Jacobian terms from B's
	907	! JVS(1) = Jac_FULL(1,1)
	908	jvs(1) = 0
	909	! JVS(2) = Jac_FULL(2,2)
	910	jvs(2) = 0
	911
	912	END SUBROUTINE jac_sp
	913
	914	elemental REAL(kind=dp)FUNCTION k_arr (k_298,tdep,temp)
	915	! arrhenius FUNCTION
	916
	917	REAL, INTENT(in):: k_298 ! k at t = 298.15k
	918	REAL, INTENT(in):: tdep ! temperature dependence
	919	REAL(kind=dp),INTENT(in):: temp ! temperature
	920
	921	intrinsic exp
	922
	923	k_arr = k_298 exp(tdep(1._dp/temp- 3.3540e-3_dp))! 1/298.15=3.3540e-3
	924
	925	END FUNCTION k_arr
	926
	927	SUBROUTINE update_rconst()
	928	INTEGER :: j,k
	929
	930	k = is
	931
	932	! Begin INLINED RCONST
	933
	934
	935	! End INLINED RCONST
	936
	937	rconst(1) = (1.0_dp)
	938	rconst(2) = (1.0_dp)
	939
	940	END SUBROUTINE update_rconst
	941
	942	REAL(kind=dp)FUNCTION arr2( a0,b0,temp)
	943	REAL(kind=dp):: temp
	944	REAL(kind=dp):: a0,b0
	945	arr2 = a0 *exp( - b0 / temp)
	946	END FUNCTION arr2
	947
	948	SUBROUTINE initialize_kpp_ctrl(status,iou,modstr)
	949
	950
	951	! i/o
	952	INTEGER, INTENT(out):: status
	953	INTEGER, INTENT(in) :: iou ! LOGICAL i/o unit
	954	CHARACTER(len=*),INTENT(in) :: modstr ! read <modstr>.nml
	955
	956	! local
	957	REAL(dp):: tsum
	958	INTEGER :: i
	959
	960	! check fixed time steps
	961	tsum = 0.0_dp
	962	DO i=1,nmaxfixsteps
	963	IF (t_steps(i)< tiny(0.0_dp))exit
	964	tsum = tsum + t_steps(i)
	965	ENDDO
	966
	967	nfsteps = i- 1
	968
	969	l_fixed_step = (nfsteps > 0).and.((tsum - 1.0)< tiny(0.0_dp))
	970
	971	IF (l_vector)THEN
	972	WRITE(,) ' MODE : VECTOR (LENGTH=',VL_DIM,')'
	973	ELSE
	974	WRITE(,) ' MODE : SCALAR'
	975	ENDIF
	976	!
	977	WRITE(,) ' DE-INDEXING MODE :',I_LU_DI
	978	!
	979	WRITE(,) ' ICNTRL : ',icntrl
	980	WRITE(,) ' RCNTRL : ',rcntrl
	981	!
	982	! note: this is only meaningful for vectorized (kp4)rosenbrock- methods
	983	IF (l_vector)THEN
	984	IF (l_fixed_step)THEN
	985	WRITE(,) ' TIME STEPS : FIXED (',t_steps(1:nfsteps),')'
	986	ELSE
	987	WRITE(,) ' TIME STEPS : AUTOMATIC'
	988	ENDIF
	989	ELSE
	990	WRITE(,) ' TIME STEPS : AUTOMATIC '//&
	991	&'(t_steps (CTRL_KPP) ignored in SCALAR MODE)'
	992	ENDIF
	993	! mz_pj_20070531-
	994
	995	status = 0
	996
	997
	998	END SUBROUTINE initialize_kpp_ctrl
	999
	1000	SUBROUTINE error_output(c,ierr,pe)
	1001
	1002
	1003	INTEGER,INTENT(in):: ierr
	1004	INTEGER,INTENT(in):: pe
	1005	REAL(dp),DIMENSION(:),INTENT(in):: c
	1006
	1007	write(6,*) 'ERROR in chem_gasphase_mod ',ierr,C(1)
	1008
	1009
	1010	END SUBROUTINE error_output
	1011
	1012	SUBROUTINE wcopy(n,x,incx,y,incy)
	1013	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	1014	! copies a vector,x,to a vector,y: y <- x
	1015	! only for incX=incY=1
	1016	! after BLAS
	1017	! replace this by the function from the optimized BLAS implementation:
	1018	! CALL SCOPY(N,X,1,Y,1) or CALL DCOPY(N,X,1,Y,1)
	1019	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	1020	! USE chem_gasphase_mod_Precision
	1021
	1022	INTEGER :: i,incx,incy,m,mp1,n
	1023	REAL(kind=dp):: x(n),y(n)
	1024
	1025	IF (n.le.0)RETURN
	1026
	1027	m = mod(n,8)
	1028	IF( m .ne. 0)THEN
	1029	DO i = 1,m
	1030	y(i) = x(i)
	1031	ENDDO
	1032	IF( n .lt. 8)RETURN
	1033	ENDIF
	1034	mp1 = m+ 1
	1035	DO i = mp1,n,8
	1036	y(i) = x(i)
	1037	y(i + 1) = x(i + 1)
	1038	y(i + 2) = x(i + 2)
	1039	y(i + 3) = x(i + 3)
	1040	y(i + 4) = x(i + 4)
	1041	y(i + 5) = x(i + 5)
	1042	y(i + 6) = x(i + 6)
	1043	y(i + 7) = x(i + 7)
	1044	ENDDO
	1045
	1046	END SUBROUTINE wcopy
	1047
	1048	SUBROUTINE wscal(n,alpha,x,incx)
	1049	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	1050	! constant times a vector: x(1:N) <- Alpha*x(1:N)
	1051	! only for incX=incY=1
	1052	! after BLAS
	1053	! replace this by the function from the optimized BLAS implementation:
	1054	! CALL SSCAL(N,Alpha,X,1) or CALL DSCAL(N,Alpha,X,1)
	1055	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	1056
	1057	INTEGER :: i,incx,m,mp1,n
	1058	REAL(kind=dp) :: x(n),alpha
	1059	REAL(kind=dp),PARAMETER :: zero=0.0_dp, one=1.0_dp
	1060
	1061	IF (alpha .eq. one)RETURN
	1062	IF (n .le. 0)RETURN
	1063
	1064	m = mod(n,5)
	1065	IF( m .ne. 0)THEN
	1066	IF (alpha .eq. (- one))THEN
	1067	DO i = 1,m
	1068	x(i) = - x(i)
	1069	ENDDO
	1070	ELSEIF (alpha .eq. zero)THEN
	1071	DO i = 1,m
	1072	x(i) = zero
	1073	ENDDO
	1074	ELSE
	1075	DO i = 1,m
	1076	x(i) = alpha*x(i)
	1077	ENDDO
	1078	ENDIF
	1079	IF( n .lt. 5)RETURN
	1080	ENDIF
	1081	mp1 = m + 1
	1082	IF (alpha .eq. (- one))THEN
	1083	DO i = mp1,n,5
	1084	x(i) = - x(i)
	1085	x(i + 1) = - x(i + 1)
	1086	x(i + 2) = - x(i + 2)
	1087	x(i + 3) = - x(i + 3)
	1088	x(i + 4) = - x(i + 4)
	1089	ENDDO
	1090	ELSEIF (alpha .eq. zero)THEN
	1091	DO i = mp1,n,5
	1092	x(i) = zero
	1093	x(i + 1) = zero
	1094	x(i + 2) = zero
	1095	x(i + 3) = zero
	1096	x(i + 4) = zero
	1097	ENDDO
	1098	ELSE
	1099	DO i = mp1,n,5
	1100	x(i) = alpha*x(i)
	1101	x(i + 1) = alpha*x(i + 1)
	1102	x(i + 2) = alpha*x(i + 2)
	1103	x(i + 3) = alpha*x(i + 3)
	1104	x(i + 4) = alpha*x(i + 4)
	1105	ENDDO
	1106	ENDIF
	1107
	1108	END SUBROUTINE wscal
	1109
	1110	SUBROUTINE waxpy(n,alpha,x,incx,y,incy)
	1111	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	1112	! constant times a vector plus a vector: y <- y + Alpha*x
	1113	! only for incX=incY=1
	1114	! after BLAS
	1115	! replace this by the function from the optimized BLAS implementation:
	1116	! CALL SAXPY(N,Alpha,X,1,Y,1) or CALL DAXPY(N,Alpha,X,1,Y,1)
	1117	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	1118
	1119	INTEGER :: i,incx,incy,m,mp1,n
	1120	REAL(kind=dp):: x(n),y(n),alpha
	1121	REAL(kind=dp),PARAMETER :: zero = 0.0_dp
	1122
	1123	IF (alpha .eq. zero)RETURN
	1124	IF (n .le. 0)RETURN
	1125
	1126	m = mod(n,4)
	1127	IF( m .ne. 0)THEN
	1128	DO i = 1,m
	1129	y(i) = y(i)+ alpha*x(i)
	1130	ENDDO
	1131	IF( n .lt. 4)RETURN
	1132	ENDIF
	1133	mp1 = m + 1
	1134	DO i = mp1,n,4
	1135	y(i) = y(i)+ alpha*x(i)
	1136	y(i + 1) = y(i + 1)+ alpha*x(i + 1)
	1137	y(i + 2) = y(i + 2)+ alpha*x(i + 2)
	1138	y(i + 3) = y(i + 3)+ alpha*x(i + 3)
	1139	ENDDO
	1140
	1141	END SUBROUTINE waxpy
	1142
	1143	SUBROUTINE rosenbrock(n,y,tstart,tend,&
	1144	abstol,reltol, &
	1145	rcntrl,icntrl,rstatus,istatus,ierr)
	1146	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1147	!
	1148	! Solves the system y'=F(t,y) using a Rosenbrock method defined by:
	1149	!
	1150	! G = 1/(H*gamma(1)) - Jac(t0,Y0)
	1151	! T_i = t0 + Alpha(i)*H
	1152	! Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
	1153	! G K_i = Fun( T_i,Y_i)+ \sum_{j=1}^S C(i,j)/H K_j +
	1154	! gamma(i)*dF/dT(t0,Y0)
	1155	! Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
	1156	!
	1157	! For details on Rosenbrock methods and their implementation consult:
	1158	! E. Hairer and G. Wanner
	1159	! "Solving ODEs II. Stiff and differential-algebraic problems".
	1160	! Springer series in computational mathematics,Springer-Verlag,1996.
	1161	! The codes contained in the book inspired this implementation.
	1162	!
	1163	! (C) Adrian Sandu,August 2004
	1164	! Virginia Polytechnic Institute and State University
	1165	! Contact: sandu@cs.vt.edu
	1166	! Revised by Philipp Miehe and Adrian Sandu,May 2006
	1167	! This implementation is part of KPP - the Kinetic PreProcessor
	1168	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1169	!
	1170	!~~~> input arguments:
	1171	!
	1172	!- y(n) = vector of initial conditions (at t=tstart)
	1173	!- [tstart,tend] = time range of integration
	1174	! (if Tstart>Tend the integration is performed backwards in time)
	1175	!- reltol,abstol = user precribed accuracy
	1176	!- SUBROUTINE fun( t,y,ydot) = ode FUNCTION,
	1177	! returns Ydot = Y' = F(T,Y)
	1178	!- SUBROUTINE jac( t,y,jcb) = jacobian of the ode FUNCTION,
	1179	! returns Jcb = dFun/dY
	1180	!- icntrl(1:20) = INTEGER inputs PARAMETERs
	1181	!- rcntrl(1:20) = REAL inputs PARAMETERs
	1182	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1183	!
	1184	!~~~> output arguments:
	1185	!
	1186	!- y(n) - > vector of final states (at t- >tEND)
	1187	!- istatus(1:20) - > INTEGER output PARAMETERs
	1188	!- rstatus(1:20) - > REAL output PARAMETERs
	1189	!- ierr - > job status upon RETURN
	1190	! success (positive value) or
	1191	! failure (negative value)
	1192	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1193	!
	1194	!~~~> input PARAMETERs:
	1195	!
	1196	! Note: For input parameters equal to zero the default values of the
	1197	! corresponding variables are used.
	1198	!
	1199	! ICNTRL(1) = 1: F = F(y) Independent of T (AUTONOMOUS)
	1200	! = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
	1201	!
	1202	! ICNTRL(2) = 0: AbsTol,RelTol are N-dimensional vectors
	1203	! = 1: AbsTol,RelTol are scalars
	1204	!
	1205	! ICNTRL(3) -> selection of a particular Rosenbrock method
	1206	! = 0 : Rodas3 (default)
	1207	! = 1 : Ros2
	1208	! = 2 : Ros3
	1209	! = 3 : Ros4
	1210	! = 4 : Rodas3
	1211	! = 5 : Rodas4
	1212	!
	1213	! ICNTRL(4) -> maximum number of integration steps
	1214	! For ICNTRL(4) =0) the default value of 100000 is used
	1215	!
	1216	! RCNTRL(1) -> Hmin,lower bound for the integration step size
	1217	! It is strongly recommended to keep Hmin = ZERO
	1218	! RCNTRL(2) -> Hmax,upper bound for the integration step size
	1219	! RCNTRL(3) -> Hstart,starting value for the integration step size
	1220	!
	1221	! RCNTRL(4) -> FacMin,lower bound on step decrease factor (default=0.2)
	1222	! RCNTRL(5) -> FacMax,upper bound on step increase factor (default=6)
	1223	! RCNTRL(6) -> FacRej,step decrease factor after multiple rejections
	1224	! (default=0.1)
	1225	! RCNTRL(7) -> FacSafe,by which the new step is slightly smaller
	1226	! than the predicted value (default=0.9)
	1227	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1228	!
	1229	!
	1230	! OUTPUT ARGUMENTS:
	1231	! -----------------
	1232	!
	1233	! T -> T value for which the solution has been computed
	1234	! (after successful return T=Tend).
	1235	!
	1236	! Y(N) -> Numerical solution at T
	1237	!
	1238	! IDID -> Reports on successfulness upon return:
	1239	! = 1 for success
	1240	! < 0 for error (value equals error code)
	1241	!
	1242	! ISTATUS(1) -> No. of function calls
	1243	! ISTATUS(2) -> No. of jacobian calls
	1244	! ISTATUS(3) -> No. of steps
	1245	! ISTATUS(4) -> No. of accepted steps
	1246	! ISTATUS(5) -> No. of rejected steps (except at very beginning)
	1247	! ISTATUS(6) -> No. of LU decompositions
	1248	! ISTATUS(7) -> No. of forward/backward substitutions
	1249	! ISTATUS(8) -> No. of singular matrix decompositions
	1250	!
	1251	! RSTATUS(1) -> Texit,the time corresponding to the
	1252	! computed Y upon return
	1253	! RSTATUS(2) -> Hexit,last accepted step before exit
	1254	! RSTATUS(3) -> Hnew,last predicted step (not yet taken)
	1255	! For multiple restarts,use Hnew as Hstart
	1256	! in the subsequent run
	1257	!
	1258	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1259
	1260
	1261	!~~~> arguments
	1262	INTEGER, INTENT(in) :: n
	1263	REAL(kind=dp),INTENT(inout):: y(n)
	1264	REAL(kind=dp),INTENT(in) :: tstart,tend
	1265	REAL(kind=dp),INTENT(in) :: abstol(n),reltol(n)
	1266	INTEGER, INTENT(in) :: icntrl(20)
	1267	REAL(kind=dp),INTENT(in) :: rcntrl(20)
	1268	INTEGER, INTENT(inout):: istatus(20)
	1269	REAL(kind=dp),INTENT(inout):: rstatus(20)
	1270	INTEGER,INTENT(out) :: ierr
	1271	!~~~> PARAMETERs of the rosenbrock method,up to 6 stages
	1272	INTEGER :: ros_s,rosmethod
	1273	INTEGER,PARAMETER :: rs2=1,rs3=2,rs4=3,rd3=4,rd4=5,rg3=6
	1274	REAL(kind=dp):: ros_a(15),ros_c(15),ros_m(6),ros_e(6),&
	1275	ros_alpha(6),ros_gamma(6),ros_elo
	1276	LOGICAL :: ros_newf(6)
	1277	CHARACTER(len=12):: ros_name
	1278	!~~~> local variables
	1279	REAL(kind=dp):: roundoff,facmin,facmax,facrej,facsafe
	1280	REAL(kind=dp):: hmin,hmax,hstart
	1281	REAL(kind=dp):: texit
	1282	INTEGER :: i,uplimtol,max_no_steps
	1283	LOGICAL :: autonomous,vectortol
	1284	!~~~> PARAMETERs
	1285	REAL(kind=dp),PARAMETER :: zero = 0.0_dp, one = 1.0_dp
	1286	REAL(kind=dp),PARAMETER :: deltamin = 1.0e-5_dp
	1287
	1288	!~~~> initialize statistics
	1289	istatus(1:8) = 0
	1290	rstatus(1:3) = zero
	1291
	1292	!~~~> autonomous or time dependent ode. default is time dependent.
	1293	autonomous = .not.(icntrl(1) == 0)
	1294
	1295	!~~~> for scalar tolerances (icntrl(2).ne.0) the code uses abstol(1)and reltol(1)
	1296	! For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:N) and RelTol(1:N)
	1297	IF (icntrl(2) == 0)THEN
	1298	vectortol = .true.
	1299	uplimtol = n
	1300	ELSE
	1301	vectortol = .false.
	1302	uplimtol = 1
	1303	ENDIF
	1304
	1305	!~~~> initialize the particular rosenbrock method selected
	1306	select CASE (icntrl(3))
	1307	CASE (1)
	1308	CALL ros2
	1309	CASE (2)
	1310	CALL ros3
	1311	CASE (3)
	1312	CALL ros4
	1313	CASE (0,4)
	1314	CALL rodas3
	1315	CASE (5)
	1316	CALL rodas4
	1317	CASE (6)
	1318	CALL rang3
	1319	CASE default
	1320	PRINT *,'Unknown Rosenbrock method: ICNTRL(3) =',ICNTRL(3)
	1321	CALL ros_errormsg(- 2,tstart,zero,ierr)
	1322	RETURN
	1323	END select
	1324
	1325	!~~~> the maximum number of steps admitted
	1326	IF (icntrl(4) == 0)THEN
	1327	max_no_steps = 200000
	1328	ELSEIF (icntrl(4)> 0)THEN
	1329	max_no_steps=icntrl(4)
	1330	ELSE
	1331	PRINT *,'User-selected max no. of steps: ICNTRL(4) =',ICNTRL(4)
	1332	CALL ros_errormsg(- 1,tstart,zero,ierr)
	1333	RETURN
	1334	ENDIF
	1335
	1336	!~~~> unit roundoff (1+ roundoff>1)
	1337	Roundoff = WLAMCH('E')
	1338
	1339	!~~~> lower bound on the step size: (positive value)
	1340	IF (rcntrl(1) == zero)THEN
	1341	hmin = zero
	1342	ELSEIF (rcntrl(1)> zero)THEN
	1343	hmin = rcntrl(1)
	1344	ELSE
	1345	PRINT *,'User-selected Hmin: RCNTRL(1) =',RCNTRL(1)
	1346	CALL ros_errormsg(- 3,tstart,zero,ierr)
	1347	RETURN
	1348	ENDIF
	1349	!~~~> upper bound on the step size: (positive value)
	1350	IF (rcntrl(2) == zero)THEN
	1351	hmax = abs(tend-tstart)
	1352	ELSEIF (rcntrl(2)> zero)THEN
	1353	hmax = min(abs(rcntrl(2)),abs(tend-tstart))
	1354	ELSE
	1355	PRINT *,'User-selected Hmax: RCNTRL(2) =',RCNTRL(2)
	1356	CALL ros_errormsg(- 3,tstart,zero,ierr)
	1357	RETURN
	1358	ENDIF
	1359	!~~~> starting step size: (positive value)
	1360	IF (rcntrl(3) == zero)THEN
	1361	hstart = max(hmin,deltamin)
	1362	ELSEIF (rcntrl(3)> zero)THEN
	1363	hstart = min(abs(rcntrl(3)),abs(tend-tstart))
	1364	ELSE
	1365	PRINT *,'User-selected Hstart: RCNTRL(3) =',RCNTRL(3)
	1366	CALL ros_errormsg(- 3,tstart,zero,ierr)
	1367	RETURN
	1368	ENDIF
	1369	!~~~> step size can be changed s.t. facmin < hnew/hold < facmax
	1370	IF (rcntrl(4) == zero)THEN
	1371	facmin = 0.2_dp
	1372	ELSEIF (rcntrl(4)> zero)THEN
	1373	facmin = rcntrl(4)
	1374	ELSE
	1375	PRINT *,'User-selected FacMin: RCNTRL(4) =',RCNTRL(4)
	1376	CALL ros_errormsg(- 4,tstart,zero,ierr)
	1377	RETURN
	1378	ENDIF
	1379	IF (rcntrl(5) == zero)THEN
	1380	facmax = 6.0_dp
	1381	ELSEIF (rcntrl(5)> zero)THEN
	1382	facmax = rcntrl(5)
	1383	ELSE
	1384	PRINT *,'User-selected FacMax: RCNTRL(5) =',RCNTRL(5)
	1385	CALL ros_errormsg(- 4,tstart,zero,ierr)
	1386	RETURN
	1387	ENDIF
	1388	!~~~> facrej: factor to decrease step after 2 succesive rejections
	1389	IF (rcntrl(6) == zero)THEN
	1390	facrej = 0.1_dp
	1391	ELSEIF (rcntrl(6)> zero)THEN
	1392	facrej = rcntrl(6)
	1393	ELSE
	1394	PRINT *,'User-selected FacRej: RCNTRL(6) =',RCNTRL(6)
	1395	CALL ros_errormsg(- 4,tstart,zero,ierr)
	1396	RETURN
	1397	ENDIF
	1398	!~~~> facsafe: safety factor in the computation of new step size
	1399	IF (rcntrl(7) == zero)THEN
	1400	facsafe = 0.9_dp
	1401	ELSEIF (rcntrl(7)> zero)THEN
	1402	facsafe = rcntrl(7)
	1403	ELSE
	1404	PRINT *,'User-selected FacSafe: RCNTRL(7) =',RCNTRL(7)
	1405	CALL ros_errormsg(- 4,tstart,zero,ierr)
	1406	RETURN
	1407	ENDIF
	1408	!~~~> check IF tolerances are reasonable
	1409	DO i=1,uplimtol
	1410	IF((abstol(i)<= zero).or. (reltol(i)<= 10.0_dp*roundoff)&
	1411	.or. (reltol(i)>= 1.0_dp))THEN
	1412	PRINT *,' AbsTol(',i,') = ',AbsTol(i)
	1413	PRINT *,' RelTol(',i,') = ',RelTol(i)
	1414	CALL ros_errormsg(- 5,tstart,zero,ierr)
	1415	RETURN
	1416	ENDIF
	1417	ENDDO
	1418
	1419
	1420	!~~~> CALL rosenbrock method
	1421	CALL ros_integrator(y,tstart,tend,texit, &
	1422	abstol,reltol, &
	1423	! Integration parameters
	1424	autonomous,vectortol,max_no_steps, &
	1425	roundoff,hmin,hmax,hstart, &
	1426	facmin,facmax,facrej,facsafe, &
	1427	! Error indicator
	1428	ierr)
	1429
	1430	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1431	CONTAINS ! SUBROUTINEs internal to rosenbrock
	1432	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1433
	1434	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1435	SUBROUTINE ros_errormsg(code,t,h,ierr)
	1436	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1437	! Handles all error messages
	1438	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1439
	1440	REAL(kind=dp),INTENT(in):: t,h
	1441	INTEGER,INTENT(in) :: code
	1442	INTEGER,INTENT(out):: ierr
	1443
	1444	ierr = code
	1445	print *,&
	1446	'Forced exit from Rosenbrock due to the following error:'
	1447
	1448	select CASE (code)
	1449	CASE (- 1)
	1450	PRINT *,'--> Improper value for maximal no of steps'
	1451	CASE (- 2)
	1452	PRINT *,'--> Selected Rosenbrock method not implemented'
	1453	CASE (- 3)
	1454	PRINT *,'--> Hmin/Hmax/Hstart must be positive'
	1455	CASE (- 4)
	1456	PRINT *,'--> FacMin/FacMax/FacRej must be positive'
	1457	CASE (- 5)
	1458	PRINT *,'--> Improper tolerance values'
	1459	CASE (- 6)
	1460	PRINT *,'--> No of steps exceeds maximum bound'
	1461	CASE (- 7)
	1462	PRINT ,'--> Step size too small: T + 10H = T',&
	1463	' or H < Roundoff'
	1464	CASE (- 8)
	1465	PRINT *,'--> Matrix is repeatedly singular'
	1466	CASE default
	1467	PRINT *,'Unknown Error code: ',Code
	1468	END select
	1469
	1470	print *,"t=",t,"and h=",h
	1471
	1472	END SUBROUTINE ros_errormsg
	1473
	1474	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1475	SUBROUTINE ros_integrator (y,tstart,tend,t, &
	1476	abstol,reltol, &
	1477	!~~~> integration PARAMETERs
	1478	autonomous,vectortol,max_no_steps, &
	1479	roundoff,hmin,hmax,hstart, &
	1480	facmin,facmax,facrej,facsafe, &
	1481	!~~~> error indicator
	1482	ierr)
	1483	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1484	! Template for the implementation of a generic Rosenbrock method
	1485	! defined by ros_S (no of stages)
	1486	! and its coefficients ros_{A,C,M,E,Alpha,Gamma}
	1487	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1488
	1489
	1490	!~~~> input: the initial condition at tstart; output: the solution at t
	1491	REAL(kind=dp),INTENT(inout):: y(n)
	1492	!~~~> input: integration interval
	1493	REAL(kind=dp),INTENT(in):: tstart,tend
	1494	!~~~> output: time at which the solution is RETURNed (t=tENDIF success)
	1495	REAL(kind=dp),INTENT(out):: t
	1496	!~~~> input: tolerances
	1497	REAL(kind=dp),INTENT(in):: abstol(n),reltol(n)
	1498	!~~~> input: integration PARAMETERs
	1499	LOGICAL,INTENT(in):: autonomous,vectortol
	1500	REAL(kind=dp),INTENT(in):: hstart,hmin,hmax
	1501	INTEGER,INTENT(in):: max_no_steps
	1502	REAL(kind=dp),INTENT(in):: roundoff,facmin,facmax,facrej,facsafe
	1503	!~~~> output: error indicator
	1504	INTEGER,INTENT(out):: ierr
	1505	! ~~~~ Local variables
	1506	REAL(kind=dp):: ynew(n),fcn0(n),fcn(n)
	1507	REAL(kind=dp):: k(n*ros_s),dfdt(n)
	1508	#ifdef full_algebra
	1509	REAL(kind=dp):: jac0(n,n),ghimj(n,n)
	1510	#else
	1511	REAL(kind=dp):: jac0(lu_nonzero),ghimj(lu_nonzero)
	1512	#endif
	1513	REAL(kind=dp):: h,hnew,hc,hg,fac,tau
	1514	REAL(kind=dp):: err,yerr(n)
	1515	INTEGER :: pivot(n),direction,ioffset,j,istage
	1516	LOGICAL :: rejectlasth,rejectmoreh,singular
	1517	!~~~> local PARAMETERs
	1518	REAL(kind=dp),PARAMETER :: zero = 0.0_dp, one = 1.0_dp
	1519	REAL(kind=dp),PARAMETER :: deltamin = 1.0e-5_dp
	1520	!~~~> locally called FUNCTIONs
	1521	! REAL(kind=dp) WLAMCH
	1522	! EXTERNAL WLAMCH
	1523	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1524
	1525
	1526	!~~~> initial preparations
	1527	t = tstart
	1528	rstatus(nhexit) = zero
	1529	h = min( max(abs(hmin),abs(hstart)),abs(hmax))
	1530	IF (abs(h)<= 10.0_dp*roundoff)h = deltamin
	1531
	1532	IF (tEND >= tstart)THEN
	1533	direction = + 1
	1534	ELSE
	1535	direction = - 1
	1536	ENDIF
	1537	h = direction*h
	1538
	1539	rejectlasth=.false.
	1540	rejectmoreh=.false.
	1541
	1542	!~~~> time loop begins below
	1543
	1544	timeloop: DO WHILE((direction > 0).and.((t- tEND)+ roundoff <= zero)&
	1545	.or. (direction < 0).and.((tend-t)+ roundoff <= zero))
	1546
	1547	IF(istatus(nstp)> max_no_steps)THEN ! too many steps
	1548	CALL ros_errormsg(- 6,t,h,ierr)
	1549	RETURN
	1550	ENDIF
	1551	IF(((t+ 0.1_dp*h) == t).or.(h <= roundoff))THEN ! step size too small
	1552	CALL ros_errormsg(- 7,t,h,ierr)
	1553	RETURN
	1554	ENDIF
	1555
	1556	!~~~> limit h IF necessary to avoid going beyond tend
	1557	h = min(h,abs(tend-t))
	1558
	1559	!~~~> compute the FUNCTION at current time
	1560	CALL funtemplate(t,y,fcn0)
	1561	istatus(nfun) = istatus(nfun)+ 1
	1562
	1563	!~~~> compute the FUNCTION derivative with respect to t
	1564	IF (.not.autonomous)THEN
	1565	CALL ros_funtimederivative(t,roundoff,y,&
	1566	fcn0,dfdt)
	1567	ENDIF
	1568
	1569	!~~~> compute the jacobian at current time
	1570	CALL jactemplate(t,y,jac0)
	1571	istatus(njac) = istatus(njac)+ 1
	1572
	1573	!~~~> repeat step calculation until current step accepted
	1574	untilaccepted: do
	1575
	1576	CALL ros_preparematrix(h,direction,ros_gamma(1),&
	1577	jac0,ghimj,pivot,singular)
	1578	IF (singular)THEN ! more than 5 consecutive failed decompositions
	1579	CALL ros_errormsg(- 8,t,h,ierr)
	1580	RETURN
	1581	ENDIF
	1582
	1583	!~~~> compute the stages
	1584	stage: DO istage = 1,ros_s
	1585
	1586	! current istage offset. current istage vector is k(ioffset+ 1:ioffset+ n)
	1587	ioffset = n*(istage-1)
	1588
	1589	! for the 1st istage the FUNCTION has been computed previously
	1590	IF(istage == 1)THEN
	1591	!slim: CALL wcopy(n,fcn0,1,fcn,1)
	1592	fcn(1:n) = fcn0(1:n)
	1593	! istage>1 and a new FUNCTION evaluation is needed at the current istage
	1594	ELSEIF(ros_newf(istage))THEN
	1595	!slim: CALL wcopy(n,y,1,ynew,1)
	1596	ynew(1:n) = y(1:n)
	1597	DO j = 1,istage-1
	1598	CALL waxpy(n,ros_a((istage-1)*(istage-2)/2+ j),&
	1599	k(n*(j- 1)+ 1),1,ynew,1)
	1600	ENDDO
	1601	tau = t + ros_alpha(istage)directionh
	1602	CALL funtemplate(tau,ynew,fcn)
	1603	istatus(nfun) = istatus(nfun)+ 1
	1604	ENDIF ! IF istage == 1 ELSEIF ros_newf(istage)
	1605	!slim: CALL wcopy(n,fcn,1,k(ioffset+ 1),1)
	1606	k(ioffset+ 1:ioffset+ n) = fcn(1:n)
	1607	DO j = 1,istage-1
	1608	hc = ros_c((istage-1)(istage-2)/2+ j)/(directionh)
	1609	CALL waxpy(n,hc,k(n*(j- 1)+ 1),1,k(ioffset+ 1),1)
	1610	ENDDO
	1611	IF ((.not. autonomous).and.(ros_gamma(istage).ne.zero))THEN
	1612	hg = directionhros_gamma(istage)
	1613	CALL waxpy(n,hg,dfdt,1,k(ioffset+ 1),1)
	1614	ENDIF
	1615	CALL ros_solve(ghimj,pivot,k(ioffset+ 1))
	1616
	1617	END DO stage
	1618
	1619
	1620	!~~~> compute the new solution
	1621	!slim: CALL wcopy(n,y,1,ynew,1)
	1622	ynew(1:n) = y(1:n)
	1623	DO j=1,ros_s
	1624	CALL waxpy(n,ros_m(j),k(n*(j- 1)+ 1),1,ynew,1)
	1625	ENDDO
	1626
	1627	!~~~> compute the error estimation
	1628	!slim: CALL wscal(n,zero,yerr,1)
	1629	yerr(1:n) = zero
	1630	DO j=1,ros_s
	1631	CALL waxpy(n,ros_e(j),k(n*(j- 1)+ 1),1,yerr,1)
	1632	ENDDO
	1633	err = ros_errornorm(y,ynew,yerr,abstol,reltol,vectortol)
	1634
	1635	!~~~> new step size is bounded by facmin <= hnew/h <= facmax
	1636	fac = min(facmax,max(facmin,facsafe/err**(one/ros_elo)))
	1637	hnew = h*fac
	1638
	1639	!~~~> check the error magnitude and adjust step size
	1640	istatus(nstp) = istatus(nstp)+ 1
	1641	IF((err <= one).or.(h <= hmin))THEN !~~~> accept step
	1642	istatus(nacc) = istatus(nacc)+ 1
	1643	!slim: CALL wcopy(n,ynew,1,y,1)
	1644	y(1:n) = ynew(1:n)
	1645	t = t + direction*h
	1646	hnew = max(hmin,min(hnew,hmax))
	1647	IF (rejectlasth)THEN ! no step size increase after a rejected step
	1648	hnew = min(hnew,h)
	1649	ENDIF
	1650	rstatus(nhexit) = h
	1651	rstatus(nhnew) = hnew
	1652	rstatus(ntexit) = t
	1653	rejectlasth = .false.
	1654	rejectmoreh = .false.
	1655	h = hnew
	1656	exit untilaccepted ! exit the loop: WHILE step not accepted
	1657	ELSE !~~~> reject step
	1658	IF (rejectmoreh)THEN
	1659	hnew = h*facrej
	1660	ENDIF
	1661	rejectmoreh = rejectlasth
	1662	rejectlasth = .true.
	1663	h = hnew
	1664	IF (istatus(nacc)>= 1) istatus(nrej) = istatus(nrej)+ 1
	1665	ENDIF ! err <= 1
	1666
	1667	END DO untilaccepted
	1668
	1669	END DO timeloop
	1670
	1671	!~~~> succesful exit
	1672	ierr = 1 !~~~> the integration was successful
	1673
	1674	END SUBROUTINE ros_integrator
	1675
	1676
	1677	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1678	REAL(kind=dp)FUNCTION ros_errornorm(y,ynew,yerr,&
	1679	abstol,reltol,vectortol)
	1680	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1681	!~~~> computes the "scaled norm" of the error vector yerr
	1682	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1683
	1684	! Input arguments
	1685	REAL(kind=dp),INTENT(in):: y(n),ynew(n),&
	1686	yerr(n),abstol(n),reltol(n)
	1687	LOGICAL,INTENT(in):: vectortol
	1688	! Local variables
	1689	REAL(kind=dp):: err,scale,ymax
	1690	INTEGER :: i
	1691	REAL(kind=dp),PARAMETER :: zero = 0.0_dp
	1692
	1693	err = zero
	1694	DO i=1,n
	1695	ymax = max(abs(y(i)),abs(ynew(i)))
	1696	IF (vectortol)THEN
	1697	scale = abstol(i)+ reltol(i)*ymax
	1698	ELSE
	1699	scale = abstol(1)+ reltol(1)*ymax
	1700	ENDIF
	1701	err = err+(yerr(i)/scale)**2
	1702	ENDDO
	1703	err = sqrt(err/n)
	1704
	1705	ros_errornorm = max(err,1.0d-10)
	1706
	1707	END FUNCTION ros_errornorm
	1708
	1709
	1710	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1711	SUBROUTINE ros_funtimederivative(t,roundoff,y,&
	1712	fcn0,dfdt)
	1713	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1714	!~~~> the time partial derivative of the FUNCTION by finite differences
	1715	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1716
	1717	!~~~> input arguments
	1718	REAL(kind=dp),INTENT(in):: t,roundoff,y(n),fcn0(n)
	1719	!~~~> output arguments
	1720	REAL(kind=dp),INTENT(out):: dfdt(n)
	1721	!~~~> local variables
	1722	REAL(kind=dp):: delta
	1723	REAL(kind=dp),PARAMETER :: one = 1.0_dp, deltamin = 1.0e-6_dp
	1724
	1725	delta = sqrt(roundoff)*max(deltamin,abs(t))
	1726	CALL funtemplate(t+ delta,y,dfdt)
	1727	istatus(nfun) = istatus(nfun)+ 1
	1728	CALL waxpy(n,(- one),fcn0,1,dfdt,1)
	1729	CALL wscal(n,(one/delta),dfdt,1)
	1730
	1731	END SUBROUTINE ros_funtimederivative
	1732
	1733
	1734	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1735	SUBROUTINE ros_preparematrix(h,direction,gam,&
	1736	jac0,ghimj,pivot,singular)
	1737	! --- --- --- --- --- --- --- --- --- --- --- --- ---
	1738	! Prepares the LHS matrix for stage calculations
	1739	! 1. Construct Ghimj = 1/(H*ham) - Jac0
	1740	! "(Gamma H) Inverse Minus Jacobian"
	1741	! 2. Repeat LU decomposition of Ghimj until successful.
	1742	! -half the step size if LU decomposition fails and retry
	1743	! -exit after 5 consecutive fails
	1744	! --- --- --- --- --- --- --- --- --- --- --- --- ---
	1745
	1746	!~~~> input arguments
	1747	#ifdef full_algebra
	1748	REAL(kind=dp),INTENT(in):: jac0(n,n)
	1749	#else
	1750	REAL(kind=dp),INTENT(in):: jac0(lu_nonzero)
	1751	#endif
	1752	REAL(kind=dp),INTENT(in):: gam
	1753	INTEGER,INTENT(in):: direction
	1754	!~~~> output arguments
	1755	#ifdef full_algebra
	1756	REAL(kind=dp),INTENT(out):: ghimj(n,n)
	1757	#else
	1758	REAL(kind=dp),INTENT(out):: ghimj(lu_nonzero)
	1759	#endif
	1760	LOGICAL,INTENT(out):: singular
	1761	INTEGER,INTENT(out):: pivot(n)
	1762	!~~~> inout arguments
	1763	REAL(kind=dp),INTENT(inout):: h ! step size is decreased when lu fails
	1764	!~~~> local variables
	1765	INTEGER :: i,ising,nconsecutive
	1766	REAL(kind=dp):: ghinv
	1767	REAL(kind=dp),PARAMETER :: one = 1.0_dp, half = 0.5_dp
	1768
	1769	nconsecutive = 0
	1770	singular = .true.
	1771
	1772	DO WHILE (singular)
	1773
	1774	!~~~> construct ghimj = 1/(h*gam)- jac0
	1775	#ifdef full_algebra
	1776	!slim: CALL wcopy(n*n,jac0,1,ghimj,1)
	1777	!slim: CALL wscal(n*n,(- one),ghimj,1)
	1778	ghimj = - jac0
	1779	ghinv = one/(directionhgam)
	1780	DO i=1,n
	1781	ghimj(i,i) = ghimj(i,i)+ ghinv
	1782	ENDDO
	1783	#else
	1784	!slim: CALL wcopy(lu_nonzero,jac0,1,ghimj,1)
	1785	!slim: CALL wscal(lu_nonzero,(- one),ghimj,1)
	1786	ghimj(1:lu_nonzero) = - jac0(1:lu_nonzero)
	1787	ghinv = one/(directionhgam)
	1788	DO i=1,n
	1789	ghimj(lu_diag(i)) = ghimj(lu_diag(i))+ ghinv
	1790	ENDDO
	1791	#endif
	1792	!~~~> compute lu decomposition
	1793	CALL ros_decomp( ghimj,pivot,ising)
	1794	IF (ising == 0)THEN
	1795	!~~~> IF successful done
	1796	singular = .false.
	1797	ELSE ! ising .ne. 0
	1798	!~~~> IF unsuccessful half the step size; IF 5 consecutive fails THEN RETURN
	1799	istatus(nsng) = istatus(nsng)+ 1
	1800	nconsecutive = nconsecutive+1
	1801	singular = .true.
	1802	PRINT*,'Warning: LU Decomposition returned ISING = ',ISING
	1803	IF (nconsecutive <= 5)THEN ! less than 5 consecutive failed decompositions
	1804	h = h*half
	1805	ELSE ! more than 5 consecutive failed decompositions
	1806	RETURN
	1807	ENDIF ! nconsecutive
	1808	ENDIF ! ising
	1809
	1810	END DO ! WHILE singular
	1811
	1812	END SUBROUTINE ros_preparematrix
	1813
	1814
	1815	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1816	SUBROUTINE ros_decomp( a,pivot,ising)
	1817	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1818	! Template for the LU decomposition
	1819	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1820	!~~~> inout variables
	1821	#ifdef full_algebra
	1822	REAL(kind=dp),INTENT(inout):: a(n,n)
	1823	#else
	1824	REAL(kind=dp),INTENT(inout):: a(lu_nonzero)
	1825	#endif
	1826	!~~~> output variables
	1827	INTEGER,INTENT(out):: pivot(n),ising
	1828
	1829	#ifdef full_algebra
	1830	CALL dgetrf( n,n,a,n,pivot,ising)
	1831	#else
	1832	CALL kppdecomp(a,ising)
	1833	pivot(1) = 1
	1834	#endif
	1835	istatus(ndec) = istatus(ndec)+ 1
	1836
	1837	END SUBROUTINE ros_decomp
	1838
	1839
	1840	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1841	SUBROUTINE ros_solve( a,pivot,b)
	1842	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1843	! Template for the forward/backward substitution (using pre-computed LU decomposition)
	1844	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1845	!~~~> input variables
	1846	#ifdef full_algebra
	1847	REAL(kind=dp),INTENT(in):: a(n,n)
	1848	INTEGER :: ising
	1849	#else
	1850	REAL(kind=dp),INTENT(in):: a(lu_nonzero)
	1851	#endif
	1852	INTEGER,INTENT(in):: pivot(n)
	1853	!~~~> inout variables
	1854	REAL(kind=dp),INTENT(inout):: b(n)
	1855
	1856	#ifdef full_algebra
	1857	CALL DGETRS( 'N',N ,1,A,N,Pivot,b,N,ISING)
	1858	IF(info < 0)THEN
	1859	print*,"error in dgetrs. ising=",ising
	1860	ENDIF
	1861	#else
	1862	CALL kppsolve( a,b)
	1863	#endif
	1864
	1865	istatus(nsol) = istatus(nsol)+ 1
	1866
	1867	END SUBROUTINE ros_solve
	1868
	1869
	1870
	1871	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1872	SUBROUTINE ros2
	1873	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1874	! --- AN L-STABLE METHOD,2 stages,order 2
	1875	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1876
	1877	double precision g
	1878
	1879	g = 1.0_dp + 1.0_dp/sqrt(2.0_dp)
	1880	rosmethod = rs2
	1881	!~~~> name of the method
	1882	ros_Name = 'ROS-2'
	1883	!~~~> number of stages
	1884	ros_s = 2
	1885
	1886	!~~~> the coefficient matrices a and c are strictly lower triangular.
	1887	! The lower triangular (subdiagonal) elements are stored in row-wise order:
	1888	! A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
	1889	! The general mapping formula is:
	1890	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
	1891	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
	1892
	1893	ros_a(1) = (1.0_dp)/g
	1894	ros_c(1) = (- 2.0_dp)/g
	1895	!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
	1896	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
	1897	ros_newf(1) = .true.
	1898	ros_newf(2) = .true.
	1899	!~~~> m_i = coefficients for new step solution
	1900	ros_m(1) = (3.0_dp)/(2.0_dp*g)
	1901	ros_m(2) = (1.0_dp)/(2.0_dp*g)
	1902	! E_i = Coefficients for error estimator
	1903	ros_e(1) = 1.0_dp/(2.0_dp*g)
	1904	ros_e(2) = 1.0_dp/(2.0_dp*g)
	1905	!~~~> ros_elo = estimator of local order - the minimum between the
	1906	! main and the embedded scheme orders plus one
	1907	ros_elo = 2.0_dp
	1908	!~~~> y_stage_i ~ y( t + h*alpha_i)
	1909	ros_alpha(1) = 0.0_dp
	1910	ros_alpha(2) = 1.0_dp
	1911	!~~~> gamma_i = \sum_j gamma_{i,j}
	1912	ros_gamma(1) = g
	1913	ros_gamma(2) =- g
	1914
	1915	END SUBROUTINE ros2
	1916
	1917
	1918	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1919	SUBROUTINE ros3
	1920	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1921	! --- AN L-STABLE METHOD,3 stages,order 3,2 function evaluations
	1922	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1923
	1924	rosmethod = rs3
	1925	!~~~> name of the method
	1926	ros_Name = 'ROS-3'
	1927	!~~~> number of stages
	1928	ros_s = 3
	1929
	1930	!~~~> the coefficient matrices a and c are strictly lower triangular.
	1931	! The lower triangular (subdiagonal) elements are stored in row-wise order:
	1932	! A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
	1933	! The general mapping formula is:
	1934	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
	1935	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
	1936
	1937	ros_a(1) = 1.0_dp
	1938	ros_a(2) = 1.0_dp
	1939	ros_a(3) = 0.0_dp
	1940
	1941	ros_c(1) = - 0.10156171083877702091975600115545e+01_dp
	1942	ros_c(2) = 0.40759956452537699824805835358067e+01_dp
	1943	ros_c(3) = 0.92076794298330791242156818474003e+01_dp
	1944	!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
	1945	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
	1946	ros_newf(1) = .true.
	1947	ros_newf(2) = .true.
	1948	ros_newf(3) = .false.
	1949	!~~~> m_i = coefficients for new step solution
	1950	ros_m(1) = 0.1e+01_dp
	1951	ros_m(2) = 0.61697947043828245592553615689730e+01_dp
	1952	ros_m(3) = - 0.42772256543218573326238373806514_dp
	1953	! E_i = Coefficients for error estimator
	1954	ros_e(1) = 0.5_dp
	1955	ros_e(2) = - 0.29079558716805469821718236208017e+01_dp
	1956	ros_e(3) = 0.22354069897811569627360909276199_dp
	1957	!~~~> ros_elo = estimator of local order - the minimum between the
	1958	! main and the embedded scheme orders plus 1
	1959	ros_elo = 3.0_dp
	1960	!~~~> y_stage_i ~ y( t + h*alpha_i)
	1961	ros_alpha(1) = 0.0_dp
	1962	ros_alpha(2) = 0.43586652150845899941601945119356_dp
	1963	ros_alpha(3) = 0.43586652150845899941601945119356_dp
	1964	!~~~> gamma_i = \sum_j gamma_{i,j}
	1965	ros_gamma(1) = 0.43586652150845899941601945119356_dp
	1966	ros_gamma(2) = 0.24291996454816804366592249683314_dp
	1967	ros_gamma(3) = 0.21851380027664058511513169485832e+01_dp
	1968
	1969	END SUBROUTINE ros3
	1970
	1971	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1972
	1973
	1974	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1975	SUBROUTINE ros4
	1976	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1977	! L-STABLE ROSENBROCK METHOD OF ORDER 4,WITH 4 STAGES
	1978	! L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
	1979	!
	1980	! E. HAIRER AND G. WANNER,SOLVING ORDINARY DIFFERENTIAL
	1981	! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
	1982	! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
	1983	! SPRINGER-VERLAG (1990)
	1984	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1985
	1986
	1987	rosmethod = rs4
	1988	!~~~> name of the method
	1989	ros_Name = 'ROS-4'
	1990	!~~~> number of stages
	1991	ros_s = 4
	1992
	1993	!~~~> the coefficient matrices a and c are strictly lower triangular.
	1994	! The lower triangular (subdiagonal) elements are stored in row-wise order:
	1995	! A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
	1996	! The general mapping formula is:
	1997	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
	1998	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
	1999
	2000	ros_a(1) = 0.2000000000000000e+01_dp
	2001	ros_a(2) = 0.1867943637803922e+01_dp
	2002	ros_a(3) = 0.2344449711399156_dp
	2003	ros_a(4) = ros_a(2)
	2004	ros_a(5) = ros_a(3)
	2005	ros_a(6) = 0.0_dp
	2006
	2007	ros_c(1) =- 0.7137615036412310e+01_dp
	2008	ros_c(2) = 0.2580708087951457e+01_dp
	2009	ros_c(3) = 0.6515950076447975_dp
	2010	ros_c(4) =- 0.2137148994382534e+01_dp
	2011	ros_c(5) =- 0.3214669691237626_dp
	2012	ros_c(6) =- 0.6949742501781779_dp
	2013	!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
	2014	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
	2015	ros_newf(1) = .true.
	2016	ros_newf(2) = .true.
	2017	ros_newf(3) = .true.
	2018	ros_newf(4) = .false.
	2019	!~~~> m_i = coefficients for new step solution
	2020	ros_m(1) = 0.2255570073418735e+01_dp
	2021	ros_m(2) = 0.2870493262186792_dp
	2022	ros_m(3) = 0.4353179431840180_dp
	2023	ros_m(4) = 0.1093502252409163e+01_dp
	2024	!~~~> e_i = coefficients for error estimator
	2025	ros_e(1) =- 0.2815431932141155_dp
	2026	ros_e(2) =- 0.7276199124938920e-01_dp
	2027	ros_e(3) =- 0.1082196201495311_dp
	2028	ros_e(4) =- 0.1093502252409163e+01_dp
	2029	!~~~> ros_elo = estimator of local order - the minimum between the
	2030	! main and the embedded scheme orders plus 1
	2031	ros_elo = 4.0_dp
	2032	!~~~> y_stage_i ~ y( t + h*alpha_i)
	2033	ros_alpha(1) = 0.0_dp
	2034	ros_alpha(2) = 0.1145640000000000e+01_dp
	2035	ros_alpha(3) = 0.6552168638155900_dp
	2036	ros_alpha(4) = ros_alpha(3)
	2037	!~~~> gamma_i = \sum_j gamma_{i,j}
	2038	ros_gamma(1) = 0.5728200000000000_dp
	2039	ros_gamma(2) =- 0.1769193891319233e+01_dp
	2040	ros_gamma(3) = 0.7592633437920482_dp
	2041	ros_gamma(4) =- 0.1049021087100450_dp
	2042
	2043	END SUBROUTINE ros4
	2044
	2045	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	2046	SUBROUTINE rodas3
	2047	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	2048	! --- A STIFFLY-STABLE METHOD,4 stages,order 3
	2049	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	2050
	2051
	2052	rosmethod = rd3
	2053	!~~~> name of the method
	2054	ros_Name = 'RODAS-3'
	2055	!~~~> number of stages
	2056	ros_s = 4
	2057
	2058	!~~~> the coefficient matrices a and c are strictly lower triangular.
	2059	! The lower triangular (subdiagonal) elements are stored in row-wise order:
	2060	! A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
	2061	! The general mapping formula is:
	2062	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
	2063	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
	2064
	2065	ros_a(1) = 0.0_dp
	2066	ros_a(2) = 2.0_dp
	2067	ros_a(3) = 0.0_dp
	2068	ros_a(4) = 2.0_dp
	2069	ros_a(5) = 0.0_dp
	2070	ros_a(6) = 1.0_dp
	2071
	2072	ros_c(1) = 4.0_dp
	2073	ros_c(2) = 1.0_dp
	2074	ros_c(3) =- 1.0_dp
	2075	ros_c(4) = 1.0_dp
	2076	ros_c(5) =- 1.0_dp
	2077	ros_c(6) =- (8.0_dp/3.0_dp)
	2078
	2079	!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
	2080	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
	2081	ros_newf(1) = .true.
	2082	ros_newf(2) = .false.
	2083	ros_newf(3) = .true.
	2084	ros_newf(4) = .true.
	2085	!~~~> m_i = coefficients for new step solution
	2086	ros_m(1) = 2.0_dp
	2087	ros_m(2) = 0.0_dp
	2088	ros_m(3) = 1.0_dp
	2089	ros_m(4) = 1.0_dp
	2090	!~~~> e_i = coefficients for error estimator
	2091	ros_e(1) = 0.0_dp
	2092	ros_e(2) = 0.0_dp
	2093	ros_e(3) = 0.0_dp
	2094	ros_e(4) = 1.0_dp
	2095	!~~~> ros_elo = estimator of local order - the minimum between the
	2096	! main and the embedded scheme orders plus 1
	2097	ros_elo = 3.0_dp
	2098	!~~~> y_stage_i ~ y( t + h*alpha_i)
	2099	ros_alpha(1) = 0.0_dp
	2100	ros_alpha(2) = 0.0_dp
	2101	ros_alpha(3) = 1.0_dp
	2102	ros_alpha(4) = 1.0_dp
	2103	!~~~> gamma_i = \sum_j gamma_{i,j}
	2104	ros_gamma(1) = 0.5_dp
	2105	ros_gamma(2) = 1.5_dp
	2106	ros_gamma(3) = 0.0_dp
	2107	ros_gamma(4) = 0.0_dp
	2108
	2109	END SUBROUTINE rodas3
	2110
	2111	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	2112	SUBROUTINE rodas4
	2113	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	2114	! STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4,WITH 6 STAGES
	2115	!
	2116	! E. HAIRER AND G. WANNER,SOLVING ORDINARY DIFFERENTIAL
	2117	! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
	2118	! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
	2119	! SPRINGER-VERLAG (1996)
	2120	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	2121
	2122
	2123	rosmethod = rd4
	2124	!~~~> name of the method
	2125	ros_Name = 'RODAS-4'
	2126	!~~~> number of stages
	2127	ros_s = 6
	2128
	2129	!~~~> y_stage_i ~ y( t + h*alpha_i)
	2130	ros_alpha(1) = 0.000_dp
	2131	ros_alpha(2) = 0.386_dp
	2132	ros_alpha(3) = 0.210_dp
	2133	ros_alpha(4) = 0.630_dp
	2134	ros_alpha(5) = 1.000_dp
	2135	ros_alpha(6) = 1.000_dp
	2136
	2137	!~~~> gamma_i = \sum_j gamma_{i,j}
	2138	ros_gamma(1) = 0.2500000000000000_dp
	2139	ros_gamma(2) =- 0.1043000000000000_dp
	2140	ros_gamma(3) = 0.1035000000000000_dp
	2141	ros_gamma(4) =- 0.3620000000000023e-01_dp
	2142	ros_gamma(5) = 0.0_dp
	2143	ros_gamma(6) = 0.0_dp
	2144
	2145	!~~~> the coefficient matrices a and c are strictly lower triangular.
	2146	! The lower triangular (subdiagonal) elements are stored in row-wise order:
	2147	! A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
	2148	! The general mapping formula is: A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
	2149	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
	2150
	2151	ros_a(1) = 0.1544000000000000e+01_dp
	2152	ros_a(2) = 0.9466785280815826_dp
	2153	ros_a(3) = 0.2557011698983284_dp
	2154	ros_a(4) = 0.3314825187068521e+01_dp
	2155	ros_a(5) = 0.2896124015972201e+01_dp
	2156	ros_a(6) = 0.9986419139977817_dp
	2157	ros_a(7) = 0.1221224509226641e+01_dp
	2158	ros_a(8) = 0.6019134481288629e+01_dp
	2159	ros_a(9) = 0.1253708332932087e+02_dp
	2160	ros_a(10) =- 0.6878860361058950_dp
	2161	ros_a(11) = ros_a(7)
	2162	ros_a(12) = ros_a(8)
	2163	ros_a(13) = ros_a(9)
	2164	ros_a(14) = ros_a(10)
	2165	ros_a(15) = 1.0_dp
	2166
	2167	ros_c(1) =- 0.5668800000000000e+01_dp
	2168	ros_c(2) =- 0.2430093356833875e+01_dp
	2169	ros_c(3) =- 0.2063599157091915_dp
	2170	ros_c(4) =- 0.1073529058151375_dp
	2171	ros_c(5) =- 0.9594562251023355e+01_dp
	2172	ros_c(6) =- 0.2047028614809616e+02_dp
	2173	ros_c(7) = 0.7496443313967647e+01_dp
	2174	ros_c(8) =- 0.1024680431464352e+02_dp
	2175	ros_c(9) =- 0.3399990352819905e+02_dp
	2176	ros_c(10) = 0.1170890893206160e+02_dp
	2177	ros_c(11) = 0.8083246795921522e+01_dp
	2178	ros_c(12) =- 0.7981132988064893e+01_dp
	2179	ros_c(13) =- 0.3152159432874371e+02_dp
	2180	ros_c(14) = 0.1631930543123136e+02_dp
	2181	ros_c(15) =- 0.6058818238834054e+01_dp
	2182
	2183	!~~~> m_i = coefficients for new step solution
	2184	ros_m(1) = ros_a(7)
	2185	ros_m(2) = ros_a(8)
	2186	ros_m(3) = ros_a(9)
	2187	ros_m(4) = ros_a(10)
	2188	ros_m(5) = 1.0_dp
	2189	ros_m(6) = 1.0_dp
	2190
	2191	!~~~> e_i = coefficients for error estimator
	2192	ros_e(1) = 0.0_dp
	2193	ros_e(2) = 0.0_dp
	2194	ros_e(3) = 0.0_dp
	2195	ros_e(4) = 0.0_dp
	2196	ros_e(5) = 0.0_dp
	2197	ros_e(6) = 1.0_dp
	2198
	2199	!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
	2200	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
	2201	ros_newf(1) = .true.
	2202	ros_newf(2) = .true.
	2203	ros_newf(3) = .true.
	2204	ros_newf(4) = .true.
	2205	ros_newf(5) = .true.
	2206	ros_newf(6) = .true.
	2207
	2208	!~~~> ros_elo = estimator of local order - the minimum between the
	2209	! main and the embedded scheme orders plus 1
	2210	ros_elo = 4.0_dp
	2211
	2212	END SUBROUTINE rodas4
	2213	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	2214	SUBROUTINE rang3
	2215	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	2216	! STIFFLY-STABLE W METHOD OF ORDER 3,WITH 4 STAGES
	2217	!
	2218	! J. RANG and L. ANGERMANN
	2219	! NEW ROSENBROCK W-METHODS OF ORDER 3
	2220	! FOR PARTIAL DIFFERENTIAL ALGEBRAIC
	2221	! EQUATIONS OF INDEX 1
	2222	! BIT Numerical Mathematics (2005) 45: 761-787
	2223	! DOI: 10.1007/s10543-005-0035-y
	2224	! Table 4.1-4.2
	2225	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	2226
	2227
	2228	rosmethod = rg3
	2229	!~~~> name of the method
	2230	ros_Name = 'RANG-3'
	2231	!~~~> number of stages
	2232	ros_s = 4
	2233
	2234	ros_a(1) = 5.09052051067020d+00;
	2235	ros_a(2) = 5.09052051067020d+00;
	2236	ros_a(3) = 0.0d0;
	2237	ros_a(4) = 4.97628111010787d+00;
	2238	ros_a(5) = 2.77268164715849d-02;
	2239	ros_a(6) = 2.29428036027904d-01;
	2240
	2241	ros_c(1) = - 1.16790812312283d+01;
	2242	ros_c(2) = - 1.64057326467367d+01;
	2243	ros_c(3) = - 2.77268164715850d-01;
	2244	ros_c(4) = - 8.38103960500476d+00;
	2245	ros_c(5) = - 8.48328409199343d-01;
	2246	ros_c(6) = 2.87009860433106d-01;
	2247
	2248	ros_m(1) = 5.22582761233094d+00;
	2249	ros_m(2) = - 5.56971148154165d-01;
	2250	ros_m(3) = 3.57979469353645d-01;
	2251	ros_m(4) = 1.72337398521064d+00;
	2252
	2253	ros_e(1) = - 5.16845212784040d+00;
	2254	ros_e(2) = - 1.26351942603842d+00;
	2255	ros_e(3) = - 1.11022302462516d-16;
	2256	ros_e(4) = 2.22044604925031d-16;
	2257
	2258	ros_alpha(1) = 0.0d00;
	2259	ros_alpha(2) = 2.21878746765329d+00;
	2260	ros_alpha(3) = 2.21878746765329d+00;
	2261	ros_alpha(4) = 1.55392337535788d+00;
	2262
	2263	ros_gamma(1) = 4.35866521508459d-01;
	2264	ros_gamma(2) = - 1.78292094614483d+00;
	2265	ros_gamma(3) = - 2.46541900496934d+00;
	2266	ros_gamma(4) = - 8.05529997906370d-01;
	2267
	2268
	2269	!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
	2270	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
	2271	ros_newf(1) = .true.
	2272	ros_newf(2) = .true.
	2273	ros_newf(3) = .true.
	2274	ros_newf(4) = .true.
	2275
	2276	!~~~> ros_elo = estimator of local order - the minimum between the
	2277	! main and the embedded scheme orders plus 1
	2278	ros_elo = 3.0_dp
	2279
	2280	END SUBROUTINE rang3
	2281	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	2282
	2283	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	2284	! End of the set of internal Rosenbrock subroutines
	2285	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	2286	END SUBROUTINE rosenbrock
	2287
	2288	SUBROUTINE funtemplate( t,y,ydot)
	2289	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	2290	! Template for the ODE function call.
	2291	! Updates the rate coefficients (and possibly the fixed species) at each call
	2292	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	2293	!~~~> input variables
	2294	REAL(kind=dp):: t,y(nvar)
	2295	!~~~> output variables
	2296	REAL(kind=dp):: ydot(nvar)
	2297	!~~~> local variables
	2298	REAL(kind=dp):: told
	2299
	2300	told = time
	2301	time = t
	2302	CALL fun( y,fix,rconst,ydot)
	2303	time = told
	2304
	2305	END SUBROUTINE funtemplate
	2306
	2307	SUBROUTINE jactemplate( t,y,jcb)
	2308	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	2309	! Template for the ODE Jacobian call.
	2310	! Updates the rate coefficients (and possibly the fixed species) at each call
	2311	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	2312	!~~~> input variables
	2313	REAL(kind=dp):: t,y(nvar)
	2314	!~~~> output variables
	2315	#ifdef full_algebra
	2316	REAL(kind=dp):: jv(lu_nonzero),jcb(nvar,nvar)
	2317	#else
	2318	REAL(kind=dp):: jcb(lu_nonzero)
	2319	#endif
	2320	!~~~> local variables
	2321	REAL(kind=dp):: told
	2322	#ifdef full_algebra
	2323	INTEGER :: i,j
	2324	#endif
	2325
	2326	told = time
	2327	time = t
	2328	#ifdef full_algebra
	2329	CALL jac_sp(y,fix,rconst,jv)
	2330	DO j=1,nvar
	2331	DO i=1,nvar
	2332	jcb(i,j) = 0.0_dp
	2333	ENDDO
	2334	ENDDO
	2335	DO i=1,lu_nonzero
	2336	jcb(lu_irow(i),lu_icol(i)) = jv(i)
	2337	ENDDO
	2338	#else
	2339	CALL jac_sp( y,fix,rconst,jcb)
	2340	#endif
	2341	time = told
	2342
	2343	END SUBROUTINE jactemplate
	2344
	2345	SUBROUTINE chem_gasphase_integrate (time_step_len,conc,tempk,photo,ierrf,xnacc,xnrej,istatus,l_debug,pe)
	2346
	2347	IMPLICIT NONE
	2348
	2349	REAL(dp), INTENT(IN) :: time_step_len
	2350	REAL(dp), DIMENSION(:,:), INTENT(INOUT) :: conc
	2351	REAL(dp), DIMENSION(:,:), INTENT(INOUT) :: photo
	2352	REAL(dp), DIMENSION(:), INTENT(IN) :: tempk
	2353	INTEGER, INTENT(OUT), OPTIONAL :: ierrf(:)
	2354	INTEGER, INTENT(OUT), OPTIONAL :: xnacc(:)
	2355	INTEGER, INTENT(OUT), OPTIONAL :: xnrej(:)
	2356	INTEGER, INTENT(INOUT), OPTIONAL :: istatus(:)
	2357	INTEGER, INTENT(IN), OPTIONAL :: pe
	2358	LOGICAL, INTENT(IN), OPTIONAL :: l_debug
	2359
	2360	INTEGER :: k ! loop variable
	2361	REAL(dp) :: dt
	2362	INTEGER, DIMENSION(20) :: istatus_u
	2363	INTEGER :: ierr_u
	2364
	2365
	2366	if (present (istatus)) istatus = 0
	2367
	2368	DO k=1,vl_glo,vl_dim
	2369	is = k
	2370	ie = min(k+vl_dim-1,vl_glo)
	2371	vl = ie-is+1
	2372
	2373	c(:) = conc(is,:)
	2374
	2375	temp = tempk(is)
	2376
	2377	phot(:) = photo(is,:)
	2378
	2379	CALL update_rconst
	2380
	2381	dt = time_step_len
	2382
	2383	! integrate from t=0 to t=dt
	2384	CALL integrate(0._dp, dt,icntrl,rcntrl,istatus_u = istatus_u,ierr_u=ierr_u)
	2385
	2386
	2387	IF (PRESENT(l_debug) .AND. PRESENT(pe)) THEN
	2388	IF (l_debug) CALL error_output(conc(is,:),ierr_u,pe)
	2389	ENDIF
	2390	conc(is,:) = c(:)
	2391
	2392	! Return Diagnostic Information
	2393
	2394	if(PRESENT(ierrf)) ierrf(is) = ierr_u
	2395	if(PRESENT(xnacc)) xnacc(is) = istatus_u(4)
	2396	if(PRESENT(xnrej)) xnrej(is) = istatus_u(5)
	2397
	2398	if (PRESENT (istatus)) then
	2399	istatus(1:8) = istatus(1:8) + istatus_u(1:8)
	2400	ENDIF
	2401
	2402	END DO
	2403
	2404
	2405	! Deallocate input arrays
	2406
	2407	if (ALLOCATED(temp)) DEALLOCATE(temp)
	2408
	2409	data_loaded = .false.
	2410
	2411	RETURN
	2412	END SUBROUTINE chem_gasphase_integrate
	2413
	2414	SUBROUTINE fill_temp(status,array)
	2415
	2416	INTEGER, INTENT(OUT) :: status
	2417	REAL(dp), INTENT(IN), DIMENSION(:) :: array
	2418
	2419	status = 0
	2420	IF (.not. ALLOCATED(temp)) &
	2421	ALLOCATE(temp(size(array)))
	2422
	2423	IF (data_loaded .AND. (vl_glo /= size(array,1))) THEN
	2424	status = 1
	2425	RETURN
	2426	END IF
	2427
	2428	vl_glo = size(array,1)
	2429	temp = array
	2430	data_loaded = .TRUE.
	2431
	2432	RETURN
	2433
	2434	END SUBROUTINE fill_temp
[2668]	2435	#endif
[2657]	2436	END MODULE chem_gasphase_mod
	2437

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |