Home

Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

chem_gasphase_mod.f90

Last change on this file was 4841, checked in by forkel, 3 years ago
updated copyright statements for chem_gasphase_mod.f90. This must be done in UTIL/chemistry/gasphase_preproc/kpp4palm/templates/module_header and NOT in SOURCE/chem_gasphase_mod.f90.
File size: 87.9 KB

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

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

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_salsa+simple/chem_gasphase_mod.f90

Download in other formats: