Home

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_salsa+phstat/chem_gasphase_mod.f90 @ 4370

Last change on this file since 4370 was 4370, checked in by raasch, 4 years ago
bugfixes for previous commit: unused variables removed, Temperton-fft usage on GPU, openacc porting of vector version of Obukhov length calculation, collective read switched off on NEC to avoid hanging; some vector directives added in prognostic equations to force vectorization on Intel19 compiler, configuration files for NEC Aurora added
File size: 82.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |