Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |