Home

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_simple/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.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |