Home

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_smog/chem_gasphase_mod.f90 @ 3780

Last change on this file since 3780 was 3780, checked in by forkel, 5 years ago
removed read from unit 10 in chemistry_model_mod.f90, added get_mechanismname
File size: 90.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |