Home

Context Navigation

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

Last change on this file since 3780 was 3780, checked in by forkel, 5 years ago
removed read from unit 10 in chemistry_model_mod.f90, added get_mechanismname
File size: 85.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |