Home

Context Navigation

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

Last change on this file since 2805 was 2766, checked in by kanani, 7 years ago
Removal of chem directive, plus minor changes
File size: 79.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |