Home

Context Navigation

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

Last change on this file since 3298 was 3298, checked in by kanani, 6 years ago
Merge chemistry branch at r3297 to trunk
File size: 90.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |