Home

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_passive1/chem_gasphase_mod.f90 @ 3458

Last change on this file since 3458 was 3458, checked in by kanani, 5 years ago
Reintegrated fixes/changes from branch chemistry
File size: 78.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |