Home

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_simplep/chem_gasphase_mod.f90 @ 3833

Last change on this file since 3833 was 3833, checked in by forkel, 5 years ago
removed USE chem_gasphase_mod from chem_modules, apply USE chem_gasphase for nvar, nspec, cs_mech and spc_names instead
File size: 86.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |