Home

Context Navigation

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

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

Line
1	MODULE chem_gasphase_mod
2
3	! Mechanism: cbm4
4	!
5	!------------------------------------------------------------------------------!
6	!
7	! ****Module chem_gasphase_mod is automatically generated by kpp4palm ****
8	!
9	! *******Please do NOT change this Code,it will be ovewritten *******
10	!
11	!------------------------------------------------------------------------------!
12	! This file was created by KPP (http://people.cs.vt.edu/asandu/Software/Kpp/)
13	! and kpp4palm (created by Klaus Ketelsen). kpp4palm is an adapted version
14	! of KP4 (JÃ¶ckel,P.,Kerkweg,A.,Pozzer,A.,Sander,R.,Tost,H.,Riede,
15	! H.,Baumgaertner,A.,Gromov,S.,and Kern,B.,2010: Development cycle 2 of
16	! the Modular Earth Submodel System (MESSy2),Geosci. Model Dev.,3,717-752,
17	! https://doi.org/10.5194/gmd-3-717-2010). KP4 is part of the Modular Earth
18	! Submodel System (MESSy),which is is available under the GNU General Public
19	! License (GPL).
20	!
21	! KPP is free software; you can redistribute it and/or modify it under the terms
22	! of the General Public Licence as published by the Free Software Foundation;
23	! either version 2 of the License,or (at your option) any later version.
24	! KPP is distributed in the hope that it will be useful,but WITHOUT ANY WARRANTY;
25	! without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
26	! PURPOSE. See the GNU General Public Licence for more details.
27	!
28	!------------------------------------------------------------------------------!
29	! This file is part of the PALM model system.
30	!
31	! PALM is free software: you can redistribute it and/or modify it under the
32	! terms of the GNU General Public License as published by the Free Software
33	! Foundation,either version 3 of the License,or (at your option) any later
34	! version.
35	!
36	! PALM is distributed in the hope that it will be useful,but WITHOUT ANY
37	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
38	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
39	!
40	! You should have received a copy of the GNU General Public License along with
41	! PALM. If not,see <http://www.gnu.org/licenses/>.
42	!
43	! Copyright 1997-2019 Leibniz Universitaet Hannover
44	!--------------------------------------------------------------------------------!
45	!
46	!
47	! MODULE HEADER TEMPLATE
48	!
49	! Initial version (Nov. 2016,ketelsen),for later modifications of module_header
50	! see comments in kpp4palm/src/create_kpp_module.C
51
52	! Set kpp Double Precision to PALM Default Precision
53
54	USE kinds, ONLY: dp=>wp
55
56	USE pegrid, ONLY: myid, threads_per_task
57
58	IMPLICIT NONE
59	PRIVATE
60	!SAVE ! note: occurs again in automatically generated code ...
61
62	! PUBLIC :: IERR_NAMES
63
64	! PUBLIC :: SPC_NAMES,EQN_NAMES,EQN_TAGS,REQ_HET,REQ_AEROSOL,REQ_PHOTRAT &
65	! ,REQ_MCFCT,IP_MAX,jname
66
67	PUBLIC :: cs_mech
68	PUBLIC :: eqn_names, phot_names, spc_names
69	PUBLIC :: nmaxfixsteps
70	PUBLIC :: atol, rtol
71	PUBLIC :: nspec, nreact
72	PUBLIC :: temp
73	PUBLIC :: qvap
74	PUBLIC :: fakt
75	PUBLIC :: phot
76	PUBLIC :: rconst
77	PUBLIC :: nvar
78	PUBLIC :: nphot
79	PUBLIC :: vl_dim ! PUBLIC to ebable other MODULEs to distiguish between scalar and vec
80
81	PUBLIC :: initialize, integrate, update_rconst
82	PUBLIC :: chem_gasphase_integrate
83	PUBLIC :: initialize_kpp_ctrl
84	PUBLIC :: get_mechanismname
85
86	! END OF MODULE HEADER TEMPLATE
87
88	! Variables used for vector mode
89
90	LOGICAL, PARAMETER :: l_vector = .FALSE.
91	INTEGER, PARAMETER :: i_lu_di = 2
92	INTEGER, PARAMETER :: vl_dim = 1
93	INTEGER :: vl
94
95	INTEGER :: vl_glo
96	INTEGER :: is, ie
97
98
99	INTEGER, DIMENSION(vl_dim) :: kacc, krej
100	INTEGER, DIMENSION(vl_dim) :: ierrv
101	LOGICAL :: data_loaded = .FALSE.
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
117	! Time : Tue Mar 5 11:50:49 2019
118	! Working directory : /home/forkel-r/palmstuff/work/trunk20190305/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
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
130	INTEGER, PARAMETER :: nspec = 33
131	! NVAR - Number of Variable species
132	INTEGER, PARAMETER :: nvar = 32
133	! NVARACT - Number of Active species
134	INTEGER, PARAMETER :: nvaract = 32
135	! NFIX - Number of Fixed species
136	INTEGER, PARAMETER :: nfix = 1
137	! NREACT - Number of reactions
138	INTEGER, PARAMETER :: nreact = 81
139	! NVARST - Starting of variables in conc. vect.
140	INTEGER, PARAMETER :: nvarst = 1
141	! NFIXST - Starting of fixed in conc. vect.
142	INTEGER, PARAMETER :: nfixst = 33
143	! NONZERO - Number of nonzero entries in Jacobian
144	INTEGER, PARAMETER :: nonzero = 276
145	! LU_NONZERO - Number of nonzero entries in LU factoriz. of Jacobian
146	INTEGER, PARAMETER :: lu_nonzero = 300
147	! CNVAR - (NVAR+1) Number of elements in compressed row format
148	INTEGER, PARAMETER :: cnvar = 33
149	! CNEQN - (NREACT+1) Number stoicm elements in compressed col format
150	INTEGER, PARAMETER :: cneqn = 82
151	! NHESS - Length of Sparse Hessian
152	INTEGER, PARAMETER :: nhess = 258
153	! NMASS - Number of atoms to check mass balance
154	INTEGER, PARAMETER :: nmass = 1
155
156	! Index declaration for variable species in C and VAR
157	! VAR(ind_spc) = C(ind_spc)
158
159	INTEGER, PARAMETER, PUBLIC :: ind_o1d_cb4 = 1
160	INTEGER, PARAMETER, PUBLIC :: ind_h2o2 = 2
161	INTEGER, PARAMETER, PUBLIC :: ind_pan = 3
162	INTEGER, PARAMETER, PUBLIC :: ind_cro = 4
163	INTEGER, PARAMETER, PUBLIC :: ind_tol = 5
164	INTEGER, PARAMETER, PUBLIC :: ind_n2o5 = 6
165	INTEGER, PARAMETER, PUBLIC :: ind_xyl = 7
166	INTEGER, PARAMETER, PUBLIC :: ind_xo2n = 8
167	INTEGER, PARAMETER, PUBLIC :: ind_hono = 9
168	INTEGER, PARAMETER, PUBLIC :: ind_pna = 10
169	INTEGER, PARAMETER, PUBLIC :: ind_to2 = 11
170	INTEGER, PARAMETER, PUBLIC :: ind_hno3 = 12
171	INTEGER, PARAMETER, PUBLIC :: ind_ror = 13
172	INTEGER, PARAMETER, PUBLIC :: ind_cres = 14
173	INTEGER, PARAMETER, PUBLIC :: ind_mgly = 15
174	INTEGER, PARAMETER, PUBLIC :: ind_co = 16
175	INTEGER, PARAMETER, PUBLIC :: ind_eth = 17
176	INTEGER, PARAMETER, PUBLIC :: ind_xo2 = 18
177	INTEGER, PARAMETER, PUBLIC :: ind_open = 19
178	INTEGER, PARAMETER, PUBLIC :: ind_par = 20
179	INTEGER, PARAMETER, PUBLIC :: ind_hcho = 21
180	INTEGER, PARAMETER, PUBLIC :: ind_isop = 22
181	INTEGER, PARAMETER, PUBLIC :: ind_ole = 23
182	INTEGER, PARAMETER, PUBLIC :: ind_ald2 = 24
183	INTEGER, PARAMETER, PUBLIC :: ind_o3 = 25
184	INTEGER, PARAMETER, PUBLIC :: ind_no2 = 26
185	INTEGER, PARAMETER, PUBLIC :: ind_ho = 27
186	INTEGER, PARAMETER, PUBLIC :: ind_ho2 = 28
187	INTEGER, PARAMETER, PUBLIC :: ind_o = 29
188	INTEGER, PARAMETER, PUBLIC :: ind_no3 = 30
189	INTEGER, PARAMETER, PUBLIC :: ind_no = 31
190	INTEGER, PARAMETER, PUBLIC :: ind_c2o3 = 32
191
192	! Index declaration for fixed species in C
193	! C(ind_spc)
194
195	INTEGER, PARAMETER, PUBLIC :: ind_h2o = 33
196
197	! Index declaration for fixed species in FIX
198	! FIX(indf_spc) = C(ind_spc) = C(NVAR+indf_spc)
199
200	INTEGER, PARAMETER :: indf_h2o = 1
201
202	! NJVRP - Length of sparse Jacobian JVRP
203	INTEGER, PARAMETER :: njvrp = 134
204
205	! NSTOICM - Length of Sparse Stoichiometric Matrix
206	INTEGER, PARAMETER :: nstoicm = 329
207
208
209	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
210	!
211	! Global Data Module File
212	!
213	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
214	! (http://www.cs.vt.edu/~asandu/Software/KPP)
215	! KPP is distributed under GPL,the general public licence
216	! (http://www.gnu.org/copyleft/gpl.html)
217	! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
218	! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
219	! With important contributions from:
220	! M. Damian,Villanova University,USA
221	! R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
222	!
223	! File : chem_gasphase_mod_Global.f90
224	! Time : Tue Mar 5 11:50:49 2019
225	! Working directory : /home/forkel-r/palmstuff/work/trunk20190305/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
226	! Equation file : chem_gasphase_mod.kpp
227	! Output root filename : chem_gasphase_mod
228	!
229	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
230
231
232
233
234
235
236	! Declaration of global variables
237
238	! C - Concentration of all species
239	REAL(kind=dp):: c(nspec)
240	! VAR - Concentrations of variable species (global)
241	REAL(kind=dp):: var(nvar)
242	! FIX - Concentrations of fixed species (global)
243	REAL(kind=dp):: fix(nfix)
244	! VAR,FIX are chunks of array C
245	EQUIVALENCE( c(1), var(1))
246	EQUIVALENCE( c(33), fix(1))
247	! RCONST - Rate constants (global)
248	REAL(kind=dp):: rconst(nreact)
249	! TIME - Current integration time
250	REAL(kind=dp):: time
251	! TEMP - Temperature
252	REAL(kind=dp):: temp
253	! TSTART - Integration start time
254	REAL(kind=dp):: tstart
255	! ATOL - Absolute tolerance
256	REAL(kind=dp):: atol(nvar)
257	! RTOL - Relative tolerance
258	REAL(kind=dp):: rtol(nvar)
259	! STEPMIN - Lower bound for integration step
260	REAL(kind=dp):: stepmin
261	! CFACTOR - Conversion factor for concentration units
262	REAL(kind=dp):: cfactor
263
264	! INLINED global variable declarations
265
266	! QVAP - Water vapor
267	REAL(kind=dp):: qvap
268	! FAKT - Conversion factor
269	REAL(kind=dp):: fakt
270
271	! CS_MECH for check of mechanism name with namelist
272	CHARACTER(len=30):: cs_mech
273
274	! INLINED global variable declarations
275
276
277
278	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
279	!
280	! Sparse Jacobian Data Structures File
281	!
282	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
283	! (http://www.cs.vt.edu/~asandu/Software/KPP)
284	! KPP is distributed under GPL,the general public licence
285	! (http://www.gnu.org/copyleft/gpl.html)
286	! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
287	! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
288	! With important contributions from:
289	! M. Damian,Villanova University,USA
290	! R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
291	!
292	! File : chem_gasphase_mod_JacobianSP.f90
293	! Time : Tue Mar 5 11:50:49 2019
294	! Working directory : /home/forkel-r/palmstuff/work/trunk20190305/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
295	! Equation file : chem_gasphase_mod.kpp
296	! Output root filename : chem_gasphase_mod
297	!
298	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
299
300
301
302
303
304
305	! Sparse Jacobian Data
306
307
308	INTEGER, PARAMETER, DIMENSION(300):: lu_irow = (/ &
309	1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, &
310	4, 5, 5, 6, 6, 6, 7, 7, 8, 8, 8, 8, &
311	8, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 10, &
312	10, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, &
313	12, 12, 13, 13, 13, 13, 14, 14, 14, 14, 14, 14, &
314	14, 15, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, &
315	16, 16, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, &
316	18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, &
317	18, 18, 18, 18, 18, 18, 18, 18, 19, 19, 19, 19, &
318	19, 19, 19, 20, 20, 20, 20, 20, 20, 20, 20, 20, &
319	20, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, &
320	21, 21, 22, 22, 22, 22, 22, 23, 23, 23, 23, 23, &
321	23, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, &
322	24, 24, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, &
323	25, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, &
324	26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 26, 27, &
325	27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, &
326	27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, &
327	28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, &
328	28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, &
329	29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, 29, &
330	29, 29, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, &
331	30, 30, 30, 30, 30, 31, 31, 31, 31, 31, 31, 31, &
332	31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 32, &
333	32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32 /)
334
335	INTEGER, PARAMETER, DIMENSION(300):: lu_icol = (/ &
336	1, 25, 2, 27, 28, 3, 26, 32, 4, 14, 26, 27, &
337	30, 5, 27, 6, 26, 30, 7, 27, 8, 13, 20, 22, &
338	23, 27, 29, 30, 31, 9, 26, 27, 31, 10, 26, 27, &
339	28, 5, 7, 11, 27, 31, 6, 12, 14, 21, 24, 26, &
340	27, 30, 13, 20, 26, 27, 5, 7, 11, 14, 27, 30, &
341	31, 7, 15, 19, 22, 25, 27, 15, 16, 17, 19, 21, &
342	22, 23, 24, 25, 27, 29, 30, 17, 22, 25, 27, 29, &
343	5, 7, 13, 14, 15, 17, 18, 19, 20, 22, 23, 24, &
344	25, 26, 27, 28, 29, 30, 31, 32, 11, 14, 19, 25, &
345	27, 30, 31, 7, 13, 20, 22, 23, 25, 26, 27, 29, &
346	30, 17, 19, 21, 22, 23, 24, 25, 27, 28, 29, 30, &
347	31, 32, 22, 25, 27, 29, 30, 22, 23, 25, 27, 29, &
348	30, 13, 17, 19, 20, 22, 23, 24, 25, 26, 27, 29, &
349	30, 31, 17, 19, 22, 23, 25, 26, 27, 28, 29, 30, &
350	31, 3, 4, 6, 9, 10, 11, 13, 14, 18, 19, 20, &
351	22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 1, &
352	2, 5, 7, 9, 10, 12, 14, 15, 16, 17, 19, 20, &
353	21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, &
354	2, 5, 7, 10, 11, 13, 14, 15, 16, 17, 19, 20, &
355	21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, &
356	1, 17, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, &
357	31, 32, 6, 12, 14, 21, 22, 23, 24, 25, 26, 27, &
358	28, 29, 30, 31, 32, 8, 9, 11, 13, 18, 19, 20, &
359	22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 3, &
360	15, 19, 22, 24, 25, 26, 27, 28, 29, 30, 31, 32 /)
361
362	INTEGER, PARAMETER, DIMENSION(33):: lu_crow = (/ &
363	1, 3, 6, 9, 14, 16, 19, 21, 30, 34, 38, 43, &
364	51, 55, 62, 68, 80, 85, 105, 112, 122, 135, 140, 146, &
365	159, 170, 192, 217, 241, 255, 270, 288, 301 /)
366
367	INTEGER, PARAMETER, DIMENSION(33):: lu_diag = (/ &
368	1, 3, 6, 9, 14, 16, 19, 21, 30, 34, 40, 44, &
369	51, 58, 63, 69, 80, 91, 107, 114, 124, 135, 141, 152, &
370	163, 185, 211, 236, 251, 267, 286, 300, 301 /)
371
372
373
374	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
375	!
376	! Utility Data Module File
377	!
378	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
379	! (http://www.cs.vt.edu/~asandu/Software/KPP)
380	! KPP is distributed under GPL,the general public licence
381	! (http://www.gnu.org/copyleft/gpl.html)
382	! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
383	! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
384	! With important contributions from:
385	! M. Damian,Villanova University,USA
386	! R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
387	!
388	! File : chem_gasphase_mod_Monitor.f90
389	! Time : Tue Mar 5 11:50:49 2019
390	! Working directory : /home/forkel-r/palmstuff/work/trunk20190305/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
391	! Equation file : chem_gasphase_mod.kpp
392	! Output root filename : chem_gasphase_mod
393	!
394	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
395
396
397
398
399
400	CHARACTER(len=15), PARAMETER, DIMENSION(33):: spc_names = (/ &
401	'O1D_CB4 ','H2O2 ','PAN ',&
402	'CRO ','TOL ','N2O5 ',&
403	'XYL ','XO2N ','HONO ',&
404	'PNA ','TO2 ','HNO3 ',&
405	'ROR ','CRES ','MGLY ',&
406	'CO ','ETH ','XO2 ',&
407	'OPEN ','PAR ','HCHO ',&
408	'ISOP ','OLE ','ALD2 ',&
409	'O3 ','NO2 ','HO ',&
410	'HO2 ','O ','NO3 ',&
411	'NO ','C2O3 ','H2O ' /)
412
413	CHARACTER(len=100), PARAMETER, DIMENSION(30):: eqn_names_0 = (/ &
414	' NO2 --> O + NO ',&
415	' O3 --> O ',&
416	' O3 --> O1D_CB4 ',&
417	' NO3 --> 0.89 NO2 + 0.89 O + 0.11 NO ',&
418	' HONO --> HO + NO ',&
419	' H2O2 --> 2 HO ',&
420	' HCHO --> CO + 2 HO2 ',&
421	' HCHO --> CO ',&
422	' ALD2 --> CO + XO2 + HCHO + 2 HO2 ',&
423	' OPEN --> CO + HO2 + C2O3 ',&
424	' MGLY --> CO + HO2 + C2O3 ',&
425	' O --> O3 ',&
426	' O3 + NO --> NO2 ',&
427	' NO2 + O --> NO ',&
428	' NO2 + O --> NO3 ',&
429	' O + NO --> NO2 ',&
430	' O3 + NO2 --> NO3 ',&
431	' O1D_CB4 --> O ',&
432	' O1D_CB4 + H2O --> 2 HO ',&
433	' O3 + HO --> HO2 ',&
434	' O3 + HO2 --> HO ',&
435	' NO3 + NO --> 2 NO2 ',&
436	' NO2 + NO3 --> NO2 + NO ',&
437	' NO2 + NO3 --> N2O5 ',&
438	' N2O5 + H2O --> 2 HNO3 ',&
439	' N2O5 --> NO2 + NO3 ',&
440	' 2 NO --> 2 NO2 ',&
441	'NO2 + NO + H2O --> 2 HONO ',&
442	' HO + NO --> HONO ',&
443	' HONO + HO --> NO2 ' /)
444	CHARACTER(len=100), PARAMETER, DIMENSION(30):: eqn_names_1 = (/ &
445	' 2 HONO --> NO2 + NO ',&
446	' NO2 + HO --> HNO3 ',&
447	' HNO3 + HO --> NO3 ',&
448	' HO2 + NO --> NO2 + HO ',&
449	' NO2 + HO2 --> PNA ',&
450	' PNA --> NO2 + HO2 ',&
451	' PNA + HO --> NO2 ',&
452	' 2 HO2 --> H2O2 ',&
453	' 2 HO2 + H2O --> H2O2 ',&
454	' H2O2 + HO --> HO2 ',&
455	' CO + HO --> HO2 ',&
456	' HCHO + HO --> CO + HO2 ',&
457	' HCHO + O --> CO + HO + HO2 ',&
458	' HCHO + NO3 --> HNO3 + CO + HO2 ',&
459	' ALD2 + O --> HO + C2O3 ',&
460	' ALD2 + HO --> C2O3 ',&
461	' ALD2 + NO3 --> HNO3 + C2O3 ',&
462	' NO + C2O3 --> XO2 + HCHO + NO2 + HO2 ',&
463	' NO2 + C2O3 --> PAN ',&
464	' PAN --> NO2 + C2O3 ',&
465	' 2 C2O3 --> 2 XO2 + 2 HCHO + 2 HO2 ',&
466	' HO2 + C2O3 --> 0.79 XO2 + 0.79 HCHO + 0.79 HO + 0.79 HO2 ',&
467	' HO --> XO2 + HCHO + HO2 ',&
468	' PAR + HO --> 0.13 XO2N + 0.76 ROR + 0.87 XO2 - -0.11 PAR + 0.11 ALD2 ... etc. ',&
469	' ROR --> 0.04 XO2N + 0.02 ROR + 0.96 XO2 - -2.1 PAR + 1.1 ALD2 ... etc. ',&
470	' ROR --> HO2 ',&
471	' ROR + NO2 --> ',&
472	' OLE + O --> 0.02 XO2N + 0.3 CO + 0.28 XO2 + 0.22 PAR + 0.2 HCHO ... etc. ',&
473	' OLE + HO --> XO2 - PAR + HCHO + ALD2 + HO2 ',&
474	' OLE + O3 --> 0.33 CO + 0.22 XO2 - PAR + 0.74 HCHO + 0.5 ALD2 + 0.1 HO ... etc. ' /)
475	CHARACTER(len=100), PARAMETER, DIMENSION(21):: eqn_names_2 = (/ &
476	' OLE + NO3 --> 0.09 XO2N + 0.91 XO2 - PAR + HCHO + ALD2 + NO2 ',&
477	' ETH + O --> CO + 0.7 XO2 + HCHO + 0.3 HO + 1.7 HO2 ',&
478	' ETH + HO --> XO2 + 1.56 HCHO + 0.22 ALD2 + HO2 ',&
479	' ETH + O3 --> 0.42 CO + HCHO + 0.12 HO2 ',&
480	' TOL + HO --> 0.56 TO2 + 0.36 CRES + 0.08 XO2 + 0.44 HO2 ',&
481	' TO2 + NO --> 0.9 OPEN + 0.9 NO2 + 0.9 HO2 ',&
482	' TO2 --> CRES + HO2 ',&
483	' CRES + HO --> 0.4 CRO + 0.6 XO2 + 0.3 OPEN + 0.6 HO2 ',&
484	' CRES + NO3 --> CRO + HNO3 ',&
485	' CRO + NO2 --> ',&
486	' XYL + HO --> 0.3 TO2 + 0.2 CRES + 0.8 MGLY + 0.5 XO2 + 1.1 PAR + 0.7 HO2 ... etc. ',&
487	' OPEN + HO --> 2 CO + XO2 + HCHO + 2 HO2 + C2O3 ',&
488	' OPEN + O3 --> 0.2 MGLY + 0.69 CO + 0.03 XO2 + 0.7 HCHO + 0.03 ALD2 ... etc. ',&
489	' MGLY + HO --> XO2 + C2O3 ',&
490	' ISOP + O --> 0.5 CO + 0.45 ETH + 0.5 XO2 + 0.9 PAR + 0.55 OLE + 0.8 ALD2 ... etc. ',&
491	' ISOP + HO --> 0.13 XO2N + 0.4 MGLY + ETH + XO2 + HCHO + 0.2 ALD2 + 0.67 HO2 ... etc. ',&
492	' ISOP + O3 --> 0.2 MGLY + 0.06 CO + 0.55 ETH + 0.1 PAR + HCHO + 0.4 ALD2 ... etc. ',&
493	' ISOP + NO3 --> XO2N ',&
494	' XO2 + NO --> NO2 ',&
495	' 2 XO2 --> ',&
496	' XO2N + NO --> ' /)
497	CHARACTER(len=100), PARAMETER, DIMENSION(81):: eqn_names = (/&
498	eqn_names_0, eqn_names_1, eqn_names_2 /)
499
500	! INLINED global variables
501
502	! inline f90_data: declaration of global variables for photolysis
503	! REAL(kind=dp):: phot(nphot)must eventually be moved to global later for
504	INTEGER, PARAMETER :: nphot = 9
505	! phot photolysis frequencies
506	REAL(kind=dp):: phot(nphot)
507
508	INTEGER, PARAMETER, PUBLIC :: j_no2 = 1
509	INTEGER, PARAMETER, PUBLIC :: j_o33p = 2
510	INTEGER, PARAMETER, PUBLIC :: j_o31d = 3
511	INTEGER, PARAMETER, PUBLIC :: j_no3o = 4
512	INTEGER, PARAMETER, PUBLIC :: j_no3o2 = 5
513	INTEGER, PARAMETER, PUBLIC :: j_hono = 6
514	INTEGER, PARAMETER, PUBLIC :: j_h2o2 = 7
515	INTEGER, PARAMETER, PUBLIC :: j_ch2or = 8
516	INTEGER, PARAMETER, PUBLIC :: j_ch2om = 9
517
518	CHARACTER(len=15), PARAMETER, DIMENSION(nphot):: phot_names = (/ &
519	'J_NO2 ','J_O33P ','J_O31D ', &
520	'J_NO3O ','J_NO3O2 ','J_HONO ', &
521	'J_H2O2 ','J_HCHO_B ','J_HCHO_A '/)
522
523	! End INLINED global variables
524
525
526	! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
527
528	! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
529
530	! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
531
532	! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
533
534
535	! variable definations from individual module headers
536
537	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
538	!
539	! Initialization File
540	!
541	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
542	! (http://www.cs.vt.edu/~asandu/Software/KPP)
543	! KPP is distributed under GPL,the general public licence
544	! (http://www.gnu.org/copyleft/gpl.html)
545	! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
546	! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
547	! With important contributions from:
548	! M. Damian,Villanova University,USA
549	! R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
550	!
551	! File : chem_gasphase_mod_Initialize.f90
552	! Time : Tue Mar 5 11:50:49 2019
553	! Working directory : /home/forkel-r/palmstuff/work/trunk20190305/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
554	! Equation file : chem_gasphase_mod.kpp
555	! Output root filename : chem_gasphase_mod
556	!
557	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
558
559
560
561
562
563	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
564	!
565	! Numerical Integrator (Time-Stepping) File
566	!
567	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
568	! (http://www.cs.vt.edu/~asandu/Software/KPP)
569	! KPP is distributed under GPL,the general public licence
570	! (http://www.gnu.org/copyleft/gpl.html)
571	! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
572	! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
573	! With important contributions from:
574	! M. Damian,Villanova University,USA
575	! R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
576	!
577	! File : chem_gasphase_mod_Integrator.f90
578	! Time : Tue Mar 5 11:50:49 2019
579	! Working directory : /home/forkel-r/palmstuff/work/trunk20190305/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
580	! Equation file : chem_gasphase_mod.kpp
581	! Output root filename : chem_gasphase_mod
582	!
583	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
584
585
586
587	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
588	!
589	! INTEGRATE - Integrator routine
590	! Arguments :
591	! TIN - Start Time for Integration
592	! TOUT - End Time for Integration
593	!
594	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
595
596	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
597	! Rosenbrock - Implementation of several Rosenbrock methods: !
598	! *Ros2 !
599	! *Ros3 !
600	! *Ros4 !
601	! *Rodas3 !
602	! *Rodas4 !
603	! By default the code employs the KPP sparse linear algebra routines !
604	! Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK) !
605	! !
606	! (C) Adrian Sandu,August 2004 !
607	! Virginia Polytechnic Institute and State University !
608	! Contact: sandu@cs.vt.edu !
609	! Revised by Philipp Miehe and Adrian Sandu,May 2006 ! !
610	! This implementation is part of KPP - the Kinetic PreProcessor !
611	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
612
613
614	SAVE
615
616	!~~~> statistics on the work performed by the rosenbrock method
617	INTEGER, PARAMETER :: nfun=1, njac=2, nstp=3, nacc=4, &
618	nrej=5, ndec=6, nsol=7, nsng=8, &
619	ntexit=1, nhexit=2, nhnew = 3
620
621	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
622	!
623	! Linear Algebra Data and Routines File
624	!
625	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
626	! (http://www.cs.vt.edu/~asandu/Software/KPP)
627	! KPP is distributed under GPL,the general public licence
628	! (http://www.gnu.org/copyleft/gpl.html)
629	! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
630	! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
631	! With important contributions from:
632	! M. Damian,Villanova University,USA
633	! R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
634	!
635	! File : chem_gasphase_mod_LinearAlgebra.f90
636	! Time : Tue Mar 5 11:50:49 2019
637	! Working directory : /home/forkel-r/palmstuff/work/trunk20190305/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
638	! Equation file : chem_gasphase_mod.kpp
639	! Output root filename : chem_gasphase_mod
640	!
641	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
642
643
644
645
646
647
648	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
649	!
650	! The ODE Jacobian of Chemical Model File
651	!
652	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
653	! (http://www.cs.vt.edu/~asandu/Software/KPP)
654	! KPP is distributed under GPL,the general public licence
655	! (http://www.gnu.org/copyleft/gpl.html)
656	! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
657	! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
658	! With important contributions from:
659	! M. Damian,Villanova University,USA
660	! R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
661	!
662	! File : chem_gasphase_mod_Jacobian.f90
663	! Time : Tue Mar 5 11:50:49 2019
664	! Working directory : /home/forkel-r/palmstuff/work/trunk20190305/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
665	! Equation file : chem_gasphase_mod.kpp
666	! Output root filename : chem_gasphase_mod
667	!
668	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
669
670
671
672
673
674
675	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
676	!
677	! The ODE Function of Chemical Model File
678	!
679	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
680	! (http://www.cs.vt.edu/~asandu/Software/KPP)
681	! KPP is distributed under GPL,the general public licence
682	! (http://www.gnu.org/copyleft/gpl.html)
683	! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
684	! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
685	! With important contributions from:
686	! M. Damian,Villanova University,USA
687	! R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
688	!
689	! File : chem_gasphase_mod_Function.f90
690	! Time : Tue Mar 5 11:50:49 2019
691	! Working directory : /home/forkel-r/palmstuff/work/trunk20190305/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
692	! Equation file : chem_gasphase_mod.kpp
693	! Output root filename : chem_gasphase_mod
694	!
695	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
696
697
698
699
700
701	! A - Rate for each equation
702	REAL(kind=dp):: a(nreact)
703
704	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
705	!
706	! The Reaction Rates File
707	!
708	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
709	! (http://www.cs.vt.edu/~asandu/Software/KPP)
710	! KPP is distributed under GPL,the general public licence
711	! (http://www.gnu.org/copyleft/gpl.html)
712	! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
713	! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
714	! With important contributions from:
715	! M. Damian,Villanova University,USA
716	! R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
717	!
718	! File : chem_gasphase_mod_Rates.f90
719	! Time : Tue Mar 5 11:50:49 2019
720	! Working directory : /home/forkel-r/palmstuff/work/trunk20190305/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
721	! Equation file : chem_gasphase_mod.kpp
722	! Output root filename : chem_gasphase_mod
723	!
724	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
725
726
727
728
729
730	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
731	!
732	! Auxiliary Routines File
733	!
734	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
735	! (http://www.cs.vt.edu/~asandu/Software/KPP)
736	! KPP is distributed under GPL,the general public licence
737	! (http://www.gnu.org/copyleft/gpl.html)
738	! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
739	! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
740	! With important contributions from:
741	! M. Damian,Villanova University,USA
742	! R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
743	!
744	! File : chem_gasphase_mod_Util.f90
745	! Time : Tue Mar 5 11:50:49 2019
746	! Working directory : /home/forkel-r/palmstuff/work/trunk20190305/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
747	! Equation file : chem_gasphase_mod.kpp
748	! Output root filename : chem_gasphase_mod
749	!
750	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
751
752
753
754
755
756
757	! header MODULE initialize_kpp_ctrl_template
758
759	! notes:
760	! - l_vector is automatically defined by kp4
761	! - vl_dim is automatically defined by kp4
762	! - i_lu_di is automatically defined by kp4
763	! - wanted is automatically defined by xmecca
764	! - icntrl rcntrl are automatically defined by kpp
765	! - "USE messy_main_tools" is in MODULE_header of messy_mecca_kpp.f90
766	! - SAVE will be automatically added by kp4
767
768	!SAVE
769
770	! for fixed time step control
771	! ... max. number of fixed time steps (sum must be 1)
772	INTEGER, PARAMETER :: nmaxfixsteps = 50
773	! ... switch for fixed time stepping
774	LOGICAL, PUBLIC :: l_fixed_step = .FALSE.
775	INTEGER, PUBLIC :: nfsteps = 1
776	! ... number of kpp control PARAMETERs
777	INTEGER, PARAMETER, PUBLIC :: nkppctrl = 20
778	!
779	INTEGER, DIMENSION(nkppctrl), PUBLIC :: icntrl = 0
780	REAL(dp), DIMENSION(nkppctrl), PUBLIC :: rcntrl = 0.0_dp
781	REAL(dp), DIMENSION(nmaxfixsteps), PUBLIC :: t_steps = 0.0_dp
782
783	! END header MODULE initialize_kpp_ctrl_template
784
785
786	! Interface Block
787
788	INTERFACE initialize
789	MODULE PROCEDURE initialize
790	END INTERFACE initialize
791
792	INTERFACE integrate
793	MODULE PROCEDURE integrate
794	END INTERFACE integrate
795
796	INTERFACE fun
797	MODULE PROCEDURE fun
798	END INTERFACE fun
799
800	INTERFACE kppsolve
801	MODULE PROCEDURE kppsolve
802	END INTERFACE kppsolve
803
804	INTERFACE jac_sp
805	MODULE PROCEDURE jac_sp
806	END INTERFACE jac_sp
807
808	INTERFACE k_arr
809	MODULE PROCEDURE k_arr
810	END INTERFACE k_arr
811
812	INTERFACE update_rconst
813	MODULE PROCEDURE update_rconst
814	END INTERFACE update_rconst
815
816	INTERFACE arr2
817	MODULE PROCEDURE arr2
818	END INTERFACE arr2
819
820	INTERFACE initialize_kpp_ctrl
821	MODULE PROCEDURE initialize_kpp_ctrl
822	END INTERFACE initialize_kpp_ctrl
823
824	INTERFACE error_output
825	MODULE PROCEDURE error_output
826	END INTERFACE error_output
827
828	INTERFACE wscal
829	MODULE PROCEDURE wscal
830	END INTERFACE wscal
831
832	!INTERFACE not working INTERFACE waxpy
833	!INTERFACE not working MODULE PROCEDURE waxpy
834	!INTERFACE not working END INTERFACE waxpy
835
836	INTERFACE rosenbrock
837	MODULE PROCEDURE rosenbrock
838	END INTERFACE rosenbrock
839
840	INTERFACE funtemplate
841	MODULE PROCEDURE funtemplate
842	END INTERFACE funtemplate
843
844	INTERFACE jactemplate
845	MODULE PROCEDURE jactemplate
846	END INTERFACE jactemplate
847
848	INTERFACE kppdecomp
849	MODULE PROCEDURE kppdecomp
850	END INTERFACE kppdecomp
851
852	INTERFACE get_mechanismname
853	MODULE PROCEDURE get_mechanismname
854	END INTERFACE get_mechanismname
855
856	INTERFACE chem_gasphase_integrate
857	MODULE PROCEDURE chem_gasphase_integrate
858	END INTERFACE chem_gasphase_integrate
859
860
861	CONTAINS
862
863	SUBROUTINE initialize()
864
865
866	INTEGER :: j, k
867
868	INTEGER :: i
869	REAL(kind=dp):: x
870	k = is
871	cfactor = 1.000000e+00_dp
872
873	x = (0.) * cfactor
874	DO i = 1 , nvar
875	ENDDO
876
877	x = (0.) * cfactor
878	DO i = 1 , nfix
879	fix(i) = x
880	ENDDO
881
882	! constant rate coefficients
883	! END constant rate coefficients
884
885	! INLINED initializations
886
887	fix(indf_h2o) = qvap
888
889	! End INLINED initializations
890
891
892	END SUBROUTINE initialize
893
894	SUBROUTINE integrate( tin, tout, &
895	icntrl_u, rcntrl_u, istatus_u, rstatus_u, ierr_u)
896
897
898	REAL(kind=dp), INTENT(IN):: tin ! start time
899	REAL(kind=dp), INTENT(IN):: tout ! END time
900	! OPTIONAL input PARAMETERs and statistics
901	INTEGER, INTENT(IN), OPTIONAL :: icntrl_u(20)
902	REAL(kind=dp), INTENT(IN), OPTIONAL :: rcntrl_u(20)
903	INTEGER, INTENT(OUT), OPTIONAL :: istatus_u(20)
904	REAL(kind=dp), INTENT(OUT), OPTIONAL :: rstatus_u(20)
905	INTEGER, INTENT(OUT), OPTIONAL :: ierr_u
906
907	REAL(kind=dp):: rcntrl(20), rstatus(20)
908	INTEGER :: icntrl(20), istatus(20), ierr
909
910	INTEGER, SAVE :: ntotal = 0
911
912	icntrl(:) = 0
913	rcntrl(:) = 0.0_dp
914	istatus(:) = 0
915	rstatus(:) = 0.0_dp
916
917	!~~~> fine-tune the integrator:
918	icntrl(1) = 0 ! 0 - non- autonomous, 1 - autonomous
919	icntrl(2) = 0 ! 0 - vector tolerances, 1 - scalars
920
921	! IF OPTIONAL PARAMETERs are given, and IF they are >0,
922	! THEN they overwrite default settings.
923	IF (PRESENT(icntrl_u))THEN
924	WHERE(icntrl_u(:)> 0)icntrl(:) = icntrl_u(:)
925	ENDIF
926	IF (PRESENT(rcntrl_u))THEN
927	WHERE(rcntrl_u(:)> 0)rcntrl(:) = rcntrl_u(:)
928	ENDIF
929
930
931	CALL rosenbrock(nvar, var, tin, tout, &
932	atol, rtol, &
933	rcntrl, icntrl, rstatus, istatus, ierr)
934
935	!~~~> debug option: show no of steps
936	! ntotal = ntotal + istatus(nstp)
937	! PRINT*,'NSTEPS=',ISTATUS(Nstp),' (',Ntotal,')',' O3=',VAR(ind_O3)
938
939	stepmin = rstatus(nhexit)
940	! IF OPTIONAL PARAMETERs are given for output they
941	! are updated with the RETURN information
942	IF (PRESENT(istatus_u))istatus_u(:) = istatus(:)
943	IF (PRESENT(rstatus_u))rstatus_u(:) = rstatus(:)
944	IF (PRESENT(ierr_u)) ierr_u = ierr
945
946	END SUBROUTINE integrate
947
948	SUBROUTINE fun(v, f, rct, vdot)
949
950	! V - Concentrations of variable species (local)
951	REAL(kind=dp):: v(nvar)
952	! F - Concentrations of fixed species (local)
953	REAL(kind=dp):: f(nfix)
954	! RCT - Rate constants (local)
955	REAL(kind=dp):: rct(nreact)
956	! Vdot - Time derivative of variable species concentrations
957	REAL(kind=dp):: vdot(nvar)
958
959
960	! Computation of equation rates
961	a(1) = rct(1) * v(26)
962	a(2) = rct(2) * v(25)
963	a(3) = rct(3) * v(25)
964	a(4) = rct(4) * v(30)
965	a(5) = rct(5) * v(9)
966	a(6) = rct(6) * v(2)
967	a(7) = rct(7) * v(21)
968	a(8) = rct(8) * v(21)
969	a(9) = rct(9) * v(24)
970	a(10) = rct(10) * v(19)
971	a(11) = rct(11) * v(15)
972	a(12) = rct(12) * v(29)
973	a(13) = rct(13) * v(25) * v(31)
974	a(14) = rct(14) * v(26) * v(29)
975	a(15) = rct(15) * v(26) * v(29)
976	a(16) = rct(16) * v(29) * v(31)
977	a(17) = rct(17) * v(25) * v(26)
978	a(18) = rct(18) * v(1)
979	a(19) = rct(19) * v(1) * f(1)
980	a(20) = rct(20) * v(25) * v(27)
981	a(21) = rct(21) * v(25) * v(28)
982	a(22) = rct(22) * v(30) * v(31)
983	a(23) = rct(23) * v(26) * v(30)
984	a(24) = rct(24) * v(26) * v(30)
985	a(25) = rct(25) * v(6) * f(1)
986	a(26) = rct(26) * v(6)
987	a(27) = rct(27) * v(31) * v(31)
988	a(28) = rct(28) * v(26) * v(31) * f(1)
989	a(29) = rct(29) * v(27) * v(31)
990	a(30) = rct(30) * v(9) * v(27)
991	a(31) = rct(31) * v(9) * v(9)
992	a(32) = rct(32) * v(26) * v(27)
993	a(33) = rct(33) * v(12) * v(27)
994	a(34) = rct(34) * v(28) * v(31)
995	a(35) = rct(35) * v(26) * v(28)
996	a(36) = rct(36) * v(10)
997	a(37) = rct(37) * v(10) * v(27)
998	a(38) = rct(38) * v(28) * v(28)
999	a(39) = rct(39) * v(28) * v(28) * f(1)
1000	a(40) = rct(40) * v(2) * v(27)
1001	a(41) = rct(41) * v(16) * v(27)
1002	a(42) = rct(42) * v(21) * v(27)
1003	a(43) = rct(43) * v(21) * v(29)
1004	a(44) = rct(44) * v(21) * v(30)
1005	a(45) = rct(45) * v(24) * v(29)
1006	a(46) = rct(46) * v(24) * v(27)
1007	a(47) = rct(47) * v(24) * v(30)
1008	a(48) = rct(48) * v(31) * v(32)
1009	a(49) = rct(49) * v(26) * v(32)
1010	a(50) = rct(50) * v(3)
1011	a(51) = rct(51) * v(32) * v(32)
1012	a(52) = rct(52) * v(28) * v(32)
1013	a(53) = rct(53) * v(27)
1014	a(54) = rct(54) * v(20) * v(27)
1015	a(55) = rct(55) * v(13)
1016	a(56) = rct(56) * v(13)
1017	a(57) = rct(57) * v(13) * v(26)
1018	a(58) = rct(58) * v(23) * v(29)
1019	a(59) = rct(59) * v(23) * v(27)
1020	a(60) = rct(60) * v(23) * v(25)
1021	a(61) = rct(61) * v(23) * v(30)
1022	a(62) = rct(62) * v(17) * v(29)
1023	a(63) = rct(63) * v(17) * v(27)
1024	a(64) = rct(64) * v(17) * v(25)
1025	a(65) = rct(65) * v(5) * v(27)
1026	a(66) = rct(66) * v(11) * v(31)
1027	a(67) = rct(67) * v(11)
1028	a(68) = rct(68) * v(14) * v(27)
1029	a(69) = rct(69) * v(14) * v(30)
1030	a(70) = rct(70) * v(4) * v(26)
1031	a(71) = rct(71) * v(7) * v(27)
1032	a(72) = rct(72) * v(19) * v(27)
1033	a(73) = rct(73) * v(19) * v(25)
1034	a(74) = rct(74) * v(15) * v(27)
1035	a(75) = rct(75) * v(22) * v(29)
1036	a(76) = rct(76) * v(22) * v(27)
1037	a(77) = rct(77) * v(22) * v(25)
1038	a(78) = rct(78) * v(22) * v(30)
1039	a(79) = rct(79) * v(18) * v(31)
1040	a(80) = rct(80) * v(18) * v(18)
1041	a(81) = rct(81) * v(8) * v(31)
1042
1043	! Aggregate function
1044	vdot(1) = a(3) - a(18) - a(19)
1045	vdot(2) = - a(6) + a(38) + a(39) - a(40)
1046	vdot(3) = a(49) - a(50)
1047	vdot(4) = 0.4* a(68) + a(69) - a(70)
1048	vdot(5) = - a(65)
1049	vdot(6) = a(24) - a(25) - a(26)
1050	vdot(7) = - a(71)
1051	vdot(8) = 0.13* a(54) + 0.04* a(55) + 0.02* a(58) + 0.09* a(61) + 0.13* a(76) + a(78) - a(81)
1052	vdot(9) = - a(5) + 2* a(28) + a(29) - a(30) - 2* a(31)
1053	vdot(10) = a(35) - a(36) - a(37)
1054	vdot(11) = 0.56* a(65) - a(66) - a(67) + 0.3* a(71)
1055	vdot(12) = 2* a(25) + a(32) - a(33) + a(44) + a(47) + a(69)
1056	vdot(13) = 0.76* a(54) - 0.98* a(55) - a(56) - a(57)
1057	vdot(14) = 0.36* a(65) + a(67) - a(68) - a(69) + 0.2* a(71)
1058	vdot(15) = - a(11) + 0.8* a(71) + 0.2* a(73) - a(74) + 0.4* a(76) + 0.2* a(77)
1059	vdot(16) = a(7) + a(8) + a(9) + a(10) + a(11) - a(41) + a(42) + a(43) + a(44) + 0.3* a(58) + 0.33* a(60) + a(62) + 0.42* a(64) &
1060	+ 2* a(72) + 0.69* a(73)&
1061	&+ 0.5* a(75) + 0.06* a(77)
1062	vdot(17) = - a(62) - a(63) - a(64) + 0.45* a(75) + a(76) + 0.55* a(77)
1063	vdot(18) = a(9) + a(48) + 2* a(51) + 0.79* a(52) + a(53) + 0.87* a(54) + 0.96* a(55) + 0.28* a(58) + a(59) + 0.22* a(60) + &
1064	0.91* a(61) + 0.7* a(62)&
1065	&+ a(63) + 0.08* a(65) + 0.6* a(68) + 0.5* a(71) + a(72) + 0.03* a(73) + a(74) + 0.5* a(75) + a(76) - a(79) - 2* &
1066	a(80)
1067	vdot(19) = - a(10) + 0.9* a(66) + 0.3* a(68) - a(72) - a(73)
1068	vdot(20) = - 1.11* a(54) - 2.1* a(55) + 0.22* a(58) - a(59) - a(60) - a(61) + 1.1* a(71) + 0.9* a(75) + 0.1* a(77)
1069	vdot(21) = - a(7) - a(8) + a(9) - a(42) - a(43) - a(44) + a(48) + 2* a(51) + 0.79* a(52) + a(53) + 0.2* a(58) + a(59) + 0.74* &
1070	a(60) + a(61) + a(62)&
1071	&+ 1.56* a(63) + a(64) + a(72) + 0.7* a(73) + a(76) + a(77)
1072	vdot(22) = - a(75) - a(76) - a(77) - a(78)
1073	vdot(23) = - a(58) - a(59) - a(60) - a(61) + 0.55* a(75)
1074	vdot(24) = - a(9) - a(45) - a(46) - a(47) + 0.11* a(54) + 1.1* a(55) + 0.63* a(58) + a(59) + 0.5* a(60) + a(61) + 0.22* a(63) + &
1075	0.03* a(73) + 0.8&
1076	&* a(75) + 0.2* a(76) + 0.4* a(77)
1077	vdot(25) = - a(2) - a(3) + a(12) - a(13) - a(17) - a(20) - a(21) - a(60) - a(64) - a(73) - a(77)
1078	vdot(26) = - a(1) + 0.89* a(4) + a(13) - a(14) - a(15) + a(16) - a(17) + 2* a(22) - a(24) + a(26) + 2* a(27) - a(28) + a(30) + &
1079	a(31) - a(32) + a(34)&
1080	&- a(35) + a(36) + a(37) + a(48) - a(49) + a(50) - a(57) + a(61) + 0.9* a(66) - a(70) + a(79)
1081	vdot(27) = a(5) + 2* a(6) + 2* a(19) - a(20) + a(21) - a(29) - a(30) - a(32) - a(33) + a(34) - a(37) - a(40) - a(41) - a(42) + &
1082	a(43) + a(45) - a(46)&
1083	&+ 0.79* a(52) - a(53) - a(54) + 0.2* a(58) - a(59) + 0.1* a(60) + 0.3* a(62) - a(63) - a(65) - a(68) - a(71) - &
1084	a(72) + 0.08* a(73) - a(74)&
1085	&- a(76) + 0.1* a(77)
1086	vdot(28) = 2* a(7) + 2* a(9) + a(10) + a(11) + a(20) - a(21) - a(34) - a(35) + a(36) - 2* a(38) - 2* a(39) + a(40) + a(41) + &
1087	a(42) + a(43) + a(44) + a(48)&
1088	&+ 2* a(51) - 0.21* a(52) + a(53) + 0.11* a(54) + 0.94* a(55) + a(56) + 0.38* a(58) + a(59) + 0.44* a(60) + 1.7* &
1089	a(62) + a(63) + 0.12* a(64)&
1090	&+ 0.44* a(65) + 0.9* a(66) + a(67) + 0.6* a(68) + 0.7* a(71) + 2* a(72) + 0.76* a(73) + 0.6* a(75) + 0.67* a(76) &
1091	+ 0.44* a(77)
1092	vdot(29) = a(1) + a(2) + 0.89* a(4) - a(12) - a(14) - a(15) - a(16) + a(18) - a(43) - a(45) - a(58) - a(62) - a(75)
1093	vdot(30) = - a(4) + a(15) + a(17) - a(22) - a(23) - a(24) + a(26) + a(33) - a(44) - a(47) - a(61) - a(69) - a(78)
1094	vdot(31) = a(1) + 0.11* a(4) + a(5) - a(13) + a(14) - a(16) - a(22) + a(23) - 2* a(27) - a(28) - a(29) + a(31) - a(34) - a(48) &
1095	- a(66) - a(79) - a(81)
1096	vdot(32) = a(10) + a(11) + a(45) + a(46) + a(47) - a(48) - a(49) + a(50) - 2* a(51) - a(52) + a(72) + 0.62* a(73) + a(74) + &
1097	0.2* a(76)
1098
1099	END SUBROUTINE fun
1100
1101	SUBROUTINE kppsolve(jvs, x)
1102
1103	! JVS - sparse Jacobian of variables
1104	REAL(kind=dp):: jvs(lu_nonzero)
1105	! X - Vector for variables
1106	REAL(kind=dp):: x(nvar)
1107
1108	x(11) = x(11) - jvs(38) * x(5) - jvs(39) * x(7)
1109	x(12) = x(12) - jvs(43) * x(6)
1110	x(14) = x(14) - jvs(55) * x(5) - jvs(56) * x(7) - jvs(57) * x(11)
1111	x(15) = x(15) - jvs(62) * x(7)
1112	x(16) = x(16) - jvs(68) * x(15)
1113	x(18) = x(18) - jvs(85) * x(5) - jvs(86) * x(7) - jvs(87) * x(13) - jvs(88) * x(14) - jvs(89) * x(15) - jvs(90) * x(17)
1114	x(19) = x(19) - jvs(105) * x(11) - jvs(106) * x(14)
1115	x(20) = x(20) - jvs(112) * x(7) - jvs(113) * x(13)
1116	x(21) = x(21) - jvs(122) * x(17) - jvs(123) * x(19)
1117	x(23) = x(23) - jvs(140) * x(22)
1118	x(24) = x(24) - jvs(146) * x(13) - jvs(147) * x(17) - jvs(148) * x(19) - jvs(149) * x(20) - jvs(150) * x(22) - jvs(151) * x(23)
1119	x(25) = x(25) - jvs(159) * x(17) - jvs(160) * x(19) - jvs(161) * x(22) - jvs(162) * x(23)
1120	x(26) = x(26) - jvs(170) * x(3) - jvs(171) * x(4) - jvs(172) * x(6) - jvs(173) * x(9) - jvs(174) * x(10) - jvs(175) * x(11) - &
1121	jvs(176) * x(13)&
1122	&- jvs(177) * x(14) - jvs(178) * x(18) - jvs(179) * x(19) - jvs(180) * x(20) - jvs(181) * x(22) - jvs(182) * x(23) - &
1123	jvs(183) * x(24)&
1124	&- jvs(184) * x(25)
1125	x(27) = x(27) - jvs(192) * x(1) - jvs(193) * x(2) - jvs(194) * x(5) - jvs(195) * x(7) - jvs(196) * x(9) - jvs(197) * x(10) - &
1126	jvs(198) * x(12)&
1127	&- jvs(199) * x(14) - jvs(200) * x(15) - jvs(201) * x(16) - jvs(202) * x(17) - jvs(203) * x(19) - jvs(204) * x(20) - &
1128	jvs(205) * x(21)&
1129	&- jvs(206) * x(22) - jvs(207) * x(23) - jvs(208) * x(24) - jvs(209) * x(25) - jvs(210) * x(26)
1130	x(28) = x(28) - jvs(217) * x(2) - jvs(218) * x(5) - jvs(219) * x(7) - jvs(220) * x(10) - jvs(221) * x(11) - jvs(222) * x(13) - &
1131	jvs(223) * x(14)&
1132	&- jvs(224) * x(15) - jvs(225) * x(16) - jvs(226) * x(17) - jvs(227) * x(19) - jvs(228) * x(20) - jvs(229) * x(21) - &
1133	jvs(230) * x(22)&
1134	&- jvs(231) * x(23) - jvs(232) * x(24) - jvs(233) * x(25) - jvs(234) * x(26) - jvs(235) * x(27)
1135	x(29) = x(29) - jvs(241) * x(1) - jvs(242) * x(17) - jvs(243) * x(21) - jvs(244) * x(22) - jvs(245) * x(23) - jvs(246) * x(24) &
1136	- jvs(247) * x(25)&
1137	&- jvs(248) * x(26) - jvs(249) * x(27) - jvs(250) * x(28)
1138	x(30) = x(30) - jvs(255) * x(6) - jvs(256) * x(12) - jvs(257) * x(14) - jvs(258) * x(21) - jvs(259) * x(22) - jvs(260) * x(23) &
1139	- jvs(261) * x(24)&
1140	&- jvs(262) * x(25) - jvs(263) * x(26) - jvs(264) * x(27) - jvs(265) * x(28) - jvs(266) * x(29)
1141	x(31) = x(31) - jvs(270) * x(8) - jvs(271) * x(9) - jvs(272) * x(11) - jvs(273) * x(13) - jvs(274) * x(18) - jvs(275) * x(19) - &
1142	jvs(276) * x(20)&
1143	&- jvs(277) * x(22) - jvs(278) * x(23) - jvs(279) * x(24) - jvs(280) * x(25) - jvs(281) * x(26) - jvs(282) * x(27) - &
1144	jvs(283) * x(28)&
1145	&- jvs(284) * x(29) - jvs(285) * x(30)
1146	x(32) = x(32) - jvs(288) * x(3) - jvs(289) * x(15) - jvs(290) * x(19) - jvs(291) * x(22) - jvs(292) * x(24) - jvs(293) * x(25) &
1147	- jvs(294) * x(26)&
1148	&- jvs(295) * x(27) - jvs(296) * x(28) - jvs(297) * x(29) - jvs(298) * x(30) - jvs(299) * x(31)
1149	x(32) = x(32) / jvs(300)
1150	x(31) = (x(31) - jvs(287) * x(32)) /(jvs(286))
1151	x(30) = (x(30) - jvs(268) * x(31) - jvs(269) * x(32)) /(jvs(267))
1152	x(29) = (x(29) - jvs(252) * x(30) - jvs(253) * x(31) - jvs(254) * x(32)) /(jvs(251))
1153	x(28) = (x(28) - jvs(237) * x(29) - jvs(238) * x(30) - jvs(239) * x(31) - jvs(240) * x(32)) /(jvs(236))
1154	x(27) = (x(27) - jvs(212) * x(28) - jvs(213) * x(29) - jvs(214) * x(30) - jvs(215) * x(31) - jvs(216) * x(32)) /(jvs(211))
1155	x(26) = (x(26) - jvs(186) * x(27) - jvs(187) * x(28) - jvs(188) * x(29) - jvs(189) * x(30) - jvs(190) * x(31) - jvs(191) * &
1156	x(32)) /(jvs(185))
1157	x(25) = (x(25) - jvs(164) * x(26) - jvs(165) * x(27) - jvs(166) * x(28) - jvs(167) * x(29) - jvs(168) * x(30) - jvs(169) * &
1158	x(31)) /(jvs(163))
1159	x(24) = (x(24) - jvs(153) * x(25) - jvs(154) * x(26) - jvs(155) * x(27) - jvs(156) * x(29) - jvs(157) * x(30) - jvs(158) * &
1160	x(31)) /(jvs(152))
1161	x(23) = (x(23) - jvs(142) * x(25) - jvs(143) * x(27) - jvs(144) * x(29) - jvs(145) * x(30)) /(jvs(141))
1162	x(22) = (x(22) - jvs(136) * x(25) - jvs(137) * x(27) - jvs(138) * x(29) - jvs(139) * x(30)) /(jvs(135))
1163	x(21) = (x(21) - jvs(125) * x(22) - jvs(126) * x(23) - jvs(127) * x(24) - jvs(128) * x(25) - jvs(129) * x(27) - jvs(130) * &
1164	x(28) - jvs(131)&
1165	&* x(29) - jvs(132) * x(30) - jvs(133) * x(31) - jvs(134) * x(32)) /(jvs(124))
1166	x(20) = (x(20) - jvs(115) * x(22) - jvs(116) * x(23) - jvs(117) * x(25) - jvs(118) * x(26) - jvs(119) * x(27) - jvs(120) * &
1167	x(29) - jvs(121)&
1168	&* x(30)) /(jvs(114))
1169	x(19) = (x(19) - jvs(108) * x(25) - jvs(109) * x(27) - jvs(110) * x(30) - jvs(111) * x(31)) /(jvs(107))
1170	x(18) = (x(18) - jvs(92) * x(19) - jvs(93) * x(20) - jvs(94) * x(22) - jvs(95) * x(23) - jvs(96) * x(24) - jvs(97) * x(25) - &
1171	jvs(98) * x(26)&
1172	&- jvs(99) * x(27) - jvs(100) * x(28) - jvs(101) * x(29) - jvs(102) * x(30) - jvs(103) * x(31) - jvs(104) * x(32)) &
1173	/(jvs(91))
1174	x(17) = (x(17) - jvs(81) * x(22) - jvs(82) * x(25) - jvs(83) * x(27) - jvs(84) * x(29)) /(jvs(80))
1175	x(16) = (x(16) - jvs(70) * x(17) - jvs(71) * x(19) - jvs(72) * x(21) - jvs(73) * x(22) - jvs(74) * x(23) - jvs(75) * x(24) - &
1176	jvs(76) * x(25)&
1177	&- jvs(77) * x(27) - jvs(78) * x(29) - jvs(79) * x(30)) /(jvs(69))
1178	x(15) = (x(15) - jvs(64) * x(19) - jvs(65) * x(22) - jvs(66) * x(25) - jvs(67) * x(27)) /(jvs(63))
1179	x(14) = (x(14) - jvs(59) * x(27) - jvs(60) * x(30) - jvs(61) * x(31)) /(jvs(58))
1180	x(13) = (x(13) - jvs(52) * x(20) - jvs(53) * x(26) - jvs(54) * x(27)) /(jvs(51))
1181	x(12) = (x(12) - jvs(45) * x(14) - jvs(46) * x(21) - jvs(47) * x(24) - jvs(48) * x(26) - jvs(49) * x(27) - jvs(50) * x(30)) &
1182	/(jvs(44))
1183	x(11) = (x(11) - jvs(41) * x(27) - jvs(42) * x(31)) /(jvs(40))
1184	x(10) = (x(10) - jvs(35) * x(26) - jvs(36) * x(27) - jvs(37) * x(28)) /(jvs(34))
1185	x(9) = (x(9) - jvs(31) * x(26) - jvs(32) * x(27) - jvs(33) * x(31)) /(jvs(30))
1186	x(8) = (x(8) - jvs(22) * x(13) - jvs(23) * x(20) - jvs(24) * x(22) - jvs(25) * x(23) - jvs(26) * x(27) - jvs(27) * x(29) - &
1187	jvs(28) * x(30) - jvs(29)&
1188	&* x(31)) /(jvs(21))
1189	x(7) = (x(7) - jvs(20) * x(27)) /(jvs(19))
1190	x(6) = (x(6) - jvs(17) * x(26) - jvs(18) * x(30)) /(jvs(16))
1191	x(5) = (x(5) - jvs(15) * x(27)) /(jvs(14))
1192	x(4) = (x(4) - jvs(10) * x(14) - jvs(11) * x(26) - jvs(12) * x(27) - jvs(13) * x(30)) /(jvs(9))
1193	x(3) = (x(3) - jvs(7) * x(26) - jvs(8) * x(32)) /(jvs(6))
1194	x(2) = (x(2) - jvs(4) * x(27) - jvs(5) * x(28)) /(jvs(3))
1195	x(1) = (x(1) - jvs(2) * x(25)) /(jvs(1))
1196
1197	END SUBROUTINE kppsolve
1198
1199	SUBROUTINE jac_sp(v, f, rct, jvs)
1200
1201	! V - Concentrations of variable species (local)
1202	REAL(kind=dp):: v(nvar)
1203	! F - Concentrations of fixed species (local)
1204	REAL(kind=dp):: f(nfix)
1205	! RCT - Rate constants (local)
1206	REAL(kind=dp):: rct(nreact)
1207	! JVS - sparse Jacobian of variables
1208	REAL(kind=dp):: jvs(lu_nonzero)
1209
1210
1211	! Local variables
1212	! B - Temporary array
1213	REAL(kind=dp):: b(138)
1214
1215	! B(1) = dA(1)/dV(26)
1216	b(1) = rct(1)
1217	! B(2) = dA(2)/dV(25)
1218	b(2) = rct(2)
1219	! B(3) = dA(3)/dV(25)
1220	b(3) = rct(3)
1221	! B(4) = dA(4)/dV(30)
1222	b(4) = rct(4)
1223	! B(5) = dA(5)/dV(9)
1224	b(5) = rct(5)
1225	! B(6) = dA(6)/dV(2)
1226	b(6) = rct(6)
1227	! B(7) = dA(7)/dV(21)
1228	b(7) = rct(7)
1229	! B(8) = dA(8)/dV(21)
1230	b(8) = rct(8)
1231	! B(9) = dA(9)/dV(24)
1232	b(9) = rct(9)
1233	! B(10) = dA(10)/dV(19)
1234	b(10) = rct(10)
1235	! B(11) = dA(11)/dV(15)
1236	b(11) = rct(11)
1237	! B(12) = dA(12)/dV(29)
1238	b(12) = rct(12)
1239	! B(13) = dA(13)/dV(25)
1240	b(13) = rct(13) * v(31)
1241	! B(14) = dA(13)/dV(31)
1242	b(14) = rct(13) * v(25)
1243	! B(15) = dA(14)/dV(26)
1244	b(15) = rct(14) * v(29)
1245	! B(16) = dA(14)/dV(29)
1246	b(16) = rct(14) * v(26)
1247	! B(17) = dA(15)/dV(26)
1248	b(17) = rct(15) * v(29)
1249	! B(18) = dA(15)/dV(29)
1250	b(18) = rct(15) * v(26)
1251	! B(19) = dA(16)/dV(29)
1252	b(19) = rct(16) * v(31)
1253	! B(20) = dA(16)/dV(31)
1254	b(20) = rct(16) * v(29)
1255	! B(21) = dA(17)/dV(25)
1256	b(21) = rct(17) * v(26)
1257	! B(22) = dA(17)/dV(26)
1258	b(22) = rct(17) * v(25)
1259	! B(23) = dA(18)/dV(1)
1260	b(23) = rct(18)
1261	! B(24) = dA(19)/dV(1)
1262	b(24) = rct(19) * f(1)
1263	! B(26) = dA(20)/dV(25)
1264	b(26) = rct(20) * v(27)
1265	! B(27) = dA(20)/dV(27)
1266	b(27) = rct(20) * v(25)
1267	! B(28) = dA(21)/dV(25)
1268	b(28) = rct(21) * v(28)
1269	! B(29) = dA(21)/dV(28)
1270	b(29) = rct(21) * v(25)
1271	! B(30) = dA(22)/dV(30)
1272	b(30) = rct(22) * v(31)
1273	! B(31) = dA(22)/dV(31)
1274	b(31) = rct(22) * v(30)
1275	! B(32) = dA(23)/dV(26)
1276	b(32) = rct(23) * v(30)
1277	! B(33) = dA(23)/dV(30)
1278	b(33) = rct(23) * v(26)
1279	! B(34) = dA(24)/dV(26)
1280	b(34) = rct(24) * v(30)
1281	! B(35) = dA(24)/dV(30)
1282	b(35) = rct(24) * v(26)
1283	! B(36) = dA(25)/dV(6)
1284	b(36) = rct(25) * f(1)
1285	! B(38) = dA(26)/dV(6)
1286	b(38) = rct(26)
1287	! B(39) = dA(27)/dV(31)
1288	b(39) = rct(27) * 2* v(31)
1289	! B(40) = dA(28)/dV(26)
1290	b(40) = rct(28) * v(31) * f(1)
1291	! B(41) = dA(28)/dV(31)
1292	b(41) = rct(28) * v(26) * f(1)
1293	! B(43) = dA(29)/dV(27)
1294	b(43) = rct(29) * v(31)
1295	! B(44) = dA(29)/dV(31)
1296	b(44) = rct(29) * v(27)
1297	! B(45) = dA(30)/dV(9)
1298	b(45) = rct(30) * v(27)
1299	! B(46) = dA(30)/dV(27)
1300	b(46) = rct(30) * v(9)
1301	! B(47) = dA(31)/dV(9)
1302	b(47) = rct(31) * 2* v(9)
1303	! B(48) = dA(32)/dV(26)
1304	b(48) = rct(32) * v(27)
1305	! B(49) = dA(32)/dV(27)
1306	b(49) = rct(32) * v(26)
1307	! B(50) = dA(33)/dV(12)
1308	b(50) = rct(33) * v(27)
1309	! B(51) = dA(33)/dV(27)
1310	b(51) = rct(33) * v(12)
1311	! B(52) = dA(34)/dV(28)
1312	b(52) = rct(34) * v(31)
1313	! B(53) = dA(34)/dV(31)
1314	b(53) = rct(34) * v(28)
1315	! B(54) = dA(35)/dV(26)
1316	b(54) = rct(35) * v(28)
1317	! B(55) = dA(35)/dV(28)
1318	b(55) = rct(35) * v(26)
1319	! B(56) = dA(36)/dV(10)
1320	b(56) = rct(36)
1321	! B(57) = dA(37)/dV(10)
1322	b(57) = rct(37) * v(27)
1323	! B(58) = dA(37)/dV(27)
1324	b(58) = rct(37) * v(10)
1325	! B(59) = dA(38)/dV(28)
1326	b(59) = rct(38) * 2* v(28)
1327	! B(60) = dA(39)/dV(28)
1328	b(60) = rct(39) * 2* v(28) * f(1)
1329	! B(62) = dA(40)/dV(2)
1330	b(62) = rct(40) * v(27)
1331	! B(63) = dA(40)/dV(27)
1332	b(63) = rct(40) * v(2)
1333	! B(64) = dA(41)/dV(16)
1334	b(64) = rct(41) * v(27)
1335	! B(65) = dA(41)/dV(27)
1336	b(65) = rct(41) * v(16)
1337	! B(66) = dA(42)/dV(21)
1338	b(66) = rct(42) * v(27)
1339	! B(67) = dA(42)/dV(27)
1340	b(67) = rct(42) * v(21)
1341	! B(68) = dA(43)/dV(21)
1342	b(68) = rct(43) * v(29)
1343	! B(69) = dA(43)/dV(29)
1344	b(69) = rct(43) * v(21)
1345	! B(70) = dA(44)/dV(21)
1346	b(70) = rct(44) * v(30)
1347	! B(71) = dA(44)/dV(30)
1348	b(71) = rct(44) * v(21)
1349	! B(72) = dA(45)/dV(24)
1350	b(72) = rct(45) * v(29)
1351	! B(73) = dA(45)/dV(29)
1352	b(73) = rct(45) * v(24)
1353	! B(74) = dA(46)/dV(24)
1354	b(74) = rct(46) * v(27)
1355	! B(75) = dA(46)/dV(27)
1356	b(75) = rct(46) * v(24)
1357	! B(76) = dA(47)/dV(24)
1358	b(76) = rct(47) * v(30)
1359	! B(77) = dA(47)/dV(30)
1360	b(77) = rct(47) * v(24)
1361	! B(78) = dA(48)/dV(31)
1362	b(78) = rct(48) * v(32)
1363	! B(79) = dA(48)/dV(32)
1364	b(79) = rct(48) * v(31)
1365	! B(80) = dA(49)/dV(26)
1366	b(80) = rct(49) * v(32)
1367	! B(81) = dA(49)/dV(32)
1368	b(81) = rct(49) * v(26)
1369	! B(82) = dA(50)/dV(3)
1370	b(82) = rct(50)
1371	! B(83) = dA(51)/dV(32)
1372	b(83) = rct(51) * 2* v(32)
1373	! B(84) = dA(52)/dV(28)
1374	b(84) = rct(52) * v(32)
1375	! B(85) = dA(52)/dV(32)
1376	b(85) = rct(52) * v(28)
1377	! B(86) = dA(53)/dV(27)
1378	b(86) = rct(53)
1379	! B(87) = dA(54)/dV(20)
1380	b(87) = rct(54) * v(27)
1381	! B(88) = dA(54)/dV(27)
1382	b(88) = rct(54) * v(20)
1383	! B(89) = dA(55)/dV(13)
1384	b(89) = rct(55)
1385	! B(90) = dA(56)/dV(13)
1386	b(90) = rct(56)
1387	! B(91) = dA(57)/dV(13)
1388	b(91) = rct(57) * v(26)
1389	! B(92) = dA(57)/dV(26)
1390	b(92) = rct(57) * v(13)
1391	! B(93) = dA(58)/dV(23)
1392	b(93) = rct(58) * v(29)
1393	! B(94) = dA(58)/dV(29)
1394	b(94) = rct(58) * v(23)
1395	! B(95) = dA(59)/dV(23)
1396	b(95) = rct(59) * v(27)
1397	! B(96) = dA(59)/dV(27)
1398	b(96) = rct(59) * v(23)
1399	! B(97) = dA(60)/dV(23)
1400	b(97) = rct(60) * v(25)
1401	! B(98) = dA(60)/dV(25)
1402	b(98) = rct(60) * v(23)
1403	! B(99) = dA(61)/dV(23)
1404	b(99) = rct(61) * v(30)
1405	! B(100) = dA(61)/dV(30)
1406	b(100) = rct(61) * v(23)
1407	! B(101) = dA(62)/dV(17)
1408	b(101) = rct(62) * v(29)
1409	! B(102) = dA(62)/dV(29)
1410	b(102) = rct(62) * v(17)
1411	! B(103) = dA(63)/dV(17)
1412	b(103) = rct(63) * v(27)
1413	! B(104) = dA(63)/dV(27)
1414	b(104) = rct(63) * v(17)
1415	! B(105) = dA(64)/dV(17)
1416	b(105) = rct(64) * v(25)
1417	! B(106) = dA(64)/dV(25)
1418	b(106) = rct(64) * v(17)
1419	! B(107) = dA(65)/dV(5)
1420	b(107) = rct(65) * v(27)
1421	! B(108) = dA(65)/dV(27)
1422	b(108) = rct(65) * v(5)
1423	! B(109) = dA(66)/dV(11)
1424	b(109) = rct(66) * v(31)
1425	! B(110) = dA(66)/dV(31)
1426	b(110) = rct(66) * v(11)
1427	! B(111) = dA(67)/dV(11)
1428	b(111) = rct(67)
1429	! B(112) = dA(68)/dV(14)
1430	b(112) = rct(68) * v(27)
1431	! B(113) = dA(68)/dV(27)
1432	b(113) = rct(68) * v(14)
1433	! B(114) = dA(69)/dV(14)
1434	b(114) = rct(69) * v(30)
1435	! B(115) = dA(69)/dV(30)
1436	b(115) = rct(69) * v(14)
1437	! B(116) = dA(70)/dV(4)
1438	b(116) = rct(70) * v(26)
1439	! B(117) = dA(70)/dV(26)
1440	b(117) = rct(70) * v(4)
1441	! B(118) = dA(71)/dV(7)
1442	b(118) = rct(71) * v(27)
1443	! B(119) = dA(71)/dV(27)
1444	b(119) = rct(71) * v(7)
1445	! B(120) = dA(72)/dV(19)
1446	b(120) = rct(72) * v(27)
1447	! B(121) = dA(72)/dV(27)
1448	b(121) = rct(72) * v(19)
1449	! B(122) = dA(73)/dV(19)
1450	b(122) = rct(73) * v(25)
1451	! B(123) = dA(73)/dV(25)
1452	b(123) = rct(73) * v(19)
1453	! B(124) = dA(74)/dV(15)
1454	b(124) = rct(74) * v(27)
1455	! B(125) = dA(74)/dV(27)
1456	b(125) = rct(74) * v(15)
1457	! B(126) = dA(75)/dV(22)
1458	b(126) = rct(75) * v(29)
1459	! B(127) = dA(75)/dV(29)
1460	b(127) = rct(75) * v(22)
1461	! B(128) = dA(76)/dV(22)
1462	b(128) = rct(76) * v(27)
1463	! B(129) = dA(76)/dV(27)
1464	b(129) = rct(76) * v(22)
1465	! B(130) = dA(77)/dV(22)
1466	b(130) = rct(77) * v(25)
1467	! B(131) = dA(77)/dV(25)
1468	b(131) = rct(77) * v(22)
1469	! B(132) = dA(78)/dV(22)
1470	b(132) = rct(78) * v(30)
1471	! B(133) = dA(78)/dV(30)
1472	b(133) = rct(78) * v(22)
1473	! B(134) = dA(79)/dV(18)
1474	b(134) = rct(79) * v(31)
1475	! B(135) = dA(79)/dV(31)
1476	b(135) = rct(79) * v(18)
1477	! B(136) = dA(80)/dV(18)
1478	b(136) = rct(80) * 2* v(18)
1479	! B(137) = dA(81)/dV(8)
1480	b(137) = rct(81) * v(31)
1481	! B(138) = dA(81)/dV(31)
1482	b(138) = rct(81) * v(8)
1483
1484	! Construct the Jacobian terms from B's
1485	! JVS(1) = Jac_FULL(1,1)
1486	jvs(1) = - b(23) - b(24)
1487	! JVS(2) = Jac_FULL(1,25)
1488	jvs(2) = b(3)
1489	! JVS(3) = Jac_FULL(2,2)
1490	jvs(3) = - b(6) - b(62)
1491	! JVS(4) = Jac_FULL(2,27)
1492	jvs(4) = - b(63)
1493	! JVS(5) = Jac_FULL(2,28)
1494	jvs(5) = b(59) + b(60)
1495	! JVS(6) = Jac_FULL(3,3)
1496	jvs(6) = - b(82)
1497	! JVS(7) = Jac_FULL(3,26)
1498	jvs(7) = b(80)
1499	! JVS(8) = Jac_FULL(3,32)
1500	jvs(8) = b(81)
1501	! JVS(9) = Jac_FULL(4,4)
1502	jvs(9) = - b(116)
1503	! JVS(10) = Jac_FULL(4,14)
1504	jvs(10) = 0.4* b(112) + b(114)
1505	! JVS(11) = Jac_FULL(4,26)
1506	jvs(11) = - b(117)
1507	! JVS(12) = Jac_FULL(4,27)
1508	jvs(12) = 0.4* b(113)
1509	! JVS(13) = Jac_FULL(4,30)
1510	jvs(13) = b(115)
1511	! JVS(14) = Jac_FULL(5,5)
1512	jvs(14) = - b(107)
1513	! JVS(15) = Jac_FULL(5,27)
1514	jvs(15) = - b(108)
1515	! JVS(16) = Jac_FULL(6,6)
1516	jvs(16) = - b(36) - b(38)
1517	! JVS(17) = Jac_FULL(6,26)
1518	jvs(17) = b(34)
1519	! JVS(18) = Jac_FULL(6,30)
1520	jvs(18) = b(35)
1521	! JVS(19) = Jac_FULL(7,7)
1522	jvs(19) = - b(118)
1523	! JVS(20) = Jac_FULL(7,27)
1524	jvs(20) = - b(119)
1525	! JVS(21) = Jac_FULL(8,8)
1526	jvs(21) = - b(137)
1527	! JVS(22) = Jac_FULL(8,13)
1528	jvs(22) = 0.04* b(89)
1529	! JVS(23) = Jac_FULL(8,20)
1530	jvs(23) = 0.13* b(87)
1531	! JVS(24) = Jac_FULL(8,22)
1532	jvs(24) = 0.13* b(128) + b(132)
1533	! JVS(25) = Jac_FULL(8,23)
1534	jvs(25) = 0.02* b(93) + 0.09* b(99)
1535	! JVS(26) = Jac_FULL(8,27)
1536	jvs(26) = 0.13* b(88) + 0.13* b(129)
1537	! JVS(27) = Jac_FULL(8,29)
1538	jvs(27) = 0.02* b(94)
1539	! JVS(28) = Jac_FULL(8,30)
1540	jvs(28) = 0.09* b(100) + b(133)
1541	! JVS(29) = Jac_FULL(8,31)
1542	jvs(29) = - b(138)
1543	! JVS(30) = Jac_FULL(9,9)
1544	jvs(30) = - b(5) - b(45) - 2* b(47)
1545	! JVS(31) = Jac_FULL(9,26)
1546	jvs(31) = 2* b(40)
1547	! JVS(32) = Jac_FULL(9,27)
1548	jvs(32) = b(43) - b(46)
1549	! JVS(33) = Jac_FULL(9,31)
1550	jvs(33) = 2* b(41) + b(44)
1551	! JVS(34) = Jac_FULL(10,10)
1552	jvs(34) = - b(56) - b(57)
1553	! JVS(35) = Jac_FULL(10,26)
1554	jvs(35) = b(54)
1555	! JVS(36) = Jac_FULL(10,27)
1556	jvs(36) = - b(58)
1557	! JVS(37) = Jac_FULL(10,28)
1558	jvs(37) = b(55)
1559	! JVS(38) = Jac_FULL(11,5)
1560	jvs(38) = 0.56* b(107)
1561	! JVS(39) = Jac_FULL(11,7)
1562	jvs(39) = 0.3* b(118)
1563	! JVS(40) = Jac_FULL(11,11)
1564	jvs(40) = - b(109) - b(111)
1565	! JVS(41) = Jac_FULL(11,27)
1566	jvs(41) = 0.56* b(108) + 0.3* b(119)
1567	! JVS(42) = Jac_FULL(11,31)
1568	jvs(42) = - b(110)
1569	! JVS(43) = Jac_FULL(12,6)
1570	jvs(43) = 2* b(36)
1571	! JVS(44) = Jac_FULL(12,12)
1572	jvs(44) = - b(50)
1573	! JVS(45) = Jac_FULL(12,14)
1574	jvs(45) = b(114)
1575	! JVS(46) = Jac_FULL(12,21)
1576	jvs(46) = b(70)
1577	! JVS(47) = Jac_FULL(12,24)
1578	jvs(47) = b(76)
1579	! JVS(48) = Jac_FULL(12,26)
1580	jvs(48) = b(48)
1581	! JVS(49) = Jac_FULL(12,27)
1582	jvs(49) = b(49) - b(51)
1583	! JVS(50) = Jac_FULL(12,30)
1584	jvs(50) = b(71) + b(77) + b(115)
1585	! JVS(51) = Jac_FULL(13,13)
1586	jvs(51) = - 0.98* b(89) - b(90) - b(91)
1587	! JVS(52) = Jac_FULL(13,20)
1588	jvs(52) = 0.76* b(87)
1589	! JVS(53) = Jac_FULL(13,26)
1590	jvs(53) = - b(92)
1591	! JVS(54) = Jac_FULL(13,27)
1592	jvs(54) = 0.76* b(88)
1593	! JVS(55) = Jac_FULL(14,5)
1594	jvs(55) = 0.36* b(107)
1595	! JVS(56) = Jac_FULL(14,7)
1596	jvs(56) = 0.2* b(118)
1597	! JVS(57) = Jac_FULL(14,11)
1598	jvs(57) = b(111)
1599	! JVS(58) = Jac_FULL(14,14)
1600	jvs(58) = - b(112) - b(114)
1601	! JVS(59) = Jac_FULL(14,27)
1602	jvs(59) = 0.36* b(108) - b(113) + 0.2* b(119)
1603	! JVS(60) = Jac_FULL(14,30)
1604	jvs(60) = - b(115)
1605	! JVS(61) = Jac_FULL(14,31)
1606	jvs(61) = 0
1607	! JVS(62) = Jac_FULL(15,7)
1608	jvs(62) = 0.8* b(118)
1609	! JVS(63) = Jac_FULL(15,15)
1610	jvs(63) = - b(11) - b(124)
1611	! JVS(64) = Jac_FULL(15,19)
1612	jvs(64) = 0.2* b(122)
1613	! JVS(65) = Jac_FULL(15,22)
1614	jvs(65) = 0.4* b(128) + 0.2* b(130)
1615	! JVS(66) = Jac_FULL(15,25)
1616	jvs(66) = 0.2* b(123) + 0.2* b(131)
1617	! JVS(67) = Jac_FULL(15,27)
1618	jvs(67) = 0.8* b(119) - b(125) + 0.4* b(129)
1619	! JVS(68) = Jac_FULL(16,15)
1620	jvs(68) = b(11)
1621	! JVS(69) = Jac_FULL(16,16)
1622	jvs(69) = - b(64)
1623	! JVS(70) = Jac_FULL(16,17)
1624	jvs(70) = b(101) + 0.42* b(105)
1625	! JVS(71) = Jac_FULL(16,19)
1626	jvs(71) = b(10) + 2* b(120) + 0.69* b(122)
1627	! JVS(72) = Jac_FULL(16,21)
1628	jvs(72) = b(7) + b(8) + b(66) + b(68) + b(70)
1629	! JVS(73) = Jac_FULL(16,22)
1630	jvs(73) = 0.5* b(126) + 0.06* b(130)
1631	! JVS(74) = Jac_FULL(16,23)
1632	jvs(74) = 0.3* b(93) + 0.33* b(97)
1633	! JVS(75) = Jac_FULL(16,24)
1634	jvs(75) = b(9)
1635	! JVS(76) = Jac_FULL(16,25)
1636	jvs(76) = 0.33* b(98) + 0.42* b(106) + 0.69* b(123) + 0.06* b(131)
1637	! JVS(77) = Jac_FULL(16,27)
1638	jvs(77) = - b(65) + b(67) + 2* b(121)
1639	! JVS(78) = Jac_FULL(16,29)
1640	jvs(78) = b(69) + 0.3* b(94) + b(102) + 0.5* b(127)
1641	! JVS(79) = Jac_FULL(16,30)
1642	jvs(79) = b(71)
1643	! JVS(80) = Jac_FULL(17,17)
1644	jvs(80) = - b(101) - b(103) - b(105)
1645	! JVS(81) = Jac_FULL(17,22)
1646	jvs(81) = 0.45* b(126) + b(128) + 0.55* b(130)
1647	! JVS(82) = Jac_FULL(17,25)
1648	jvs(82) = - b(106) + 0.55* b(131)
1649	! JVS(83) = Jac_FULL(17,27)
1650	jvs(83) = - b(104) + b(129)
1651	! JVS(84) = Jac_FULL(17,29)
1652	jvs(84) = - b(102) + 0.45* b(127)
1653	! JVS(85) = Jac_FULL(18,5)
1654	jvs(85) = 0.08* b(107)
1655	! JVS(86) = Jac_FULL(18,7)
1656	jvs(86) = 0.5* b(118)
1657	! JVS(87) = Jac_FULL(18,13)
1658	jvs(87) = 0.96* b(89)
1659	! JVS(88) = Jac_FULL(18,14)
1660	jvs(88) = 0.6* b(112)
1661	! JVS(89) = Jac_FULL(18,15)
1662	jvs(89) = b(124)
1663	! JVS(90) = Jac_FULL(18,17)
1664	jvs(90) = 0.7* b(101) + b(103)
1665	! JVS(91) = Jac_FULL(18,18)
1666	jvs(91) = - b(134) - 2* b(136)
1667	! JVS(92) = Jac_FULL(18,19)
1668	jvs(92) = b(120) + 0.03* b(122)
1669	! JVS(93) = Jac_FULL(18,20)
1670	jvs(93) = 0.87* b(87)
1671	! JVS(94) = Jac_FULL(18,22)
1672	jvs(94) = 0.5* b(126) + b(128)
1673	! JVS(95) = Jac_FULL(18,23)
1674	jvs(95) = 0.28* b(93) + b(95) + 0.22* b(97) + 0.91* b(99)
1675	! JVS(96) = Jac_FULL(18,24)
1676	jvs(96) = b(9)
1677	! JVS(97) = Jac_FULL(18,25)
1678	jvs(97) = 0.22* b(98) + 0.03* b(123)
1679	! JVS(98) = Jac_FULL(18,26)
1680	jvs(98) = 0
1681	! JVS(99) = Jac_FULL(18,27)
1682	jvs(99) = b(86) + 0.87* b(88) + b(96) + b(104) + 0.08* b(108) + 0.6* b(113) + 0.5* b(119) + b(121) + b(125) + b(129)
1683	! JVS(100) = Jac_FULL(18,28)
1684	jvs(100) = 0.79* b(84)
1685	! JVS(101) = Jac_FULL(18,29)
1686	jvs(101) = 0.28* b(94) + 0.7* b(102) + 0.5* b(127)
1687	! JVS(102) = Jac_FULL(18,30)
1688	jvs(102) = 0.91* b(100)
1689	! JVS(103) = Jac_FULL(18,31)
1690	jvs(103) = b(78) - b(135)
1691	! JVS(104) = Jac_FULL(18,32)
1692	jvs(104) = b(79) + 2* b(83) + 0.79* b(85)
1693	! JVS(105) = Jac_FULL(19,11)
1694	jvs(105) = 0.9* b(109)
1695	! JVS(106) = Jac_FULL(19,14)
1696	jvs(106) = 0.3* b(112)
1697	! JVS(107) = Jac_FULL(19,19)
1698	jvs(107) = - b(10) - b(120) - b(122)
1699	! JVS(108) = Jac_FULL(19,25)
1700	jvs(108) = - b(123)
1701	! JVS(109) = Jac_FULL(19,27)
1702	jvs(109) = 0.3* b(113) - b(121)
1703	! JVS(110) = Jac_FULL(19,30)
1704	jvs(110) = 0
1705	! JVS(111) = Jac_FULL(19,31)
1706	jvs(111) = 0.9* b(110)
1707	! JVS(112) = Jac_FULL(20,7)
1708	jvs(112) = 1.1* b(118)
1709	! JVS(113) = Jac_FULL(20,13)
1710	jvs(113) = - 2.1* b(89)
1711	! JVS(114) = Jac_FULL(20,20)
1712	jvs(114) = - 1.11* b(87)
1713	! JVS(115) = Jac_FULL(20,22)
1714	jvs(115) = 0.9* b(126) + 0.1* b(130)
1715	! JVS(116) = Jac_FULL(20,23)
1716	jvs(116) = 0.22* b(93) - b(95) - b(97) - b(99)
1717	! JVS(117) = Jac_FULL(20,25)
1718	jvs(117) = - b(98) + 0.1* b(131)
1719	! JVS(118) = Jac_FULL(20,26)
1720	jvs(118) = 0
1721	! JVS(119) = Jac_FULL(20,27)
1722	jvs(119) = - 1.11* b(88) - b(96) + 1.1* b(119)
1723	! JVS(120) = Jac_FULL(20,29)
1724	jvs(120) = 0.22* b(94) + 0.9* b(127)
1725	! JVS(121) = Jac_FULL(20,30)
1726	jvs(121) = - b(100)
1727	! JVS(122) = Jac_FULL(21,17)
1728	jvs(122) = b(101) + 1.56* b(103) + b(105)
1729	! JVS(123) = Jac_FULL(21,19)
1730	jvs(123) = b(120) + 0.7* b(122)
1731	! JVS(124) = Jac_FULL(21,21)
1732	jvs(124) = - b(7) - b(8) - b(66) - b(68) - b(70)
1733	! JVS(125) = Jac_FULL(21,22)
1734	jvs(125) = b(128) + b(130)
1735	! JVS(126) = Jac_FULL(21,23)
1736	jvs(126) = 0.2* b(93) + b(95) + 0.74* b(97) + b(99)
1737	! JVS(127) = Jac_FULL(21,24)
1738	jvs(127) = b(9)
1739	! JVS(128) = Jac_FULL(21,25)
1740	jvs(128) = 0.74* b(98) + b(106) + 0.7* b(123) + b(131)
1741	! JVS(129) = Jac_FULL(21,27)
1742	jvs(129) = - b(67) + b(86) + b(96) + 1.56* b(104) + b(121) + b(129)
1743	! JVS(130) = Jac_FULL(21,28)
1744	jvs(130) = 0.79* b(84)
1745	! JVS(131) = Jac_FULL(21,29)
1746	jvs(131) = - b(69) + 0.2* b(94) + b(102)
1747	! JVS(132) = Jac_FULL(21,30)
1748	jvs(132) = - b(71) + b(100)
1749	! JVS(133) = Jac_FULL(21,31)
1750	jvs(133) = b(78)
1751	! JVS(134) = Jac_FULL(21,32)
1752	jvs(134) = b(79) + 2* b(83) + 0.79* b(85)
1753	! JVS(135) = Jac_FULL(22,22)
1754	jvs(135) = - b(126) - b(128) - b(130) - b(132)
1755	! JVS(136) = Jac_FULL(22,25)
1756	jvs(136) = - b(131)
1757	! JVS(137) = Jac_FULL(22,27)
1758	jvs(137) = - b(129)
1759	! JVS(138) = Jac_FULL(22,29)
1760	jvs(138) = - b(127)
1761	! JVS(139) = Jac_FULL(22,30)
1762	jvs(139) = - b(133)
1763	! JVS(140) = Jac_FULL(23,22)
1764	jvs(140) = 0.55* b(126)
1765	! JVS(141) = Jac_FULL(23,23)
1766	jvs(141) = - b(93) - b(95) - b(97) - b(99)
1767	! JVS(142) = Jac_FULL(23,25)
1768	jvs(142) = - b(98)
1769	! JVS(143) = Jac_FULL(23,27)
1770	jvs(143) = - b(96)
1771	! JVS(144) = Jac_FULL(23,29)
1772	jvs(144) = - b(94) + 0.55* b(127)
1773	! JVS(145) = Jac_FULL(23,30)
1774	jvs(145) = - b(100)
1775	! JVS(146) = Jac_FULL(24,13)
1776	jvs(146) = 1.1* b(89)
1777	! JVS(147) = Jac_FULL(24,17)
1778	jvs(147) = 0.22* b(103)
1779	! JVS(148) = Jac_FULL(24,19)
1780	jvs(148) = 0.03* b(122)
1781	! JVS(149) = Jac_FULL(24,20)
1782	jvs(149) = 0.11* b(87)
1783	! JVS(150) = Jac_FULL(24,22)
1784	jvs(150) = 0.8* b(126) + 0.2* b(128) + 0.4* b(130)
1785	! JVS(151) = Jac_FULL(24,23)
1786	jvs(151) = 0.63* b(93) + b(95) + 0.5* b(97) + b(99)
1787	! JVS(152) = Jac_FULL(24,24)
1788	jvs(152) = - b(9) - b(72) - b(74) - b(76)
1789	! JVS(153) = Jac_FULL(24,25)
1790	jvs(153) = 0.5* b(98) + 0.03* b(123) + 0.4* b(131)
1791	! JVS(154) = Jac_FULL(24,26)
1792	jvs(154) = 0
1793	! JVS(155) = Jac_FULL(24,27)
1794	jvs(155) = - b(75) + 0.11* b(88) + b(96) + 0.22* b(104) + 0.2* b(129)
1795	! JVS(156) = Jac_FULL(24,29)
1796	jvs(156) = - b(73) + 0.63* b(94) + 0.8* b(127)
1797	! JVS(157) = Jac_FULL(24,30)
1798	jvs(157) = - b(77) + b(100)
1799	! JVS(158) = Jac_FULL(24,31)
1800	jvs(158) = 0
1801	! JVS(159) = Jac_FULL(25,17)
1802	jvs(159) = - b(105)
1803	! JVS(160) = Jac_FULL(25,19)
1804	jvs(160) = - b(122)
1805	! JVS(161) = Jac_FULL(25,22)
1806	jvs(161) = - b(130)
1807	! JVS(162) = Jac_FULL(25,23)
1808	jvs(162) = - b(97)
1809	! JVS(163) = Jac_FULL(25,25)
1810	jvs(163) = - b(2) - b(3) - b(13) - b(21) - b(26) - b(28) - b(98) - b(106) - b(123) - b(131)
1811	! JVS(164) = Jac_FULL(25,26)
1812	jvs(164) = - b(22)
1813	! JVS(165) = Jac_FULL(25,27)
1814	jvs(165) = - b(27)
1815	! JVS(166) = Jac_FULL(25,28)
1816	jvs(166) = - b(29)
1817	! JVS(167) = Jac_FULL(25,29)
1818	jvs(167) = b(12)
1819	! JVS(168) = Jac_FULL(25,30)
1820	jvs(168) = 0
1821	! JVS(169) = Jac_FULL(25,31)
1822	jvs(169) = - b(14)
1823	! JVS(170) = Jac_FULL(26,3)
1824	jvs(170) = b(82)
1825	! JVS(171) = Jac_FULL(26,4)
1826	jvs(171) = - b(116)
1827	! JVS(172) = Jac_FULL(26,6)
1828	jvs(172) = b(38)
1829	! JVS(173) = Jac_FULL(26,9)
1830	jvs(173) = b(45) + b(47)
1831	! JVS(174) = Jac_FULL(26,10)
1832	jvs(174) = b(56) + b(57)
1833	! JVS(175) = Jac_FULL(26,11)
1834	jvs(175) = 0.9* b(109)
1835	! JVS(176) = Jac_FULL(26,13)
1836	jvs(176) = - b(91)
1837	! JVS(177) = Jac_FULL(26,14)
1838	jvs(177) = 0
1839	! JVS(178) = Jac_FULL(26,18)
1840	jvs(178) = b(134)
1841	! JVS(179) = Jac_FULL(26,19)
1842	jvs(179) = 0
1843	! JVS(180) = Jac_FULL(26,20)
1844	jvs(180) = 0
1845	! JVS(181) = Jac_FULL(26,22)
1846	jvs(181) = 0
1847	! JVS(182) = Jac_FULL(26,23)
1848	jvs(182) = b(99)
1849	! JVS(183) = Jac_FULL(26,24)
1850	jvs(183) = 0
1851	! JVS(184) = Jac_FULL(26,25)
1852	jvs(184) = b(13) - b(21)
1853	! JVS(185) = Jac_FULL(26,26)
1854	jvs(185) = - b(1) - b(15) - b(17) - b(22) - b(34) - b(40) - b(48) - b(54) - b(80) - b(92) - b(117)
1855	! JVS(186) = Jac_FULL(26,27)
1856	jvs(186) = b(46) - b(49) + b(58)
1857	! JVS(187) = Jac_FULL(26,28)
1858	jvs(187) = b(52) - b(55)
1859	! JVS(188) = Jac_FULL(26,29)
1860	jvs(188) = - b(16) - b(18) + b(19)
1861	! JVS(189) = Jac_FULL(26,30)
1862	jvs(189) = 0.89* b(4) + 2* b(30) - b(35) + b(100)
1863	! JVS(190) = Jac_FULL(26,31)
1864	jvs(190) = b(14) + b(20) + 2* b(31) + 2* b(39) - b(41) + b(53) + b(78) + 0.9* b(110) + b(135)
1865	! JVS(191) = Jac_FULL(26,32)
1866	jvs(191) = b(79) - b(81)
1867	! JVS(192) = Jac_FULL(27,1)
1868	jvs(192) = 2* b(24)
1869	! JVS(193) = Jac_FULL(27,2)
1870	jvs(193) = 2* b(6) - b(62)
1871	! JVS(194) = Jac_FULL(27,5)
1872	jvs(194) = - b(107)
1873	! JVS(195) = Jac_FULL(27,7)
1874	jvs(195) = - b(118)
1875	! JVS(196) = Jac_FULL(27,9)
1876	jvs(196) = b(5) - b(45)
1877	! JVS(197) = Jac_FULL(27,10)
1878	jvs(197) = - b(57)
1879	! JVS(198) = Jac_FULL(27,12)
1880	jvs(198) = - b(50)
1881	! JVS(199) = Jac_FULL(27,14)
1882	jvs(199) = - b(112)
1883	! JVS(200) = Jac_FULL(27,15)
1884	jvs(200) = - b(124)
1885	! JVS(201) = Jac_FULL(27,16)
1886	jvs(201) = - b(64)
1887	! JVS(202) = Jac_FULL(27,17)
1888	jvs(202) = 0.3* b(101) - b(103)
1889	! JVS(203) = Jac_FULL(27,19)
1890	jvs(203) = - b(120) + 0.08* b(122)
1891	! JVS(204) = Jac_FULL(27,20)
1892	jvs(204) = - b(87)
1893	! JVS(205) = Jac_FULL(27,21)
1894	jvs(205) = - b(66) + b(68)
1895	! JVS(206) = Jac_FULL(27,22)
1896	jvs(206) = - b(128) + 0.1* b(130)
1897	! JVS(207) = Jac_FULL(27,23)
1898	jvs(207) = 0.2* b(93) - b(95) + 0.1* b(97)
1899	! JVS(208) = Jac_FULL(27,24)
1900	jvs(208) = b(72) - b(74)
1901	! JVS(209) = Jac_FULL(27,25)
1902	jvs(209) = - b(26) + b(28) + 0.1* b(98) + 0.08* b(123) + 0.1* b(131)
1903	! JVS(210) = Jac_FULL(27,26)
1904	jvs(210) = - b(48)
1905	! JVS(211) = Jac_FULL(27,27)
1906	jvs(211) = - b(27) - b(43) - b(46) - b(49) - b(51) - b(58) - b(63) - b(65) - b(67) - b(75) - b(86) - b(88) - b(96) - b(104) - &
1907	b(108) - b(113) - b(119)&
1908	&- b(121) - b(125) - b(129)
1909	! JVS(212) = Jac_FULL(27,28)
1910	jvs(212) = b(29) + b(52) + 0.79* b(84)
1911	! JVS(213) = Jac_FULL(27,29)
1912	jvs(213) = b(69) + b(73) + 0.2* b(94) + 0.3* b(102)
1913	! JVS(214) = Jac_FULL(27,30)
1914	jvs(214) = 0
1915	! JVS(215) = Jac_FULL(27,31)
1916	jvs(215) = - b(44) + b(53)
1917	! JVS(216) = Jac_FULL(27,32)
1918	jvs(216) = 0.79* b(85)
1919	! JVS(217) = Jac_FULL(28,2)
1920	jvs(217) = b(62)
1921	! JVS(218) = Jac_FULL(28,5)
1922	jvs(218) = 0.44* b(107)
1923	! JVS(219) = Jac_FULL(28,7)
1924	jvs(219) = 0.7* b(118)
1925	! JVS(220) = Jac_FULL(28,10)
1926	jvs(220) = b(56)
1927	! JVS(221) = Jac_FULL(28,11)
1928	jvs(221) = 0.9* b(109) + b(111)
1929	! JVS(222) = Jac_FULL(28,13)
1930	jvs(222) = 0.94* b(89) + b(90)
1931	! JVS(223) = Jac_FULL(28,14)
1932	jvs(223) = 0.6* b(112)
1933	! JVS(224) = Jac_FULL(28,15)
1934	jvs(224) = b(11)
1935	! JVS(225) = Jac_FULL(28,16)
1936	jvs(225) = b(64)
1937	! JVS(226) = Jac_FULL(28,17)
1938	jvs(226) = 1.7* b(101) + b(103) + 0.12* b(105)
1939	! JVS(227) = Jac_FULL(28,19)
1940	jvs(227) = b(10) + 2* b(120) + 0.76* b(122)
1941	! JVS(228) = Jac_FULL(28,20)
1942	jvs(228) = 0.11* b(87)
1943	! JVS(229) = Jac_FULL(28,21)
1944	jvs(229) = 2* b(7) + b(66) + b(68) + b(70)
1945	! JVS(230) = Jac_FULL(28,22)
1946	jvs(230) = 0.6* b(126) + 0.67* b(128) + 0.44* b(130)
1947	! JVS(231) = Jac_FULL(28,23)
1948	jvs(231) = 0.38* b(93) + b(95) + 0.44* b(97)
1949	! JVS(232) = Jac_FULL(28,24)
1950	jvs(232) = 2* b(9)
1951	! JVS(233) = Jac_FULL(28,25)
1952	jvs(233) = b(26) - b(28) + 0.44* b(98) + 0.12* b(106) + 0.76* b(123) + 0.44* b(131)
1953	! JVS(234) = Jac_FULL(28,26)
1954	jvs(234) = - b(54)
1955	! JVS(235) = Jac_FULL(28,27)
1956	jvs(235) = b(27) + b(63) + b(65) + b(67) + b(86) + 0.11* b(88) + b(96) + b(104) + 0.44* b(108) + 0.6* b(113) + 0.7* b(119) + 2* &
1957	b(121) + 0.67&
1958	&* b(129)
1959	! JVS(236) = Jac_FULL(28,28)
1960	jvs(236) = - b(29) - b(52) - b(55) - 2* b(59) - 2* b(60) - 0.21* b(84)
1961	! JVS(237) = Jac_FULL(28,29)
1962	jvs(237) = b(69) + 0.38* b(94) + 1.7* b(102) + 0.6* b(127)
1963	! JVS(238) = Jac_FULL(28,30)
1964	jvs(238) = b(71)
1965	! JVS(239) = Jac_FULL(28,31)
1966	jvs(239) = - b(53) + b(78) + 0.9* b(110)
1967	! JVS(240) = Jac_FULL(28,32)
1968	jvs(240) = b(79) + 2* b(83) - 0.21* b(85)
1969	! JVS(241) = Jac_FULL(29,1)
1970	jvs(241) = b(23)
1971	! JVS(242) = Jac_FULL(29,17)
1972	jvs(242) = - b(101)
1973	! JVS(243) = Jac_FULL(29,21)
1974	jvs(243) = - b(68)
1975	! JVS(244) = Jac_FULL(29,22)
1976	jvs(244) = - b(126)
1977	! JVS(245) = Jac_FULL(29,23)
1978	jvs(245) = - b(93)
1979	! JVS(246) = Jac_FULL(29,24)
1980	jvs(246) = - b(72)
1981	! JVS(247) = Jac_FULL(29,25)
1982	jvs(247) = b(2)
1983	! JVS(248) = Jac_FULL(29,26)
1984	jvs(248) = b(1) - b(15) - b(17)
1985	! JVS(249) = Jac_FULL(29,27)
1986	jvs(249) = 0
1987	! JVS(250) = Jac_FULL(29,28)
1988	jvs(250) = 0
1989	! JVS(251) = Jac_FULL(29,29)
1990	jvs(251) = - b(12) - b(16) - b(18) - b(19) - b(69) - b(73) - b(94) - b(102) - b(127)
1991	! JVS(252) = Jac_FULL(29,30)
1992	jvs(252) = 0.89* b(4)
1993	! JVS(253) = Jac_FULL(29,31)
1994	jvs(253) = - b(20)
1995	! JVS(254) = Jac_FULL(29,32)
1996	jvs(254) = 0
1997	! JVS(255) = Jac_FULL(30,6)
1998	jvs(255) = b(38)
1999	! JVS(256) = Jac_FULL(30,12)
2000	jvs(256) = b(50)
2001	! JVS(257) = Jac_FULL(30,14)
2002	jvs(257) = - b(114)
2003	! JVS(258) = Jac_FULL(30,21)
2004	jvs(258) = - b(70)
2005	! JVS(259) = Jac_FULL(30,22)
2006	jvs(259) = - b(132)
2007	! JVS(260) = Jac_FULL(30,23)
2008	jvs(260) = - b(99)
2009	! JVS(261) = Jac_FULL(30,24)
2010	jvs(261) = - b(76)
2011	! JVS(262) = Jac_FULL(30,25)
2012	jvs(262) = b(21)
2013	! JVS(263) = Jac_FULL(30,26)
2014	jvs(263) = b(17) + b(22) - b(32) - b(34)
2015	! JVS(264) = Jac_FULL(30,27)
2016	jvs(264) = b(51)
2017	! JVS(265) = Jac_FULL(30,28)
2018	jvs(265) = 0
2019	! JVS(266) = Jac_FULL(30,29)
2020	jvs(266) = b(18)
2021	! JVS(267) = Jac_FULL(30,30)
2022	jvs(267) = - b(4) - b(30) - b(33) - b(35) - b(71) - b(77) - b(100) - b(115) - b(133)
2023	! JVS(268) = Jac_FULL(30,31)
2024	jvs(268) = - b(31)
2025	! JVS(269) = Jac_FULL(30,32)
2026	jvs(269) = 0
2027	! JVS(270) = Jac_FULL(31,8)
2028	jvs(270) = - b(137)
2029	! JVS(271) = Jac_FULL(31,9)
2030	jvs(271) = b(5) + b(47)
2031	! JVS(272) = Jac_FULL(31,11)
2032	jvs(272) = - b(109)
2033	! JVS(273) = Jac_FULL(31,13)
2034	jvs(273) = 0
2035	! JVS(274) = Jac_FULL(31,18)
2036	jvs(274) = - b(134)
2037	! JVS(275) = Jac_FULL(31,19)
2038	jvs(275) = 0
2039	! JVS(276) = Jac_FULL(31,20)
2040	jvs(276) = 0
2041	! JVS(277) = Jac_FULL(31,22)
2042	jvs(277) = 0
2043	! JVS(278) = Jac_FULL(31,23)
2044	jvs(278) = 0
2045	! JVS(279) = Jac_FULL(31,24)
2046	jvs(279) = 0
2047	! JVS(280) = Jac_FULL(31,25)
2048	jvs(280) = - b(13)
2049	! JVS(281) = Jac_FULL(31,26)
2050	jvs(281) = b(1) + b(15) + b(32) - b(40)
2051	! JVS(282) = Jac_FULL(31,27)
2052	jvs(282) = - b(43)
2053	! JVS(283) = Jac_FULL(31,28)
2054	jvs(283) = - b(52)
2055	! JVS(284) = Jac_FULL(31,29)
2056	jvs(284) = b(16) - b(19)
2057	! JVS(285) = Jac_FULL(31,30)
2058	jvs(285) = 0.11* b(4) - b(30) + b(33)
2059	! JVS(286) = Jac_FULL(31,31)
2060	jvs(286) = - b(14) - b(20) - b(31) - 2* b(39) - b(41) - b(44) - b(53) - b(78) - b(110) - b(135) - b(138)
2061	! JVS(287) = Jac_FULL(31,32)
2062	jvs(287) = - b(79)
2063	! JVS(288) = Jac_FULL(32,3)
2064	jvs(288) = b(82)
2065	! JVS(289) = Jac_FULL(32,15)
2066	jvs(289) = b(11) + b(124)
2067	! JVS(290) = Jac_FULL(32,19)
2068	jvs(290) = b(10) + b(120) + 0.62* b(122)
2069	! JVS(291) = Jac_FULL(32,22)
2070	jvs(291) = 0.2* b(128)
2071	! JVS(292) = Jac_FULL(32,24)
2072	jvs(292) = b(72) + b(74) + b(76)
2073	! JVS(293) = Jac_FULL(32,25)
2074	jvs(293) = 0.62* b(123)
2075	! JVS(294) = Jac_FULL(32,26)
2076	jvs(294) = - b(80)
2077	! JVS(295) = Jac_FULL(32,27)
2078	jvs(295) = b(75) + b(121) + b(125) + 0.2* b(129)
2079	! JVS(296) = Jac_FULL(32,28)
2080	jvs(296) = - b(84)
2081	! JVS(297) = Jac_FULL(32,29)
2082	jvs(297) = b(73)
2083	! JVS(298) = Jac_FULL(32,30)
2084	jvs(298) = b(77)
2085	! JVS(299) = Jac_FULL(32,31)
2086	jvs(299) = - b(78)
2087	! JVS(300) = Jac_FULL(32,32)
2088	jvs(300) = - b(79) - b(81) - 2* b(83) - b(85)
2089
2090	END SUBROUTINE jac_sp
2091
2092	elemental REAL(kind=dp)FUNCTION k_arr (k_298, tdep, temp)
2093	! arrhenius FUNCTION
2094
2095	REAL, INTENT(IN):: k_298 ! k at t = 298.15k
2096	REAL, INTENT(IN):: tdep ! temperature dependence
2097	REAL(kind=dp), INTENT(IN):: temp ! temperature
2098
2099	intrinsic exp
2100
2101	k_arr = k_298 * exp(tdep* (1._dp/temp- 3.3540e-3_dp))! 1/298.15=3.3540e-3
2102
2103	END FUNCTION k_arr
2104
2105	SUBROUTINE update_rconst()
2106	INTEGER :: k
2107
2108	k = is
2109
2110	! Begin INLINED RCONST
2111
2112
2113	! End INLINED RCONST
2114
2115	rconst(1) = (phot(j_no2))
2116	rconst(2) = (phot(j_o33p))
2117	rconst(3) = (phot(j_o31d))
2118	rconst(4) = (phot(j_no3o) + phot(j_no3o2))
2119	rconst(5) = (phot(j_hono))
2120	rconst(6) = (phot(j_h2o2))
2121	rconst(7) = (phot(j_ch2or))
2122	rconst(8) = (phot(j_ch2om))
2123	rconst(9) = (4.6e-4_dp * phot(j_no2))
2124	rconst(10) = (9.04_dp * phot(j_ch2or))
2125	rconst(11) = (9.64_dp * phot(j_ch2or))
2126	rconst(12) = (arr2(1.4e+3_dp , -1175.0_dp , temp))
2127	rconst(13) = (arr2(1.8e-12_dp , + 1370.0_dp , temp))
2128	rconst(14) = (9.3e-12_dp)
2129	rconst(15) = (arr2(1.6e-13_dp , -687.0_dp , temp))
2130	rconst(16) = (arr2(2.2e-13_dp , -602.0_dp , temp))
2131	rconst(17) = (arr2(1.2e-13_dp , + 2450.0_dp , temp))
2132	rconst(18) = (arr2(1.9e+8_dp , -390.0_dp , temp))
2133	rconst(19) = (2.2e-10_dp)
2134	rconst(20) = (arr2(1.6e-12_dp , + 940.0_dp , temp))
2135	rconst(21) = (arr2(1.4e-14_dp , + 580.0_dp , temp))
2136	rconst(22) = (arr2(1.3e-11_dp , -250.0_dp , temp))
2137	rconst(23) = (arr2(2.5e-14_dp , + 1230.0_dp , temp))
2138	rconst(24) = (arr2(5.3e-13_dp , -256.0_dp , temp))
2139	rconst(25) = (1.3e-21_dp)
2140	rconst(26) = (arr2(3.5e+14_dp , + 10897.0_dp , temp))
2141	rconst(27) = (arr2(1.8e-20_dp , -530.0_dp , temp))
2142	rconst(28) = (4.4e-40_dp)
2143	rconst(29) = (arr2(4.5e-13_dp , -806.0_dp , temp))
2144	rconst(30) = (6.6e-12_dp)
2145	rconst(31) = (1.0e-20_dp)
2146	rconst(32) = (arr2(1.0e-12_dp , -713.0_dp , temp))
2147	rconst(33) = (arr2(5.1e-15_dp , -1000.0_dp , temp))
2148	rconst(34) = (arr2(3.7e-12_dp , -240.0_dp , temp))
2149	rconst(35) = (arr2(1.2e-13_dp , -749.0_dp , temp))
2150	rconst(36) = (arr2(4.8e+13_dp , + 10121.0_dp , temp))
2151	rconst(37) = (arr2(1.3e-12_dp , -380.0_dp , temp))
2152	rconst(38) = (arr2(5.9e-14_dp , -1150.0_dp , temp))
2153	rconst(39) = (arr2(2.2e-38_dp , -5800.0_dp , temp))
2154	rconst(40) = (arr2(3.1e-12_dp , + 187.0_dp , temp))
2155	rconst(41) = (2.2e-13_dp)
2156	rconst(42) = (1.0e-11_dp)
2157	rconst(43) = (arr2(3.0e-11_dp , + 1550.0_dp , temp))
2158	rconst(44) = (6.3e-16_dp)
2159	rconst(45) = (arr2(1.2e-11_dp , + 986.0_dp , temp))
2160	rconst(46) = (arr2(7.0e-12_dp , -250.0_dp , temp))
2161	rconst(47) = (2.5e-15_dp)
2162	rconst(48) = (arr2(5.4e-12_dp , -250.0_dp , temp))
2163	rconst(49) = (arr2(8.0e-20_dp , -5500.0_dp , temp))
2164	rconst(50) = (arr2(9.4e+16_dp , + 14000.0_dp , temp))
2165	rconst(51) = (2.0e-12_dp)
2166	rconst(52) = (6.5e-12_dp)
2167	rconst(53) = (arr2(1.1e+2_dp , + 1710.0_dp , temp))
2168	rconst(54) = (8.1e-13_dp)
2169	rconst(55) = (arr2(1.0e+15_dp , + 8000.0_dp , temp))
2170	rconst(56) = (1.6e+03_dp)
2171	rconst(57) = (1.5e-11_dp)
2172	rconst(58) = (arr2(1.2e-11_dp , + 324.0_dp , temp))
2173	rconst(59) = (arr2(5.2e-12_dp , -504.0_dp , temp))
2174	rconst(60) = (arr2(1.4e-14_dp , + 2105.0_dp , temp))
2175	rconst(61) = (7.7e-15_dp)
2176	rconst(62) = (arr2(1.0e-11_dp , + 792.0_dp , temp))
2177	rconst(63) = (arr2(2.0e-12_dp , -411.0_dp , temp))
2178	rconst(64) = (arr2(1.3e-14_dp , + 2633.0_dp , temp))
2179	rconst(65) = (arr2(2.1e-12_dp , -322.0_dp , temp))
2180	rconst(66) = (8.1e-12_dp)
2181	rconst(67) = (4.20_dp)
2182	rconst(68) = (4.1e-11_dp)
2183	rconst(69) = (2.2e-11_dp)
2184	rconst(70) = (1.4e-11_dp)
2185	rconst(71) = (arr2(1.7e-11_dp , -116.0_dp , temp))
2186	rconst(72) = (3.0e-11_dp)
2187	rconst(73) = (arr2(5.4e-17_dp , + 500.0_dp , temp))
2188	rconst(74) = (1.70e-11_dp)
2189	rconst(75) = (1.80e-11_dp)
2190	rconst(76) = (9.6e-11_dp)
2191	rconst(77) = (1.2e-17_dp)
2192	rconst(78) = (3.2e-13_dp)
2193	rconst(79) = (8.1e-12_dp)
2194	rconst(80) = (arr2(1.7e-14_dp , -1300.0_dp , temp))
2195	rconst(81) = (6.8e-13_dp)
2196
2197	END SUBROUTINE update_rconst
2198
2199	! END FUNCTION ARR2
2200	REAL(kind=dp)FUNCTION arr2( a0, b0, temp)
2201	REAL(kind=dp):: temp
2202	REAL(kind=dp):: a0, b0
2203	arr2 = a0 * exp( - b0 / temp)
2204	END FUNCTION arr2
2205
2206	SUBROUTINE initialize_kpp_ctrl(status)
2207
2208
2209	! i/o
2210	INTEGER, INTENT(OUT):: status
2211
2212	! local
2213	REAL(dp):: tsum
2214	INTEGER :: i
2215
2216	! check fixed time steps
2217	tsum = 0.0_dp
2218	DO i=1, nmaxfixsteps
2219	IF (t_steps(i)< tiny(0.0_dp))exit
2220	tsum = tsum + t_steps(i)
2221	ENDDO
2222
2223	nfsteps = i- 1
2224
2225	l_fixed_step = (nfsteps > 0).and.((tsum - 1.0)< tiny(0.0_dp))
2226
2227	IF (l_vector)THEN
2228	WRITE(,) ' MODE : VECTOR (LENGTH=',VL_DIM,')'
2229	ELSE
2230	WRITE(,) ' MODE : SCALAR'
2231	ENDIF
2232	!
2233	WRITE(,) ' DE-INDEXING MODE :',I_LU_DI
2234	!
2235	WRITE(,) ' ICNTRL : ',icntrl
2236	WRITE(,) ' RCNTRL : ',rcntrl
2237	!
2238	! note: this is ONLY meaningful for vectorized (kp4)rosenbrock- methods
2239	IF (l_vector)THEN
2240	IF (l_fixed_step)THEN
2241	WRITE(,) ' TIME STEPS : FIXED (',t_steps(1:nfsteps),')'
2242	ELSE
2243	WRITE(,) ' TIME STEPS : AUTOMATIC'
2244	ENDIF
2245	ELSE
2246	WRITE(,) ' TIME STEPS : AUTOMATIC '//&
2247	&'(t_steps (CTRL_KPP) ignored in SCALAR MODE)'
2248	ENDIF
2249	! mz_pj_20070531-
2250
2251	status = 0
2252
2253
2254	END SUBROUTINE initialize_kpp_ctrl
2255
2256	SUBROUTINE error_output(c, ierr, pe)
2257
2258
2259	INTEGER, INTENT(IN):: ierr
2260	INTEGER, INTENT(IN):: pe
2261	REAL(dp), DIMENSION(:), INTENT(IN):: c
2262
2263	write(6,*) 'ERROR in chem_gasphase_mod ',ierr,C(1)
2264
2265
2266	END SUBROUTINE error_output
2267
2268	SUBROUTINE wscal(n, alpha, x, incx)
2269	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2270	! constant times a vector: x(1:N) <- Alpha*x(1:N)
2271	! only for incX=incY=1
2272	! after BLAS
2273	! replace this by the function from the optimized BLAS implementation:
2274	! CALL SSCAL(N,Alpha,X,1) or CALL DSCAL(N,Alpha,X,1)
2275	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2276
2277	INTEGER :: i, incx, m, mp1, n
2278	REAL(kind=dp) :: x(n), alpha
2279	REAL(kind=dp), PARAMETER :: zero=0.0_dp, one=1.0_dp
2280
2281	IF (alpha .eq. one)RETURN
2282	IF (n .le. 0)RETURN
2283
2284	m = mod(n, 5)
2285	IF ( m .ne. 0)THEN
2286	IF (alpha .eq. (- one))THEN
2287	DO i = 1, m
2288	x(i) = - x(i)
2289	ENDDO
2290	ELSEIF (alpha .eq. zero)THEN
2291	DO i = 1, m
2292	x(i) = zero
2293	ENDDO
2294	ELSE
2295	DO i = 1, m
2296	x(i) = alpha* x(i)
2297	ENDDO
2298	ENDIF
2299	IF ( n .lt. 5)RETURN
2300	ENDIF
2301	mp1 = m + 1
2302	IF (alpha .eq. (- one))THEN
2303	DO i = mp1, n, 5
2304	x(i) = - x(i)
2305	x(i + 1) = - x(i + 1)
2306	x(i + 2) = - x(i + 2)
2307	x(i + 3) = - x(i + 3)
2308	x(i + 4) = - x(i + 4)
2309	ENDDO
2310	ELSEIF (alpha .eq. zero)THEN
2311	DO i = mp1, n, 5
2312	x(i) = zero
2313	x(i + 1) = zero
2314	x(i + 2) = zero
2315	x(i + 3) = zero
2316	x(i + 4) = zero
2317	ENDDO
2318	ELSE
2319	DO i = mp1, n, 5
2320	x(i) = alpha* x(i)
2321	x(i + 1) = alpha* x(i + 1)
2322	x(i + 2) = alpha* x(i + 2)
2323	x(i + 3) = alpha* x(i + 3)
2324	x(i + 4) = alpha* x(i + 4)
2325	ENDDO
2326	ENDIF
2327
2328	END SUBROUTINE wscal
2329
2330	SUBROUTINE waxpy(n, alpha, x, incx, y, incy)
2331	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2332	! constant times a vector plus a vector: y <- y + Alpha*x
2333	! only for incX=incY=1
2334	! after BLAS
2335	! replace this by the function from the optimized BLAS implementation:
2336	! CALL SAXPY(N,Alpha,X,1,Y,1) or CALL DAXPY(N,Alpha,X,1,Y,1)
2337	!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2338
2339	INTEGER :: i, incx, incy, m, mp1, n
2340	REAL(kind=dp):: x(n), y(n), alpha
2341	REAL(kind=dp), PARAMETER :: zero = 0.0_dp
2342
2343	IF (alpha .eq. zero)RETURN
2344	IF (n .le. 0)RETURN
2345
2346	m = mod(n, 4)
2347	IF ( m .ne. 0)THEN
2348	DO i = 1, m
2349	y(i) = y(i) + alpha* x(i)
2350	ENDDO
2351	IF ( n .lt. 4)RETURN
2352	ENDIF
2353	mp1 = m + 1
2354	DO i = mp1, n, 4
2355	y(i) = y(i) + alpha* x(i)
2356	y(i + 1) = y(i + 1) + alpha* x(i + 1)
2357	y(i + 2) = y(i + 2) + alpha* x(i + 2)
2358	y(i + 3) = y(i + 3) + alpha* x(i + 3)
2359	ENDDO
2360
2361	END SUBROUTINE waxpy
2362
2363	SUBROUTINE rosenbrock(n, y, tstart, tend, &
2364	abstol, reltol, &
2365	rcntrl, icntrl, rstatus, istatus, ierr)
2366	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2367	!
2368	! Solves the system y'=F(t,y) using a Rosenbrock method defined by:
2369	!
2370	! G = 1/(H*gamma(1)) - Jac(t0,Y0)
2371	! T_i = t0 + Alpha(i)*H
2372	! Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
2373	! G K_i = Fun( T_i,Y_i)+ \sum_{j=1}^S C(i,j)/H K_j +
2374	! gamma(i)*dF/dT(t0,Y0)
2375	! Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
2376	!
2377	! For details on Rosenbrock methods and their implementation consult:
2378	! E. Hairer and G. Wanner
2379	! "Solving ODEs II. Stiff and differential-algebraic problems".
2380	! Springer series in computational mathematics,Springer-Verlag,1996.
2381	! The codes contained in the book inspired this implementation.
2382	!
2383	! (C) Adrian Sandu,August 2004
2384	! Virginia Polytechnic Institute and State University
2385	! Contact: sandu@cs.vt.edu
2386	! Revised by Philipp Miehe and Adrian Sandu,May 2006
2387	! This implementation is part of KPP - the Kinetic PreProcessor
2388	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2389	!
2390	!~~~> input arguments:
2391	!
2392	!- y(n) = vector of initial conditions (at t=tstart)
2393	!- [tstart, tend] = time range of integration
2394	! (if Tstart>Tend the integration is performed backwards in time)
2395	!- reltol, abstol = user precribed accuracy
2396	!- SUBROUTINE fun( t, y, ydot) = ode FUNCTION,
2397	! returns Ydot = Y' = F(T,Y)
2398	!- SUBROUTINE jac( t, y, jcb) = jacobian of the ode FUNCTION,
2399	! returns Jcb = dFun/dY
2400	!- icntrl(1:20) = INTEGER inputs PARAMETERs
2401	!- rcntrl(1:20) = REAL inputs PARAMETERs
2402	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2403	!
2404	!~~~> output arguments:
2405	!
2406	!- y(n) - > vector of final states (at t- >tend)
2407	!- istatus(1:20) - > INTEGER output PARAMETERs
2408	!- rstatus(1:20) - > REAL output PARAMETERs
2409	!- ierr - > job status upon RETURN
2410	! success (positive value) or
2411	! failure (negative value)
2412	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2413	!
2414	!~~~> input PARAMETERs:
2415	!
2416	! Note: For input parameters equal to zero the default values of the
2417	! corresponding variables are used.
2418	!
2419	! ICNTRL(1) = 1: F = F(y) Independent of T (AUTONOMOUS)
2420	! = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
2421	!
2422	! ICNTRL(2) = 0: AbsTol,RelTol are N-dimensional vectors
2423	! = 1: AbsTol,RelTol are scalars
2424	!
2425	! ICNTRL(3) -> selection of a particular Rosenbrock method
2426	! = 0 : Rodas3 (default)
2427	! = 1 : Ros2
2428	! = 2 : Ros3
2429	! = 3 : Ros4
2430	! = 4 : Rodas3
2431	! = 5 : Rodas4
2432	!
2433	! ICNTRL(4) -> maximum number of integration steps
2434	! For ICNTRL(4) =0) the default value of 100000 is used
2435	!
2436	! RCNTRL(1) -> Hmin,lower bound for the integration step size
2437	! It is strongly recommended to keep Hmin = ZERO
2438	! RCNTRL(2) -> Hmax,upper bound for the integration step size
2439	! RCNTRL(3) -> Hstart,starting value for the integration step size
2440	!
2441	! RCNTRL(4) -> FacMin,lower bound on step decrease factor (default=0.2)
2442	! RCNTRL(5) -> FacMax,upper bound on step increase factor (default=6)
2443	! RCNTRL(6) -> FacRej,step decrease factor after multiple rejections
2444	! (default=0.1)
2445	! RCNTRL(7) -> FacSafe,by which the new step is slightly smaller
2446	! than the predicted value (default=0.9)
2447	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2448	!
2449	!
2450	! OUTPUT ARGUMENTS:
2451	! -----------------
2452	!
2453	! T -> T value for which the solution has been computed
2454	! (after successful return T=Tend).
2455	!
2456	! Y(N) -> Numerical solution at T
2457	!
2458	! IDID -> Reports on successfulness upon return:
2459	! = 1 for success
2460	! < 0 for error (value equals error code)
2461	!
2462	! ISTATUS(1) -> No. of function calls
2463	! ISTATUS(2) -> No. of jacobian calls
2464	! ISTATUS(3) -> No. of steps
2465	! ISTATUS(4) -> No. of accepted steps
2466	! ISTATUS(5) -> No. of rejected steps (except at very beginning)
2467	! ISTATUS(6) -> No. of LU decompositions
2468	! ISTATUS(7) -> No. of forward/backward substitutions
2469	! ISTATUS(8) -> No. of singular matrix decompositions
2470	!
2471	! RSTATUS(1) -> Texit,the time corresponding to the
2472	! computed Y upon return
2473	! RSTATUS(2) -> Hexit,last accepted step before exit
2474	! RSTATUS(3) -> Hnew,last predicted step (not yet taken)
2475	! For multiple restarts,use Hnew as Hstart
2476	! in the subsequent run
2477	!
2478	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2479
2480
2481	!~~~> arguments
2482	INTEGER, INTENT(IN) :: n
2483	REAL(kind=dp), INTENT(INOUT):: y(n)
2484	REAL(kind=dp), INTENT(IN) :: tstart, tend
2485	REAL(kind=dp), INTENT(IN) :: abstol(n), reltol(n)
2486	INTEGER, INTENT(IN) :: icntrl(20)
2487	REAL(kind=dp), INTENT(IN) :: rcntrl(20)
2488	INTEGER, INTENT(INOUT):: istatus(20)
2489	REAL(kind=dp), INTENT(INOUT):: rstatus(20)
2490	INTEGER, INTENT(OUT) :: ierr
2491	!~~~> PARAMETERs of the rosenbrock method, up to 6 stages
2492	INTEGER :: ros_s, rosmethod
2493	INTEGER, PARAMETER :: rs2=1, rs3=2, rs4=3, rd3=4, rd4=5, rg3=6
2494	REAL(kind=dp):: ros_a(15), ros_c(15), ros_m(6), ros_e(6), &
2495	ros_alpha(6), ros_gamma(6), ros_elo
2496	LOGICAL :: ros_newf(6)
2497	CHARACTER(len=12):: ros_name
2498	!~~~> local variables
2499	REAL(kind=dp):: roundoff, facmin, facmax, facrej, facsafe
2500	REAL(kind=dp):: hmin, hmax, hstart
2501	REAL(kind=dp):: texit
2502	INTEGER :: i, uplimtol, max_no_steps
2503	LOGICAL :: autonomous, vectortol
2504	!~~~> PARAMETERs
2505	REAL(kind=dp), PARAMETER :: zero = 0.0_dp, one = 1.0_dp
2506	REAL(kind=dp), PARAMETER :: deltamin = 1.0e-5_dp
2507
2508	!~~~> initialize statistics
2509	istatus(1:8) = 0
2510	rstatus(1:3) = zero
2511
2512	!~~~> autonomous or time dependent ode. default is time dependent.
2513	autonomous = .not.(icntrl(1) == 0)
2514
2515	!~~~> for scalar tolerances (icntrl(2).ne.0) the code uses abstol(1)and reltol(1)
2516	! For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:N) and RelTol(1:N)
2517	IF (icntrl(2) == 0)THEN
2518	vectortol = .TRUE.
2519	uplimtol = n
2520	ELSE
2521	vectortol = .FALSE.
2522	uplimtol = 1
2523	ENDIF
2524
2525	!~~~> initialize the particular rosenbrock method selected
2526	select CASE (icntrl(3))
2527	CASE (1)
2528	CALL ros2
2529	CASE (2)
2530	CALL ros3
2531	CASE (3)
2532	CALL ros4
2533	CASE (0, 4)
2534	CALL rodas3
2535	CASE (5)
2536	CALL rodas4
2537	CASE (6)
2538	CALL rang3
2539	CASE default
2540	PRINT *,'Unknown Rosenbrock method: ICNTRL(3) =',ICNTRL(3)
2541	CALL ros_errormsg(- 2, tstart, zero, ierr)
2542	RETURN
2543	END select
2544
2545	!~~~> the maximum number of steps admitted
2546	IF (icntrl(4) == 0)THEN
2547	max_no_steps = 200000
2548	ELSEIF (icntrl(4)> 0)THEN
2549	max_no_steps=icntrl(4)
2550	ELSE
2551	PRINT *,'User-selected max no. of steps: ICNTRL(4) =',ICNTRL(4)
2552	CALL ros_errormsg(- 1, tstart, zero, ierr)
2553	RETURN
2554	ENDIF
2555
2556	!~~~> unit roundoff (1+ roundoff>1)
2557	roundoff = epsilon(one)
2558
2559	!~~~> lower bound on the step size: (positive value)
2560	IF (rcntrl(1) == zero)THEN
2561	hmin = zero
2562	ELSEIF (rcntrl(1)> zero)THEN
2563	hmin = rcntrl(1)
2564	ELSE
2565	PRINT *,'User-selected Hmin: RCNTRL(1) =',RCNTRL(1)
2566	CALL ros_errormsg(- 3, tstart, zero, ierr)
2567	RETURN
2568	ENDIF
2569	!~~~> upper bound on the step size: (positive value)
2570	IF (rcntrl(2) == zero)THEN
2571	hmax = abs(tend-tstart)
2572	ELSEIF (rcntrl(2)> zero)THEN
2573	hmax = min(abs(rcntrl(2)), abs(tend-tstart))
2574	ELSE
2575	PRINT *,'User-selected Hmax: RCNTRL(2) =',RCNTRL(2)
2576	CALL ros_errormsg(- 3, tstart, zero, ierr)
2577	RETURN
2578	ENDIF
2579	!~~~> starting step size: (positive value)
2580	IF (rcntrl(3) == zero)THEN
2581	hstart = max(hmin, deltamin)
2582	ELSEIF (rcntrl(3)> zero)THEN
2583	hstart = min(abs(rcntrl(3)), abs(tend-tstart))
2584	ELSE
2585	PRINT *,'User-selected Hstart: RCNTRL(3) =',RCNTRL(3)
2586	CALL ros_errormsg(- 3, tstart, zero, ierr)
2587	RETURN
2588	ENDIF
2589	!~~~> step size can be changed s.t. facmin < hnew/hold < facmax
2590	IF (rcntrl(4) == zero)THEN
2591	facmin = 0.2_dp
2592	ELSEIF (rcntrl(4)> zero)THEN
2593	facmin = rcntrl(4)
2594	ELSE
2595	PRINT *,'User-selected FacMin: RCNTRL(4) =',RCNTRL(4)
2596	CALL ros_errormsg(- 4, tstart, zero, ierr)
2597	RETURN
2598	ENDIF
2599	IF (rcntrl(5) == zero)THEN
2600	facmax = 6.0_dp
2601	ELSEIF (rcntrl(5)> zero)THEN
2602	facmax = rcntrl(5)
2603	ELSE
2604	PRINT *,'User-selected FacMax: RCNTRL(5) =',RCNTRL(5)
2605	CALL ros_errormsg(- 4, tstart, zero, ierr)
2606	RETURN
2607	ENDIF
2608	!~~~> facrej: factor to decrease step after 2 succesive rejections
2609	IF (rcntrl(6) == zero)THEN
2610	facrej = 0.1_dp
2611	ELSEIF (rcntrl(6)> zero)THEN
2612	facrej = rcntrl(6)
2613	ELSE
2614	PRINT *,'User-selected FacRej: RCNTRL(6) =',RCNTRL(6)
2615	CALL ros_errormsg(- 4, tstart, zero, ierr)
2616	RETURN
2617	ENDIF
2618	!~~~> facsafe: safety factor in the computation of new step size
2619	IF (rcntrl(7) == zero)THEN
2620	facsafe = 0.9_dp
2621	ELSEIF (rcntrl(7)> zero)THEN
2622	facsafe = rcntrl(7)
2623	ELSE
2624	PRINT *,'User-selected FacSafe: RCNTRL(7) =',RCNTRL(7)
2625	CALL ros_errormsg(- 4, tstart, zero, ierr)
2626	RETURN
2627	ENDIF
2628	!~~~> check IF tolerances are reasonable
2629	DO i=1, uplimtol
2630	IF ((abstol(i)<= zero).or. (reltol(i)<= 10.0_dp* roundoff)&
2631	.or. (reltol(i)>= 1.0_dp))THEN
2632	PRINT *,' AbsTol(',i,') = ',AbsTol(i)
2633	PRINT *,' RelTol(',i,') = ',RelTol(i)
2634	CALL ros_errormsg(- 5, tstart, zero, ierr)
2635	RETURN
2636	ENDIF
2637	ENDDO
2638
2639
2640	!~~~> CALL rosenbrock method
2641	CALL ros_integrator(y, tstart, tend, texit, &
2642	abstol, reltol, &
2643	! Integration parameters
2644	autonomous, vectortol, max_no_steps, &
2645	roundoff, hmin, hmax, hstart, &
2646	facmin, facmax, facrej, facsafe, &
2647	! Error indicator
2648	ierr)
2649
2650	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2651	CONTAINS ! SUBROUTINEs internal to rosenbrock
2652	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2653
2654	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2655	SUBROUTINE ros_errormsg(code, t, h, ierr)
2656	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2657	! Handles all error messages
2658	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2659
2660	REAL(kind=dp), INTENT(IN):: t, h
2661	INTEGER, INTENT(IN) :: code
2662	INTEGER, INTENT(OUT):: ierr
2663
2664	ierr = code
2665	print * , &
2666	'Forced exit from Rosenbrock due to the following error:'
2667
2668	select CASE (code)
2669	CASE (- 1)
2670	PRINT *,'--> Improper value for maximal no of steps'
2671	CASE (- 2)
2672	PRINT *,'--> Selected Rosenbrock method not implemented'
2673	CASE (- 3)
2674	PRINT *,'--> Hmin/Hmax/Hstart must be positive'
2675	CASE (- 4)
2676	PRINT *,'--> FacMin/FacMax/FacRej must be positive'
2677	CASE (- 5)
2678	PRINT *,'--> Improper tolerance values'
2679	CASE (- 6)
2680	PRINT *,'--> No of steps exceeds maximum bound'
2681	CASE (- 7)
2682	PRINT ,'--> Step size too small: T + 10H = T',&
2683	' or H < Roundoff'
2684	CASE (- 8)
2685	PRINT *,'--> Matrix is repeatedly singular'
2686	CASE default
2687	PRINT *,'Unknown Error code: ',Code
2688	END select
2689
2690	print * , "t=", t, "and h=", h
2691
2692	END SUBROUTINE ros_errormsg
2693
2694	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2695	SUBROUTINE ros_integrator (y, tstart, tend, t, &
2696	abstol, reltol, &
2697	!~~~> integration PARAMETERs
2698	autonomous, vectortol, max_no_steps, &
2699	roundoff, hmin, hmax, hstart, &
2700	facmin, facmax, facrej, facsafe, &
2701	!~~~> error indicator
2702	ierr)
2703	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2704	! Template for the implementation of a generic Rosenbrock method
2705	! defined by ros_S (no of stages)
2706	! and its coefficients ros_{A,C,M,E,Alpha,Gamma}
2707	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2708
2709
2710	!~~~> input: the initial condition at tstart; output: the solution at t
2711	REAL(kind=dp), INTENT(INOUT):: y(n)
2712	!~~~> input: integration interval
2713	REAL(kind=dp), INTENT(IN):: tstart, tend
2714	!~~~> output: time at which the solution is RETURNed (t=tendIF success)
2715	REAL(kind=dp), INTENT(OUT):: t
2716	!~~~> input: tolerances
2717	REAL(kind=dp), INTENT(IN):: abstol(n), reltol(n)
2718	!~~~> input: integration PARAMETERs
2719	LOGICAL, INTENT(IN):: autonomous, vectortol
2720	REAL(kind=dp), INTENT(IN):: hstart, hmin, hmax
2721	INTEGER, INTENT(IN):: max_no_steps
2722	REAL(kind=dp), INTENT(IN):: roundoff, facmin, facmax, facrej, facsafe
2723	!~~~> output: error indicator
2724	INTEGER, INTENT(OUT):: ierr
2725	! ~~~~ Local variables
2726	REAL(kind=dp):: ynew(n), fcn0(n), fcn(n)
2727	REAL(kind=dp):: k(n* ros_s), dfdt(n)
2728	#ifdef full_algebra
2729	REAL(kind=dp):: jac0(n, n), ghimj(n, n)
2730	#else
2731	REAL(kind=dp):: jac0(lu_nonzero), ghimj(lu_nonzero)
2732	#endif
2733	REAL(kind=dp):: h, hnew, hc, hg, fac, tau
2734	REAL(kind=dp):: err, yerr(n)
2735	INTEGER :: pivot(n), direction, ioffset, j, istage
2736	LOGICAL :: rejectlasth, rejectmoreh, singular
2737	!~~~> local PARAMETERs
2738	REAL(kind=dp), PARAMETER :: zero = 0.0_dp, one = 1.0_dp
2739	REAL(kind=dp), PARAMETER :: deltamin = 1.0e-5_dp
2740	!~~~> locally called FUNCTIONs
2741	! REAL(kind=dp) WLAMCH
2742	! EXTERNAL WLAMCH
2743	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2744
2745
2746	!~~~> initial preparations
2747	t = tstart
2748	rstatus(nhexit) = zero
2749	h = min( max(abs(hmin), abs(hstart)), abs(hmax))
2750	IF (abs(h)<= 10.0_dp* roundoff)h = deltamin
2751
2752	IF (tend >= tstart)THEN
2753	direction = + 1
2754	ELSE
2755	direction = - 1
2756	ENDIF
2757	h = direction* h
2758
2759	rejectlasth=.FALSE.
2760	rejectmoreh=.FALSE.
2761
2762	!~~~> time loop begins below
2763
2764	timeloop: DO WHILE((direction > 0).and.((t- tend) + roundoff <= zero)&
2765	.or. (direction < 0).and.((tend-t) + roundoff <= zero))
2766
2767	IF (istatus(nstp)> max_no_steps)THEN ! too many steps
2768	CALL ros_errormsg(- 6, t, h, ierr)
2769	RETURN
2770	ENDIF
2771	IF (((t+ 0.1_dp* h) == t).or.(h <= roundoff))THEN ! step size too small
2772	CALL ros_errormsg(- 7, t, h, ierr)
2773	RETURN
2774	ENDIF
2775
2776	!~~~> limit h IF necessary to avoid going beyond tend
2777	h = min(h, abs(tend-t))
2778
2779	!~~~> compute the FUNCTION at current time
2780	CALL funtemplate(t, y, fcn0)
2781	istatus(nfun) = istatus(nfun) + 1
2782
2783	!~~~> compute the FUNCTION derivative with respect to t
2784	IF (.not.autonomous)THEN
2785	CALL ros_funtimederivative(t, roundoff, y, &
2786	fcn0, dfdt)
2787	ENDIF
2788
2789	!~~~> compute the jacobian at current time
2790	CALL jactemplate(t, y, jac0)
2791	istatus(njac) = istatus(njac) + 1
2792
2793	!~~~> repeat step calculation until current step accepted
2794	untilaccepted: do
2795
2796	CALL ros_preparematrix(h, direction, ros_gamma(1), &
2797	jac0, ghimj, pivot, singular)
2798	IF (singular)THEN ! more than 5 consecutive failed decompositions
2799	CALL ros_errormsg(- 8, t, h, ierr)
2800	RETURN
2801	ENDIF
2802
2803	!~~~> compute the stages
2804	stage: DO istage = 1, ros_s
2805
2806	! current istage offset. current istage vector is k(ioffset+ 1:ioffset+ n)
2807	ioffset = n* (istage-1)
2808
2809	! for the 1st istage the FUNCTION has been computed previously
2810	IF (istage == 1)THEN
2811	!slim: CALL wcopy(n, fcn0, 1, fcn, 1)
2812	fcn(1:n) = fcn0(1:n)
2813	! istage>1 and a new FUNCTION evaluation is needed at the current istage
2814	ELSEIF(ros_newf(istage))THEN
2815	!slim: CALL wcopy(n, y, 1, ynew, 1)
2816	ynew(1:n) = y(1:n)
2817	DO j = 1, istage-1
2818	CALL waxpy(n, ros_a((istage-1) * (istage-2) /2+ j), &
2819	k(n* (j- 1) + 1), 1, ynew, 1)
2820	ENDDO
2821	tau = t + ros_alpha(istage) * direction* h
2822	CALL funtemplate(tau, ynew, fcn)
2823	istatus(nfun) = istatus(nfun) + 1
2824	ENDIF ! IF istage == 1 ELSEIF ros_newf(istage)
2825	!slim: CALL wcopy(n, fcn, 1, k(ioffset+ 1), 1)
2826	k(ioffset+ 1:ioffset+ n) = fcn(1:n)
2827	DO j = 1, istage-1
2828	hc = ros_c((istage-1) * (istage-2) /2+ j) /(direction* h)
2829	CALL waxpy(n, hc, k(n* (j- 1) + 1), 1, k(ioffset+ 1), 1)
2830	ENDDO
2831	IF ((.not. autonomous).and.(ros_gamma(istage).ne.zero))THEN
2832	hg = direction* h* ros_gamma(istage)
2833	CALL waxpy(n, hg, dfdt, 1, k(ioffset+ 1), 1)
2834	ENDIF
2835	CALL ros_solve(ghimj, pivot, k(ioffset+ 1))
2836
2837	END DO stage
2838
2839
2840	!~~~> compute the new solution
2841	!slim: CALL wcopy(n, y, 1, ynew, 1)
2842	ynew(1:n) = y(1:n)
2843	DO j=1, ros_s
2844	CALL waxpy(n, ros_m(j), k(n* (j- 1) + 1), 1, ynew, 1)
2845	ENDDO
2846
2847	!~~~> compute the error estimation
2848	!slim: CALL wscal(n, zero, yerr, 1)
2849	yerr(1:n) = zero
2850	DO j=1, ros_s
2851	CALL waxpy(n, ros_e(j), k(n* (j- 1) + 1), 1, yerr, 1)
2852	ENDDO
2853	err = ros_errornorm(y, ynew, yerr, abstol, reltol, vectortol)
2854
2855	!~~~> new step size is bounded by facmin <= hnew/h <= facmax
2856	fac = min(facmax, max(facmin, facsafe/err** (one/ros_elo)))
2857	hnew = h* fac
2858
2859	!~~~> check the error magnitude and adjust step size
2860	istatus(nstp) = istatus(nstp) + 1
2861	IF ((err <= one).or.(h <= hmin))THEN !~~~> accept step
2862	istatus(nacc) = istatus(nacc) + 1
2863	!slim: CALL wcopy(n, ynew, 1, y, 1)
2864	y(1:n) = ynew(1:n)
2865	t = t + direction* h
2866	hnew = max(hmin, min(hnew, hmax))
2867	IF (rejectlasth)THEN ! no step size increase after a rejected step
2868	hnew = min(hnew, h)
2869	ENDIF
2870	rstatus(nhexit) = h
2871	rstatus(nhnew) = hnew
2872	rstatus(ntexit) = t
2873	rejectlasth = .FALSE.
2874	rejectmoreh = .FALSE.
2875	h = hnew
2876	exit untilaccepted ! exit the loop: WHILE step not accepted
2877	ELSE !~~~> reject step
2878	IF (rejectmoreh)THEN
2879	hnew = h* facrej
2880	ENDIF
2881	rejectmoreh = rejectlasth
2882	rejectlasth = .TRUE.
2883	h = hnew
2884	IF (istatus(nacc)>= 1) istatus(nrej) = istatus(nrej) + 1
2885	ENDIF ! err <= 1
2886
2887	END DO untilaccepted
2888
2889	END DO timeloop
2890
2891	!~~~> succesful exit
2892	ierr = 1 !~~~> the integration was successful
2893
2894	END SUBROUTINE ros_integrator
2895
2896
2897	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2898	REAL(kind=dp)FUNCTION ros_errornorm(y, ynew, yerr, &
2899	abstol, reltol, vectortol)
2900	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2901	!~~~> computes the "scaled norm" of the error vector yerr
2902	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2903
2904	! Input arguments
2905	REAL(kind=dp), INTENT(IN):: y(n), ynew(n), &
2906	yerr(n), abstol(n), reltol(n)
2907	LOGICAL, INTENT(IN):: vectortol
2908	! Local variables
2909	REAL(kind=dp):: err, scale, ymax
2910	INTEGER :: i
2911	REAL(kind=dp), PARAMETER :: zero = 0.0_dp
2912
2913	err = zero
2914	DO i=1, n
2915	ymax = max(abs(y(i)), abs(ynew(i)))
2916	IF (vectortol)THEN
2917	scale = abstol(i) + reltol(i) * ymax
2918	ELSE
2919	scale = abstol(1) + reltol(1) * ymax
2920	ENDIF
2921	err = err+ (yerr(i) /scale) ** 2
2922	ENDDO
2923	err = sqrt(err/n)
2924
2925	ros_errornorm = max(err, 1.0d-10)
2926
2927	END FUNCTION ros_errornorm
2928
2929
2930	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2931	SUBROUTINE ros_funtimederivative(t, roundoff, y, &
2932	fcn0, dfdt)
2933	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2934	!~~~> the time partial derivative of the FUNCTION by finite differences
2935	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2936
2937	!~~~> input arguments
2938	REAL(kind=dp), INTENT(IN):: t, roundoff, y(n), fcn0(n)
2939	!~~~> output arguments
2940	REAL(kind=dp), INTENT(OUT):: dfdt(n)
2941	!~~~> local variables
2942	REAL(kind=dp):: delta
2943	REAL(kind=dp), PARAMETER :: one = 1.0_dp, deltamin = 1.0e-6_dp
2944
2945	delta = sqrt(roundoff) * max(deltamin, abs(t))
2946	CALL funtemplate(t+ delta, y, dfdt)
2947	istatus(nfun) = istatus(nfun) + 1
2948	CALL waxpy(n, (- one), fcn0, 1, dfdt, 1)
2949	CALL wscal(n, (one/delta), dfdt, 1)
2950
2951	END SUBROUTINE ros_funtimederivative
2952
2953
2954	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2955	SUBROUTINE ros_preparematrix(h, direction, gam, &
2956	jac0, ghimj, pivot, singular)
2957	! --- --- --- --- --- --- --- --- --- --- --- --- ---
2958	! Prepares the LHS matrix for stage calculations
2959	! 1. Construct Ghimj = 1/(H*ham) - Jac0
2960	! "(Gamma H) Inverse Minus Jacobian"
2961	! 2. Repeat LU decomposition of Ghimj until successful.
2962	! -half the step size if LU decomposition fails and retry
2963	! -exit after 5 consecutive fails
2964	! --- --- --- --- --- --- --- --- --- --- --- --- ---
2965
2966	!~~~> input arguments
2967	#ifdef full_algebra
2968	REAL(kind=dp), INTENT(IN):: jac0(n, n)
2969	#else
2970	REAL(kind=dp), INTENT(IN):: jac0(lu_nonzero)
2971	#endif
2972	REAL(kind=dp), INTENT(IN):: gam
2973	INTEGER, INTENT(IN):: direction
2974	!~~~> output arguments
2975	#ifdef full_algebra
2976	REAL(kind=dp), INTENT(OUT):: ghimj(n, n)
2977	#else
2978	REAL(kind=dp), INTENT(OUT):: ghimj(lu_nonzero)
2979	#endif
2980	LOGICAL, INTENT(OUT):: singular
2981	INTEGER, INTENT(OUT):: pivot(n)
2982	!~~~> inout arguments
2983	REAL(kind=dp), INTENT(INOUT):: h ! step size is decreased when lu fails
2984	!~~~> local variables
2985	INTEGER :: i, ising, nconsecutive
2986	REAL(kind=dp):: ghinv
2987	REAL(kind=dp), PARAMETER :: one = 1.0_dp, half = 0.5_dp
2988
2989	nconsecutive = 0
2990	singular = .TRUE.
2991
2992	DO WHILE (singular)
2993
2994	!~~~> construct ghimj = 1/(h* gam) - jac0
2995	#ifdef full_algebra
2996	!slim: CALL wcopy(n* n, jac0, 1, ghimj, 1)
2997	!slim: CALL wscal(n* n, (- one), ghimj, 1)
2998	ghimj = - jac0
2999	ghinv = one/(direction* h* gam)
3000	DO i=1, n
3001	ghimj(i, i) = ghimj(i, i) + ghinv
3002	ENDDO
3003	#else
3004	!slim: CALL wcopy(lu_nonzero, jac0, 1, ghimj, 1)
3005	!slim: CALL wscal(lu_nonzero, (- one), ghimj, 1)
3006	ghimj(1:lu_nonzero) = - jac0(1:lu_nonzero)
3007	ghinv = one/(direction* h* gam)
3008	DO i=1, n
3009	ghimj(lu_diag(i)) = ghimj(lu_diag(i)) + ghinv
3010	ENDDO
3011	#endif
3012	!~~~> compute lu decomposition
3013	CALL ros_decomp( ghimj, pivot, ising)
3014	IF (ising == 0)THEN
3015	!~~~> IF successful done
3016	singular = .FALSE.
3017	ELSE ! ising .ne. 0
3018	!~~~> IF unsuccessful half the step size; IF 5 consecutive fails THEN RETURN
3019	istatus(nsng) = istatus(nsng) + 1
3020	nconsecutive = nconsecutive+1
3021	singular = .TRUE.
3022	PRINT*,'Warning: LU Decomposition returned ISING = ',ISING
3023	IF (nconsecutive <= 5)THEN ! less than 5 consecutive failed decompositions
3024	h = h* half
3025	ELSE ! more than 5 consecutive failed decompositions
3026	RETURN
3027	ENDIF ! nconsecutive
3028	ENDIF ! ising
3029
3030	END DO ! WHILE singular
3031
3032	END SUBROUTINE ros_preparematrix
3033
3034
3035	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3036	SUBROUTINE ros_decomp( a, pivot, ising)
3037	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3038	! Template for the LU decomposition
3039	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3040	!~~~> inout variables
3041	#ifdef full_algebra
3042	REAL(kind=dp), INTENT(INOUT):: a(n, n)
3043	#else
3044	REAL(kind=dp), INTENT(INOUT):: a(lu_nonzero)
3045	#endif
3046	!~~~> output variables
3047	INTEGER, INTENT(OUT):: pivot(n), ising
3048
3049	#ifdef full_algebra
3050	CALL dgetrf( n, n, a, n, pivot, ising)
3051	#else
3052	CALL kppdecomp(a, ising)
3053	pivot(1) = 1
3054	#endif
3055	istatus(ndec) = istatus(ndec) + 1
3056
3057	END SUBROUTINE ros_decomp
3058
3059
3060	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3061	SUBROUTINE ros_solve( a, pivot, b)
3062	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3063	! Template for the forward/backward substitution (using pre-computed LU decomposition)
3064	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3065	!~~~> input variables
3066	#ifdef full_algebra
3067	REAL(kind=dp), INTENT(IN):: a(n, n)
3068	INTEGER :: ising
3069	#else
3070	REAL(kind=dp), INTENT(IN):: a(lu_nonzero)
3071	#endif
3072	INTEGER, INTENT(IN):: pivot(n)
3073	!~~~> inout variables
3074	REAL(kind=dp), INTENT(INOUT):: b(n)
3075
3076	#ifdef full_algebra
3077	CALL DGETRS( 'N',N ,1,A,N,Pivot,b,N,ISING)
3078	IF (info < 0)THEN
3079	print* , "error in dgetrs. ising=", ising
3080	ENDIF
3081	#else
3082	CALL kppsolve( a, b)
3083	#endif
3084
3085	istatus(nsol) = istatus(nsol) + 1
3086
3087	END SUBROUTINE ros_solve
3088
3089
3090
3091	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3092	SUBROUTINE ros2
3093	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3094	! --- AN L-STABLE METHOD,2 stages,order 2
3095	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3096
3097	double precision g
3098
3099	g = 1.0_dp + 1.0_dp/sqrt(2.0_dp)
3100	rosmethod = rs2
3101	!~~~> name of the method
3102	ros_Name = 'ROS-2'
3103	!~~~> number of stages
3104	ros_s = 2
3105
3106	!~~~> the coefficient matrices a and c are strictly lower triangular.
3107	! The lower triangular (subdiagonal) elements are stored in row-wise order:
3108	! A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
3109	! The general mapping formula is:
3110	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
3111	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
3112
3113	ros_a(1) = (1.0_dp) /g
3114	ros_c(1) = (- 2.0_dp) /g
3115	!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
3116	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
3117	ros_newf(1) = .TRUE.
3118	ros_newf(2) = .TRUE.
3119	!~~~> m_i = coefficients for new step solution
3120	ros_m(1) = (3.0_dp) /(2.0_dp* g)
3121	ros_m(2) = (1.0_dp) /(2.0_dp* g)
3122	! E_i = Coefficients for error estimator
3123	ros_e(1) = 1.0_dp/(2.0_dp* g)
3124	ros_e(2) = 1.0_dp/(2.0_dp* g)
3125	!~~~> ros_elo = estimator of local order - the minimum between the
3126	! main and the embedded scheme orders plus one
3127	ros_elo = 2.0_dp
3128	!~~~> y_stage_i ~ y( t + h* alpha_i)
3129	ros_alpha(1) = 0.0_dp
3130	ros_alpha(2) = 1.0_dp
3131	!~~~> gamma_i = \sum_j gamma_{i, j}
3132	ros_gamma(1) = g
3133	ros_gamma(2) = -g
3134
3135	END SUBROUTINE ros2
3136
3137
3138	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3139	SUBROUTINE ros3
3140	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3141	! --- AN L-STABLE METHOD,3 stages,order 3,2 function evaluations
3142	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3143
3144	rosmethod = rs3
3145	!~~~> name of the method
3146	ros_Name = 'ROS-3'
3147	!~~~> number of stages
3148	ros_s = 3
3149
3150	!~~~> the coefficient matrices a and c are strictly lower triangular.
3151	! The lower triangular (subdiagonal) elements are stored in row-wise order:
3152	! A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
3153	! The general mapping formula is:
3154	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
3155	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
3156
3157	ros_a(1) = 1.0_dp
3158	ros_a(2) = 1.0_dp
3159	ros_a(3) = 0.0_dp
3160
3161	ros_c(1) = - 0.10156171083877702091975600115545e+01_dp
3162	ros_c(2) = 0.40759956452537699824805835358067e+01_dp
3163	ros_c(3) = 0.92076794298330791242156818474003e+01_dp
3164	!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
3165	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
3166	ros_newf(1) = .TRUE.
3167	ros_newf(2) = .TRUE.
3168	ros_newf(3) = .FALSE.
3169	!~~~> m_i = coefficients for new step solution
3170	ros_m(1) = 0.1e+01_dp
3171	ros_m(2) = 0.61697947043828245592553615689730e+01_dp
3172	ros_m(3) = - 0.42772256543218573326238373806514_dp
3173	! E_i = Coefficients for error estimator
3174	ros_e(1) = 0.5_dp
3175	ros_e(2) = - 0.29079558716805469821718236208017e+01_dp
3176	ros_e(3) = 0.22354069897811569627360909276199_dp
3177	!~~~> ros_elo = estimator of local order - the minimum between the
3178	! main and the embedded scheme orders plus 1
3179	ros_elo = 3.0_dp
3180	!~~~> y_stage_i ~ y( t + h* alpha_i)
3181	ros_alpha(1) = 0.0_dp
3182	ros_alpha(2) = 0.43586652150845899941601945119356_dp
3183	ros_alpha(3) = 0.43586652150845899941601945119356_dp
3184	!~~~> gamma_i = \sum_j gamma_{i, j}
3185	ros_gamma(1) = 0.43586652150845899941601945119356_dp
3186	ros_gamma(2) = 0.24291996454816804366592249683314_dp
3187	ros_gamma(3) = 0.21851380027664058511513169485832e+01_dp
3188
3189	END SUBROUTINE ros3
3190
3191	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3192
3193
3194	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3195	SUBROUTINE ros4
3196	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3197	! L-STABLE ROSENBROCK METHOD OF ORDER 4,WITH 4 STAGES
3198	! L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
3199	!
3200	! E. HAIRER AND G. WANNER,SOLVING ORDINARY DIFFERENTIAL
3201	! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
3202	! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
3203	! SPRINGER-VERLAG (1990)
3204	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3205
3206
3207	rosmethod = rs4
3208	!~~~> name of the method
3209	ros_Name = 'ROS-4'
3210	!~~~> number of stages
3211	ros_s = 4
3212
3213	!~~~> the coefficient matrices a and c are strictly lower triangular.
3214	! The lower triangular (subdiagonal) elements are stored in row-wise order:
3215	! A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
3216	! The general mapping formula is:
3217	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
3218	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
3219
3220	ros_a(1) = 0.2000000000000000e+01_dp
3221	ros_a(2) = 0.1867943637803922e+01_dp
3222	ros_a(3) = 0.2344449711399156_dp
3223	ros_a(4) = ros_a(2)
3224	ros_a(5) = ros_a(3)
3225	ros_a(6) = 0.0_dp
3226
3227	ros_c(1) = -0.7137615036412310e+01_dp
3228	ros_c(2) = 0.2580708087951457e+01_dp
3229	ros_c(3) = 0.6515950076447975_dp
3230	ros_c(4) = -0.2137148994382534e+01_dp
3231	ros_c(5) = -0.3214669691237626_dp
3232	ros_c(6) = -0.6949742501781779_dp
3233	!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
3234	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
3235	ros_newf(1) = .TRUE.
3236	ros_newf(2) = .TRUE.
3237	ros_newf(3) = .TRUE.
3238	ros_newf(4) = .FALSE.
3239	!~~~> m_i = coefficients for new step solution
3240	ros_m(1) = 0.2255570073418735e+01_dp
3241	ros_m(2) = 0.2870493262186792_dp
3242	ros_m(3) = 0.4353179431840180_dp
3243	ros_m(4) = 0.1093502252409163e+01_dp
3244	!~~~> e_i = coefficients for error estimator
3245	ros_e(1) = -0.2815431932141155_dp
3246	ros_e(2) = -0.7276199124938920e-01_dp
3247	ros_e(3) = -0.1082196201495311_dp
3248	ros_e(4) = -0.1093502252409163e+01_dp
3249	!~~~> ros_elo = estimator of local order - the minimum between the
3250	! main and the embedded scheme orders plus 1
3251	ros_elo = 4.0_dp
3252	!~~~> y_stage_i ~ y( t + h* alpha_i)
3253	ros_alpha(1) = 0.0_dp
3254	ros_alpha(2) = 0.1145640000000000e+01_dp
3255	ros_alpha(3) = 0.6552168638155900_dp
3256	ros_alpha(4) = ros_alpha(3)
3257	!~~~> gamma_i = \sum_j gamma_{i, j}
3258	ros_gamma(1) = 0.5728200000000000_dp
3259	ros_gamma(2) = -0.1769193891319233e+01_dp
3260	ros_gamma(3) = 0.7592633437920482_dp
3261	ros_gamma(4) = -0.1049021087100450_dp
3262
3263	END SUBROUTINE ros4
3264
3265	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3266	SUBROUTINE rodas3
3267	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3268	! --- A STIFFLY-STABLE METHOD,4 stages,order 3
3269	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3270
3271
3272	rosmethod = rd3
3273	!~~~> name of the method
3274	ros_Name = 'RODAS-3'
3275	!~~~> number of stages
3276	ros_s = 4
3277
3278	!~~~> the coefficient matrices a and c are strictly lower triangular.
3279	! The lower triangular (subdiagonal) elements are stored in row-wise order:
3280	! A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
3281	! The general mapping formula is:
3282	! A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
3283	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
3284
3285	ros_a(1) = 0.0_dp
3286	ros_a(2) = 2.0_dp
3287	ros_a(3) = 0.0_dp
3288	ros_a(4) = 2.0_dp
3289	ros_a(5) = 0.0_dp
3290	ros_a(6) = 1.0_dp
3291
3292	ros_c(1) = 4.0_dp
3293	ros_c(2) = 1.0_dp
3294	ros_c(3) = -1.0_dp
3295	ros_c(4) = 1.0_dp
3296	ros_c(5) = -1.0_dp
3297	ros_c(6) = -(8.0_dp/3.0_dp)
3298
3299	!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
3300	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
3301	ros_newf(1) = .TRUE.
3302	ros_newf(2) = .FALSE.
3303	ros_newf(3) = .TRUE.
3304	ros_newf(4) = .TRUE.
3305	!~~~> m_i = coefficients for new step solution
3306	ros_m(1) = 2.0_dp
3307	ros_m(2) = 0.0_dp
3308	ros_m(3) = 1.0_dp
3309	ros_m(4) = 1.0_dp
3310	!~~~> e_i = coefficients for error estimator
3311	ros_e(1) = 0.0_dp
3312	ros_e(2) = 0.0_dp
3313	ros_e(3) = 0.0_dp
3314	ros_e(4) = 1.0_dp
3315	!~~~> ros_elo = estimator of local order - the minimum between the
3316	! main and the embedded scheme orders plus 1
3317	ros_elo = 3.0_dp
3318	!~~~> y_stage_i ~ y( t + h* alpha_i)
3319	ros_alpha(1) = 0.0_dp
3320	ros_alpha(2) = 0.0_dp
3321	ros_alpha(3) = 1.0_dp
3322	ros_alpha(4) = 1.0_dp
3323	!~~~> gamma_i = \sum_j gamma_{i, j}
3324	ros_gamma(1) = 0.5_dp
3325	ros_gamma(2) = 1.5_dp
3326	ros_gamma(3) = 0.0_dp
3327	ros_gamma(4) = 0.0_dp
3328
3329	END SUBROUTINE rodas3
3330
3331	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3332	SUBROUTINE rodas4
3333	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3334	! STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4,WITH 6 STAGES
3335	!
3336	! E. HAIRER AND G. WANNER,SOLVING ORDINARY DIFFERENTIAL
3337	! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
3338	! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
3339	! SPRINGER-VERLAG (1996)
3340	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3341
3342
3343	rosmethod = rd4
3344	!~~~> name of the method
3345	ros_Name = 'RODAS-4'
3346	!~~~> number of stages
3347	ros_s = 6
3348
3349	!~~~> y_stage_i ~ y( t + h* alpha_i)
3350	ros_alpha(1) = 0.000_dp
3351	ros_alpha(2) = 0.386_dp
3352	ros_alpha(3) = 0.210_dp
3353	ros_alpha(4) = 0.630_dp
3354	ros_alpha(5) = 1.000_dp
3355	ros_alpha(6) = 1.000_dp
3356
3357	!~~~> gamma_i = \sum_j gamma_{i, j}
3358	ros_gamma(1) = 0.2500000000000000_dp
3359	ros_gamma(2) = -0.1043000000000000_dp
3360	ros_gamma(3) = 0.1035000000000000_dp
3361	ros_gamma(4) = -0.3620000000000023e-01_dp
3362	ros_gamma(5) = 0.0_dp
3363	ros_gamma(6) = 0.0_dp
3364
3365	!~~~> the coefficient matrices a and c are strictly lower triangular.
3366	! The lower triangular (subdiagonal) elements are stored in row-wise order:
3367	! A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
3368	! The general mapping formula is: A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
3369	! C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
3370
3371	ros_a(1) = 0.1544000000000000e+01_dp
3372	ros_a(2) = 0.9466785280815826_dp
3373	ros_a(3) = 0.2557011698983284_dp
3374	ros_a(4) = 0.3314825187068521e+01_dp
3375	ros_a(5) = 0.2896124015972201e+01_dp
3376	ros_a(6) = 0.9986419139977817_dp
3377	ros_a(7) = 0.1221224509226641e+01_dp
3378	ros_a(8) = 0.6019134481288629e+01_dp
3379	ros_a(9) = 0.1253708332932087e+02_dp
3380	ros_a(10) = -0.6878860361058950_dp
3381	ros_a(11) = ros_a(7)
3382	ros_a(12) = ros_a(8)
3383	ros_a(13) = ros_a(9)
3384	ros_a(14) = ros_a(10)
3385	ros_a(15) = 1.0_dp
3386
3387	ros_c(1) = -0.5668800000000000e+01_dp
3388	ros_c(2) = -0.2430093356833875e+01_dp
3389	ros_c(3) = -0.2063599157091915_dp
3390	ros_c(4) = -0.1073529058151375_dp
3391	ros_c(5) = -0.9594562251023355e+01_dp
3392	ros_c(6) = -0.2047028614809616e+02_dp
3393	ros_c(7) = 0.7496443313967647e+01_dp
3394	ros_c(8) = -0.1024680431464352e+02_dp
3395	ros_c(9) = -0.3399990352819905e+02_dp
3396	ros_c(10) = 0.1170890893206160e+02_dp
3397	ros_c(11) = 0.8083246795921522e+01_dp
3398	ros_c(12) = -0.7981132988064893e+01_dp
3399	ros_c(13) = -0.3152159432874371e+02_dp
3400	ros_c(14) = 0.1631930543123136e+02_dp
3401	ros_c(15) = -0.6058818238834054e+01_dp
3402
3403	!~~~> m_i = coefficients for new step solution
3404	ros_m(1) = ros_a(7)
3405	ros_m(2) = ros_a(8)
3406	ros_m(3) = ros_a(9)
3407	ros_m(4) = ros_a(10)
3408	ros_m(5) = 1.0_dp
3409	ros_m(6) = 1.0_dp
3410
3411	!~~~> e_i = coefficients for error estimator
3412	ros_e(1) = 0.0_dp
3413	ros_e(2) = 0.0_dp
3414	ros_e(3) = 0.0_dp
3415	ros_e(4) = 0.0_dp
3416	ros_e(5) = 0.0_dp
3417	ros_e(6) = 1.0_dp
3418
3419	!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
3420	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
3421	ros_newf(1) = .TRUE.
3422	ros_newf(2) = .TRUE.
3423	ros_newf(3) = .TRUE.
3424	ros_newf(4) = .TRUE.
3425	ros_newf(5) = .TRUE.
3426	ros_newf(6) = .TRUE.
3427
3428	!~~~> ros_elo = estimator of local order - the minimum between the
3429	! main and the embedded scheme orders plus 1
3430	ros_elo = 4.0_dp
3431
3432	END SUBROUTINE rodas4
3433	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3434	SUBROUTINE rang3
3435	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3436	! STIFFLY-STABLE W METHOD OF ORDER 3,WITH 4 STAGES
3437	!
3438	! J. RANG and L. ANGERMANN
3439	! NEW ROSENBROCK W-METHODS OF ORDER 3
3440	! FOR PARTIAL DIFFERENTIAL ALGEBRAIC
3441	! EQUATIONS OF INDEX 1
3442	! BIT Numerical Mathematics (2005) 45: 761-787
3443	! DOI: 10.1007/s10543-005-0035-y
3444	! Table 4.1-4.2
3445	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3446
3447
3448	rosmethod = rg3
3449	!~~~> name of the method
3450	ros_Name = 'RANG-3'
3451	!~~~> number of stages
3452	ros_s = 4
3453
3454	ros_a(1) = 5.09052051067020d+00;
3455	ros_a(2) = 5.09052051067020d+00;
3456	ros_a(3) = 0.0d0;
3457	ros_a(4) = 4.97628111010787d+00;
3458	ros_a(5) = 2.77268164715849d-02;
3459	ros_a(6) = 2.29428036027904d-01;
3460
3461	ros_c(1) = - 1.16790812312283d+01;
3462	ros_c(2) = - 1.64057326467367d+01;
3463	ros_c(3) = - 2.77268164715850d-01;
3464	ros_c(4) = - 8.38103960500476d+00;
3465	ros_c(5) = - 8.48328409199343d-01;
3466	ros_c(6) = 2.87009860433106d-01;
3467
3468	ros_m(1) = 5.22582761233094d+00;
3469	ros_m(2) = - 5.56971148154165d-01;
3470	ros_m(3) = 3.57979469353645d-01;
3471	ros_m(4) = 1.72337398521064d+00;
3472
3473	ros_e(1) = - 5.16845212784040d+00;
3474	ros_e(2) = - 1.26351942603842d+00;
3475	ros_e(3) = - 1.11022302462516d-16;
3476	ros_e(4) = 2.22044604925031d-16;
3477
3478	ros_alpha(1) = 0.0d00;
3479	ros_alpha(2) = 2.21878746765329d+00;
3480	ros_alpha(3) = 2.21878746765329d+00;
3481	ros_alpha(4) = 1.55392337535788d+00;
3482
3483	ros_gamma(1) = 4.35866521508459d-01;
3484	ros_gamma(2) = - 1.78292094614483d+00;
3485	ros_gamma(3) = - 2.46541900496934d+00;
3486	ros_gamma(4) = - 8.05529997906370d-01;
3487
3488
3489	!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
3490	! or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
3491	ros_newf(1) = .TRUE.
3492	ros_newf(2) = .TRUE.
3493	ros_newf(3) = .TRUE.
3494	ros_newf(4) = .TRUE.
3495
3496	!~~~> ros_elo = estimator of local order - the minimum between the
3497	! main and the embedded scheme orders plus 1
3498	ros_elo = 3.0_dp
3499
3500	END SUBROUTINE rang3
3501	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3502
3503	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3504	! End of the set of internal Rosenbrock subroutines
3505	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3506	END SUBROUTINE rosenbrock
3507
3508	SUBROUTINE funtemplate( t, y, ydot)
3509	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3510	! Template for the ODE function call.
3511	! Updates the rate coefficients (and possibly the fixed species) at each call
3512	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3513	!~~~> input variables
3514	REAL(kind=dp):: t, y(nvar)
3515	!~~~> output variables
3516	REAL(kind=dp):: ydot(nvar)
3517	!~~~> local variables
3518	REAL(kind=dp):: told
3519
3520	told = time
3521	time = t
3522	CALL fun( y, fix, rconst, ydot)
3523	time = told
3524
3525	END SUBROUTINE funtemplate
3526
3527	SUBROUTINE jactemplate( t, y, jcb)
3528	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3529	! Template for the ODE Jacobian call.
3530	! Updates the rate coefficients (and possibly the fixed species) at each call
3531	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3532	!~~~> input variables
3533	REAL(kind=dp):: t, y(nvar)
3534	!~~~> output variables
3535	#ifdef full_algebra
3536	REAL(kind=dp):: jv(lu_nonzero), jcb(nvar, nvar)
3537	#else
3538	REAL(kind=dp):: jcb(lu_nonzero)
3539	#endif
3540	!~~~> local variables
3541	REAL(kind=dp):: told
3542	#ifdef full_algebra
3543	INTEGER :: i, j
3544	#endif
3545
3546	told = time
3547	time = t
3548	#ifdef full_algebra
3549	CALL jac_sp(y, fix, rconst, jv)
3550	DO j=1, nvar
3551	DO i=1, nvar
3552	jcb(i, j) = 0.0_dp
3553	ENDDO
3554	ENDDO
3555	DO i=1, lu_nonzero
3556	jcb(lu_irow(i), lu_icol(i)) = jv(i)
3557	ENDDO
3558	#else
3559	CALL jac_sp( y, fix, rconst, jcb)
3560	#endif
3561	time = told
3562
3563	END SUBROUTINE jactemplate
3564
3565	SUBROUTINE kppdecomp( jvs, ier)
3566	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3567	! sparse lu factorization
3568	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3569	! loop expansion generated by kp4
3570
3571	INTEGER :: ier
3572	REAL(kind=dp):: jvs(lu_nonzero), w(nvar), a
3573	INTEGER :: k, kk, j, jj
3574
3575	a = 0.
3576	ier = 0
3577
3578	! i = 1
3579	! i = 2
3580	! i = 3
3581	! i = 4
3582	! i = 5
3583	! i = 6
3584	! i = 7
3585	! i = 8
3586	! i = 9
3587	! i = 10
3588	! i = 11
3589	jvs(38) = (jvs(38)) / jvs(14)
3590	jvs(39) = (jvs(39)) / jvs(19)
3591	jvs(41) = jvs(41) - jvs(15) * jvs(38) - jvs(20) * jvs(39)
3592	! i = 12
3593	jvs(43) = (jvs(43)) / jvs(16)
3594	jvs(48) = jvs(48) - jvs(17) * jvs(43)
3595	jvs(50) = jvs(50) - jvs(18) * jvs(43)
3596	! i = 13
3597	! i = 14
3598	jvs(55) = (jvs(55)) / jvs(14)
3599	jvs(56) = (jvs(56)) / jvs(19)
3600	jvs(57) = (jvs(57)) / jvs(40)
3601	jvs(59) = jvs(59) - jvs(15) * jvs(55) - jvs(20) * jvs(56) - jvs(41) * jvs(57)
3602	jvs(61) = jvs(61) - jvs(42) * jvs(57)
3603	! i = 15
3604	jvs(62) = (jvs(62)) / jvs(19)
3605	jvs(67) = jvs(67) - jvs(20) * jvs(62)
3606	! i = 16
3607	jvs(68) = (jvs(68)) / jvs(63)
3608	jvs(71) = jvs(71) - jvs(64) * jvs(68)
3609	jvs(73) = jvs(73) - jvs(65) * jvs(68)
3610	jvs(76) = jvs(76) - jvs(66) * jvs(68)
3611	jvs(77) = jvs(77) - jvs(67) * jvs(68)
3612	! i = 17
3613	! i = 18
3614	jvs(85) = (jvs(85)) / jvs(14)
3615	jvs(86) = (jvs(86)) / jvs(19)
3616	jvs(87) = (jvs(87)) / jvs(51)
3617	jvs(88) = (jvs(88)) / jvs(58)
3618	jvs(89) = (jvs(89)) / jvs(63)
3619	jvs(90) = (jvs(90)) / jvs(80)
3620	jvs(92) = jvs(92) - jvs(64) * jvs(89)
3621	jvs(93) = jvs(93) - jvs(52) * jvs(87)
3622	jvs(94) = jvs(94) - jvs(65) * jvs(89) - jvs(81) * jvs(90)
3623	jvs(97) = jvs(97) - jvs(66) * jvs(89) - jvs(82) * jvs(90)
3624	jvs(98) = jvs(98) - jvs(53) * jvs(87)
3625	jvs(99) = jvs(99) - jvs(15) * jvs(85) - jvs(20) * jvs(86) - jvs(54) * jvs(87) - jvs(59) * jvs(88)&
3626	- jvs(67) * jvs(89) - jvs(83) * jvs(90)
3627	jvs(101) = jvs(101) - jvs(84) * jvs(90)
3628	jvs(102) = jvs(102) - jvs(60) * jvs(88)
3629	jvs(103) = jvs(103) - jvs(61) * jvs(88)
3630	! i = 19
3631	jvs(105) = (jvs(105)) / jvs(40)
3632	jvs(106) = (jvs(106)) / jvs(58)
3633	jvs(109) = jvs(109) - jvs(41) * jvs(105) - jvs(59) * jvs(106)
3634	jvs(110) = jvs(110) - jvs(60) * jvs(106)
3635	jvs(111) = jvs(111) - jvs(42) * jvs(105) - jvs(61) * jvs(106)
3636	! i = 20
3637	jvs(112) = (jvs(112)) / jvs(19)
3638	jvs(113) = (jvs(113)) / jvs(51)
3639	jvs(114) = jvs(114) - jvs(52) * jvs(113)
3640	jvs(118) = jvs(118) - jvs(53) * jvs(113)
3641	jvs(119) = jvs(119) - jvs(20) * jvs(112) - jvs(54) * jvs(113)
3642	! i = 21
3643	jvs(122) = (jvs(122)) / jvs(80)
3644	jvs(123) = (jvs(123)) / jvs(107)
3645	jvs(125) = jvs(125) - jvs(81) * jvs(122)
3646	jvs(128) = jvs(128) - jvs(82) * jvs(122) - jvs(108) * jvs(123)
3647	jvs(129) = jvs(129) - jvs(83) * jvs(122) - jvs(109) * jvs(123)
3648	jvs(131) = jvs(131) - jvs(84) * jvs(122)
3649	jvs(132) = jvs(132) - jvs(110) * jvs(123)
3650	jvs(133) = jvs(133) - jvs(111) * jvs(123)
3651	! i = 22
3652	! i = 23
3653	jvs(140) = (jvs(140)) / jvs(135)
3654	jvs(142) = jvs(142) - jvs(136) * jvs(140)
3655	jvs(143) = jvs(143) - jvs(137) * jvs(140)
3656	jvs(144) = jvs(144) - jvs(138) * jvs(140)
3657	jvs(145) = jvs(145) - jvs(139) * jvs(140)
3658	! i = 24
3659	jvs(146) = (jvs(146)) / jvs(51)
3660	jvs(147) = (jvs(147)) / jvs(80)
3661	jvs(148) = (jvs(148)) / jvs(107)
3662	a = 0.0; a = a - jvs(52) * jvs(146)
3663	jvs(149) = (jvs(149) + a) / jvs(114)
3664	a = 0.0; a = a - jvs(81) * jvs(147) - jvs(115) * jvs(149)
3665	jvs(150) = (jvs(150) + a) / jvs(135)
3666	a = 0.0; a = a - jvs(116) * jvs(149)
3667	jvs(151) = (jvs(151) + a) / jvs(141)
3668	jvs(153) = jvs(153) - jvs(82) * jvs(147) - jvs(108) * jvs(148) - jvs(117) * jvs(149) - jvs(136) * jvs(150)&
3669	- jvs(142) * jvs(151)
3670	jvs(154) = jvs(154) - jvs(53) * jvs(146) - jvs(118) * jvs(149)
3671	jvs(155) = jvs(155) - jvs(54) * jvs(146) - jvs(83) * jvs(147) - jvs(109) * jvs(148) - jvs(119) * jvs(149)&
3672	- jvs(137) * jvs(150) - jvs(143) * jvs(151)
3673	jvs(156) = jvs(156) - jvs(84) * jvs(147) - jvs(120) * jvs(149) - jvs(138) * jvs(150) - jvs(144) * jvs(151)
3674	jvs(157) = jvs(157) - jvs(110) * jvs(148) - jvs(121) * jvs(149) - jvs(139) * jvs(150) - jvs(145) * jvs(151)
3675	jvs(158) = jvs(158) - jvs(111) * jvs(148)
3676	! i = 25
3677	jvs(159) = (jvs(159)) / jvs(80)
3678	jvs(160) = (jvs(160)) / jvs(107)
3679	a = 0.0; a = a - jvs(81) * jvs(159)
3680	jvs(161) = (jvs(161) + a) / jvs(135)
3681	jvs(162) = (jvs(162)) / jvs(141)
3682	jvs(163) = jvs(163) - jvs(82) * jvs(159) - jvs(108) * jvs(160) - jvs(136) * jvs(161) - jvs(142) * jvs(162)
3683	jvs(165) = jvs(165) - jvs(83) * jvs(159) - jvs(109) * jvs(160) - jvs(137) * jvs(161) - jvs(143) * jvs(162)
3684	jvs(167) = jvs(167) - jvs(84) * jvs(159) - jvs(138) * jvs(161) - jvs(144) * jvs(162)
3685	jvs(168) = jvs(168) - jvs(110) * jvs(160) - jvs(139) * jvs(161) - jvs(145) * jvs(162)
3686	jvs(169) = jvs(169) - jvs(111) * jvs(160)
3687	! i = 26
3688	jvs(170) = (jvs(170)) / jvs(6)
3689	jvs(171) = (jvs(171)) / jvs(9)
3690	jvs(172) = (jvs(172)) / jvs(16)
3691	jvs(173) = (jvs(173)) / jvs(30)
3692	jvs(174) = (jvs(174)) / jvs(34)
3693	jvs(175) = (jvs(175)) / jvs(40)
3694	jvs(176) = (jvs(176)) / jvs(51)
3695	a = 0.0; a = a - jvs(10) * jvs(171)
3696	jvs(177) = (jvs(177) + a) / jvs(58)
3697	jvs(178) = (jvs(178)) / jvs(91)
3698	a = 0.0; a = a - jvs(92) * jvs(178)
3699	jvs(179) = (jvs(179) + a) / jvs(107)
3700	a = 0.0; a = a - jvs(52) * jvs(176) - jvs(93) * jvs(178)
3701	jvs(180) = (jvs(180) + a) / jvs(114)
3702	a = 0.0; a = a - jvs(94) * jvs(178) - jvs(115) * jvs(180)
3703	jvs(181) = (jvs(181) + a) / jvs(135)
3704	a = 0.0; a = a - jvs(95) * jvs(178) - jvs(116) * jvs(180)
3705	jvs(182) = (jvs(182) + a) / jvs(141)
3706	a = 0.0; a = a - jvs(96) * jvs(178)
3707	jvs(183) = (jvs(183) + a) / jvs(152)
3708	a = 0.0; a = a - jvs(97) * jvs(178) - jvs(108) * jvs(179) - jvs(117) * jvs(180) - jvs(136) * jvs(181)&
3709	- jvs(142) * jvs(182) - jvs(153) * jvs(183)
3710	jvs(184) = (jvs(184) + a) / jvs(163)
3711	jvs(185) = jvs(185) - jvs(7) * jvs(170) - jvs(11) * jvs(171) - jvs(17) * jvs(172) - jvs(31) * jvs(173)&
3712	- jvs(35) * jvs(174) - jvs(53) * jvs(176) - jvs(98) * jvs(178) - jvs(118) * jvs(180) - jvs(154) * jvs(183)&
3713	- jvs(164) * jvs(184)
3714	jvs(186) = jvs(186) - jvs(12) * jvs(171) - jvs(32) * jvs(173) - jvs(36) * jvs(174) - jvs(41) * jvs(175)&
3715	- jvs(54) * jvs(176) - jvs(59) * jvs(177) - jvs(99) * jvs(178) - jvs(109) * jvs(179) - jvs(119) * jvs(180)&
3716	- jvs(137) * jvs(181) - jvs(143) * jvs(182) - jvs(155) * jvs(183) - jvs(165) * jvs(184)
3717	jvs(187) = jvs(187) - jvs(37) * jvs(174) - jvs(100) * jvs(178) - jvs(166) * jvs(184)
3718	jvs(188) = jvs(188) - jvs(101) * jvs(178) - jvs(120) * jvs(180) - jvs(138) * jvs(181) - jvs(144) * jvs(182)&
3719	- jvs(156) * jvs(183) - jvs(167) * jvs(184)
3720	jvs(189) = jvs(189) - jvs(13) * jvs(171) - jvs(18) * jvs(172) - jvs(60) * jvs(177) - jvs(102) * jvs(178)&
3721	- jvs(110) * jvs(179) - jvs(121) * jvs(180) - jvs(139) * jvs(181) - jvs(145) * jvs(182) - jvs(157) * jvs(183)&
3722	- jvs(168) * jvs(184)
3723	jvs(190) = jvs(190) - jvs(33) * jvs(173) - jvs(42) * jvs(175) - jvs(61) * jvs(177) - jvs(103) * jvs(178)&
3724	- jvs(111) * jvs(179) - jvs(158) * jvs(183) - jvs(169) * jvs(184)
3725	jvs(191) = jvs(191) - jvs(8) * jvs(170) - jvs(104) * jvs(178)
3726	! i = 27
3727	jvs(192) = (jvs(192)) / jvs(1)
3728	jvs(193) = (jvs(193)) / jvs(3)
3729	jvs(194) = (jvs(194)) / jvs(14)
3730	jvs(195) = (jvs(195)) / jvs(19)
3731	jvs(196) = (jvs(196)) / jvs(30)
3732	jvs(197) = (jvs(197)) / jvs(34)
3733	jvs(198) = (jvs(198)) / jvs(44)
3734	a = 0.0; a = a - jvs(45) * jvs(198)
3735	jvs(199) = (jvs(199) + a) / jvs(58)
3736	jvs(200) = (jvs(200)) / jvs(63)
3737	jvs(201) = (jvs(201)) / jvs(69)
3738	a = 0.0; a = a - jvs(70) * jvs(201)
3739	jvs(202) = (jvs(202) + a) / jvs(80)
3740	a = 0.0; a = a - jvs(64) * jvs(200) - jvs(71) * jvs(201)
3741	jvs(203) = (jvs(203) + a) / jvs(107)
3742	jvs(204) = (jvs(204)) / jvs(114)
3743	a = 0.0; a = a - jvs(46) * jvs(198) - jvs(72) * jvs(201)
3744	jvs(205) = (jvs(205) + a) / jvs(124)
3745	a = 0.0; a = a - jvs(65) * jvs(200) - jvs(73) * jvs(201) - jvs(81) * jvs(202) - jvs(115) * jvs(204)&
3746	- jvs(125) * jvs(205)
3747	jvs(206) = (jvs(206) + a) / jvs(135)
3748	a = 0.0; a = a - jvs(74) * jvs(201) - jvs(116) * jvs(204) - jvs(126) * jvs(205)
3749	jvs(207) = (jvs(207) + a) / jvs(141)
3750	a = 0.0; a = a - jvs(47) * jvs(198) - jvs(75) * jvs(201) - jvs(127) * jvs(205)
3751	jvs(208) = (jvs(208) + a) / jvs(152)
3752	a = 0.0; a = a - jvs(2) * jvs(192) - jvs(66) * jvs(200) - jvs(76) * jvs(201) - jvs(82) * jvs(202)&
3753	- jvs(108) * jvs(203) - jvs(117) * jvs(204) - jvs(128) * jvs(205) - jvs(136) * jvs(206) - jvs(142) * jvs(207)&
3754	- jvs(153) * jvs(208)
3755	jvs(209) = (jvs(209) + a) / jvs(163)
3756	a = 0.0; a = a - jvs(31) * jvs(196) - jvs(35) * jvs(197) - jvs(48) * jvs(198) - jvs(118) * jvs(204)&
3757	- jvs(154) * jvs(208) - jvs(164) * jvs(209)
3758	jvs(210) = (jvs(210) + a) / jvs(185)
3759	jvs(211) = jvs(211) - jvs(4) * jvs(193) - jvs(15) * jvs(194) - jvs(20) * jvs(195) - jvs(32) * jvs(196)&
3760	- jvs(36) * jvs(197) - jvs(49) * jvs(198) - jvs(59) * jvs(199) - jvs(67) * jvs(200) - jvs(77) * jvs(201)&
3761	- jvs(83) * jvs(202) - jvs(109) * jvs(203) - jvs(119) * jvs(204) - jvs(129) * jvs(205) - jvs(137) * jvs(206)&
3762	- jvs(143) * jvs(207) - jvs(155) * jvs(208) - jvs(165) * jvs(209) - jvs(186) * jvs(210)
3763	jvs(212) = jvs(212) - jvs(5) * jvs(193) - jvs(37) * jvs(197) - jvs(130) * jvs(205) - jvs(166) * jvs(209)&
3764	- jvs(187) * jvs(210)
3765	jvs(213) = jvs(213) - jvs(78) * jvs(201) - jvs(84) * jvs(202) - jvs(120) * jvs(204) - jvs(131) * jvs(205)&
3766	- jvs(138) * jvs(206) - jvs(144) * jvs(207) - jvs(156) * jvs(208) - jvs(167) * jvs(209) - jvs(188) * jvs(210)
3767	jvs(214) = jvs(214) - jvs(50) * jvs(198) - jvs(60) * jvs(199) - jvs(79) * jvs(201) - jvs(110) * jvs(203)&
3768	- jvs(121) * jvs(204) - jvs(132) * jvs(205) - jvs(139) * jvs(206) - jvs(145) * jvs(207) - jvs(157) * jvs(208)&
3769	- jvs(168) * jvs(209) - jvs(189) * jvs(210)
3770	jvs(215) = jvs(215) - jvs(33) * jvs(196) - jvs(61) * jvs(199) - jvs(111) * jvs(203) - jvs(133) * jvs(205)&
3771	- jvs(158) * jvs(208) - jvs(169) * jvs(209) - jvs(190) * jvs(210)
3772	jvs(216) = jvs(216) - jvs(134) * jvs(205) - jvs(191) * jvs(210)
3773	! i = 28
3774	jvs(217) = (jvs(217)) / jvs(3)
3775	jvs(218) = (jvs(218)) / jvs(14)
3776	jvs(219) = (jvs(219)) / jvs(19)
3777	jvs(220) = (jvs(220)) / jvs(34)
3778	jvs(221) = (jvs(221)) / jvs(40)
3779	jvs(222) = (jvs(222)) / jvs(51)
3780	jvs(223) = (jvs(223)) / jvs(58)
3781	jvs(224) = (jvs(224)) / jvs(63)
3782	jvs(225) = (jvs(225)) / jvs(69)
3783	a = 0.0; a = a - jvs(70) * jvs(225)
3784	jvs(226) = (jvs(226) + a) / jvs(80)
3785	a = 0.0; a = a - jvs(64) * jvs(224) - jvs(71) * jvs(225)
3786	jvs(227) = (jvs(227) + a) / jvs(107)
3787	a = 0.0; a = a - jvs(52) * jvs(222)
3788	jvs(228) = (jvs(228) + a) / jvs(114)
3789	a = 0.0; a = a - jvs(72) * jvs(225)
3790	jvs(229) = (jvs(229) + a) / jvs(124)
3791	a = 0.0; a = a - jvs(65) * jvs(224) - jvs(73) * jvs(225) - jvs(81) * jvs(226) - jvs(115) * jvs(228)&
3792	- jvs(125) * jvs(229)
3793	jvs(230) = (jvs(230) + a) / jvs(135)
3794	a = 0.0; a = a - jvs(74) * jvs(225) - jvs(116) * jvs(228) - jvs(126) * jvs(229)
3795	jvs(231) = (jvs(231) + a) / jvs(141)
3796	a = 0.0; a = a - jvs(75) * jvs(225) - jvs(127) * jvs(229)
3797	jvs(232) = (jvs(232) + a) / jvs(152)
3798	a = 0.0; a = a - jvs(66) * jvs(224) - jvs(76) * jvs(225) - jvs(82) * jvs(226) - jvs(108) * jvs(227)&
3799	- jvs(117) * jvs(228) - jvs(128) * jvs(229) - jvs(136) * jvs(230) - jvs(142) * jvs(231) - jvs(153) * jvs(232)
3800	jvs(233) = (jvs(233) + a) / jvs(163)
3801	a = 0.0; a = a - jvs(35) * jvs(220) - jvs(53) * jvs(222) - jvs(118) * jvs(228) - jvs(154) * jvs(232)&
3802	- jvs(164) * jvs(233)
3803	jvs(234) = (jvs(234) + a) / jvs(185)
3804	a = 0.0; a = a - jvs(4) * jvs(217) - jvs(15) * jvs(218) - jvs(20) * jvs(219) - jvs(36) * jvs(220)&
3805	- jvs(41) * jvs(221) - jvs(54) * jvs(222) - jvs(59) * jvs(223) - jvs(67) * jvs(224) - jvs(77) * jvs(225)&
3806	- jvs(83) * jvs(226) - jvs(109) * jvs(227) - jvs(119) * jvs(228) - jvs(129) * jvs(229) - jvs(137) * jvs(230)&
3807	- jvs(143) * jvs(231) - jvs(155) * jvs(232) - jvs(165) * jvs(233) - jvs(186) * jvs(234)
3808	jvs(235) = (jvs(235) + a) / jvs(211)
3809	jvs(236) = jvs(236) - jvs(5) * jvs(217) - jvs(37) * jvs(220) - jvs(130) * jvs(229) - jvs(166) * jvs(233)&
3810	- jvs(187) * jvs(234) - jvs(212) * jvs(235)
3811	jvs(237) = jvs(237) - jvs(78) * jvs(225) - jvs(84) * jvs(226) - jvs(120) * jvs(228) - jvs(131) * jvs(229)&
3812	- jvs(138) * jvs(230) - jvs(144) * jvs(231) - jvs(156) * jvs(232) - jvs(167) * jvs(233) - jvs(188) * jvs(234)&
3813	- jvs(213) * jvs(235)
3814	jvs(238) = jvs(238) - jvs(60) * jvs(223) - jvs(79) * jvs(225) - jvs(110) * jvs(227) - jvs(121) * jvs(228)&
3815	- jvs(132) * jvs(229) - jvs(139) * jvs(230) - jvs(145) * jvs(231) - jvs(157) * jvs(232) - jvs(168) * jvs(233)&
3816	- jvs(189) * jvs(234) - jvs(214) * jvs(235)
3817	jvs(239) = jvs(239) - jvs(42) * jvs(221) - jvs(61) * jvs(223) - jvs(111) * jvs(227) - jvs(133) * jvs(229)&
3818	- jvs(158) * jvs(232) - jvs(169) * jvs(233) - jvs(190) * jvs(234) - jvs(215) * jvs(235)
3819	jvs(240) = jvs(240) - jvs(134) * jvs(229) - jvs(191) * jvs(234) - jvs(216) * jvs(235)
3820	! i = 29
3821	jvs(241) = (jvs(241)) / jvs(1)
3822	jvs(242) = (jvs(242)) / jvs(80)
3823	jvs(243) = (jvs(243)) / jvs(124)
3824	a = 0.0; a = a - jvs(81) * jvs(242) - jvs(125) * jvs(243)
3825	jvs(244) = (jvs(244) + a) / jvs(135)
3826	a = 0.0; a = a - jvs(126) * jvs(243)
3827	jvs(245) = (jvs(245) + a) / jvs(141)
3828	a = 0.0; a = a - jvs(127) * jvs(243)
3829	jvs(246) = (jvs(246) + a) / jvs(152)
3830	a = 0.0; a = a - jvs(2) * jvs(241) - jvs(82) * jvs(242) - jvs(128) * jvs(243) - jvs(136) * jvs(244)&
3831	- jvs(142) * jvs(245) - jvs(153) * jvs(246)
3832	jvs(247) = (jvs(247) + a) / jvs(163)
3833	a = 0.0; a = a - jvs(154) * jvs(246) - jvs(164) * jvs(247)
3834	jvs(248) = (jvs(248) + a) / jvs(185)
3835	a = 0.0; a = a - jvs(83) * jvs(242) - jvs(129) * jvs(243) - jvs(137) * jvs(244) - jvs(143) * jvs(245)&
3836	- jvs(155) * jvs(246) - jvs(165) * jvs(247) - jvs(186) * jvs(248)
3837	jvs(249) = (jvs(249) + a) / jvs(211)
3838	a = 0.0; a = a - jvs(130) * jvs(243) - jvs(166) * jvs(247) - jvs(187) * jvs(248) - jvs(212) * jvs(249)
3839	jvs(250) = (jvs(250) + a) / jvs(236)
3840	jvs(251) = jvs(251) - jvs(84) * jvs(242) - jvs(131) * jvs(243) - jvs(138) * jvs(244) - jvs(144) * jvs(245)&
3841	- jvs(156) * jvs(246) - jvs(167) * jvs(247) - jvs(188) * jvs(248) - jvs(213) * jvs(249) - jvs(237) * jvs(250)
3842	jvs(252) = jvs(252) - jvs(132) * jvs(243) - jvs(139) * jvs(244) - jvs(145) * jvs(245) - jvs(157) * jvs(246)&
3843	- jvs(168) * jvs(247) - jvs(189) * jvs(248) - jvs(214) * jvs(249) - jvs(238) * jvs(250)
3844	jvs(253) = jvs(253) - jvs(133) * jvs(243) - jvs(158) * jvs(246) - jvs(169) * jvs(247) - jvs(190) * jvs(248)&
3845	- jvs(215) * jvs(249) - jvs(239) * jvs(250)
3846	jvs(254) = jvs(254) - jvs(134) * jvs(243) - jvs(191) * jvs(248) - jvs(216) * jvs(249) - jvs(240) * jvs(250)
3847	! i = 30
3848	jvs(255) = (jvs(255)) / jvs(16)
3849	jvs(256) = (jvs(256)) / jvs(44)
3850	a = 0.0; a = a - jvs(45) * jvs(256)
3851	jvs(257) = (jvs(257) + a) / jvs(58)
3852	a = 0.0; a = a - jvs(46) * jvs(256)
3853	jvs(258) = (jvs(258) + a) / jvs(124)
3854	a = 0.0; a = a - jvs(125) * jvs(258)
3855	jvs(259) = (jvs(259) + a) / jvs(135)
3856	a = 0.0; a = a - jvs(126) * jvs(258)
3857	jvs(260) = (jvs(260) + a) / jvs(141)
3858	a = 0.0; a = a - jvs(47) * jvs(256) - jvs(127) * jvs(258)
3859	jvs(261) = (jvs(261) + a) / jvs(152)
3860	a = 0.0; a = a - jvs(128) * jvs(258) - jvs(136) * jvs(259) - jvs(142) * jvs(260) - jvs(153) * jvs(261)
3861	jvs(262) = (jvs(262) + a) / jvs(163)
3862	a = 0.0; a = a - jvs(17) * jvs(255) - jvs(48) * jvs(256) - jvs(154) * jvs(261) - jvs(164) * jvs(262)
3863	jvs(263) = (jvs(263) + a) / jvs(185)
3864	a = 0.0; a = a - jvs(49) * jvs(256) - jvs(59) * jvs(257) - jvs(129) * jvs(258) - jvs(137) * jvs(259)&
3865	- jvs(143) * jvs(260) - jvs(155) * jvs(261) - jvs(165) * jvs(262) - jvs(186) * jvs(263)
3866	jvs(264) = (jvs(264) + a) / jvs(211)
3867	a = 0.0; a = a - jvs(130) * jvs(258) - jvs(166) * jvs(262) - jvs(187) * jvs(263) - jvs(212) * jvs(264)
3868	jvs(265) = (jvs(265) + a) / jvs(236)
3869	a = 0.0; a = a - jvs(131) * jvs(258) - jvs(138) * jvs(259) - jvs(144) * jvs(260) - jvs(156) * jvs(261)&
3870	- jvs(167) * jvs(262) - jvs(188) * jvs(263) - jvs(213) * jvs(264) - jvs(237) * jvs(265)
3871	jvs(266) = (jvs(266) + a) / jvs(251)
3872	jvs(267) = jvs(267) - jvs(18) * jvs(255) - jvs(50) * jvs(256) - jvs(60) * jvs(257) - jvs(132) * jvs(258)&
3873	- jvs(139) * jvs(259) - jvs(145) * jvs(260) - jvs(157) * jvs(261) - jvs(168) * jvs(262) - jvs(189) * jvs(263)&
3874	- jvs(214) * jvs(264) - jvs(238) * jvs(265) - jvs(252) * jvs(266)
3875	jvs(268) = jvs(268) - jvs(61) * jvs(257) - jvs(133) * jvs(258) - jvs(158) * jvs(261) - jvs(169) * jvs(262)&
3876	- jvs(190) * jvs(263) - jvs(215) * jvs(264) - jvs(239) * jvs(265) - jvs(253) * jvs(266)
3877	jvs(269) = jvs(269) - jvs(134) * jvs(258) - jvs(191) * jvs(263) - jvs(216) * jvs(264) - jvs(240) * jvs(265)&
3878	- jvs(254) * jvs(266)
3879	! i = 31
3880	jvs(270) = (jvs(270)) / jvs(21)
3881	jvs(271) = (jvs(271)) / jvs(30)
3882	jvs(272) = (jvs(272)) / jvs(40)
3883	a = 0.0; a = a - jvs(22) * jvs(270)
3884	jvs(273) = (jvs(273) + a) / jvs(51)
3885	jvs(274) = (jvs(274)) / jvs(91)
3886	a = 0.0; a = a - jvs(92) * jvs(274)
3887	jvs(275) = (jvs(275) + a) / jvs(107)
3888	a = 0.0; a = a - jvs(23) * jvs(270) - jvs(52) * jvs(273) - jvs(93) * jvs(274)
3889	jvs(276) = (jvs(276) + a) / jvs(114)
3890	a = 0.0; a = a - jvs(24) * jvs(270) - jvs(94) * jvs(274) - jvs(115) * jvs(276)
3891	jvs(277) = (jvs(277) + a) / jvs(135)
3892	a = 0.0; a = a - jvs(25) * jvs(270) - jvs(95) * jvs(274) - jvs(116) * jvs(276)
3893	jvs(278) = (jvs(278) + a) / jvs(141)
3894	a = 0.0; a = a - jvs(96) * jvs(274)
3895	jvs(279) = (jvs(279) + a) / jvs(152)
3896	a = 0.0; a = a - jvs(97) * jvs(274) - jvs(108) * jvs(275) - jvs(117) * jvs(276) - jvs(136) * jvs(277)&
3897	- jvs(142) * jvs(278) - jvs(153) * jvs(279)
3898	jvs(280) = (jvs(280) + a) / jvs(163)
3899	a = 0.0; a = a - jvs(31) * jvs(271) - jvs(53) * jvs(273) - jvs(98) * jvs(274) - jvs(118) * jvs(276)&
3900	- jvs(154) * jvs(279) - jvs(164) * jvs(280)
3901	jvs(281) = (jvs(281) + a) / jvs(185)
3902	a = 0.0; a = a - jvs(26) * jvs(270) - jvs(32) * jvs(271) - jvs(41) * jvs(272) - jvs(54) * jvs(273)&
3903	- jvs(99) * jvs(274) - jvs(109) * jvs(275) - jvs(119) * jvs(276) - jvs(137) * jvs(277) - jvs(143) * jvs(278)&
3904	- jvs(155) * jvs(279) - jvs(165) * jvs(280) - jvs(186) * jvs(281)
3905	jvs(282) = (jvs(282) + a) / jvs(211)
3906	a = 0.0; a = a - jvs(100) * jvs(274) - jvs(166) * jvs(280) - jvs(187) * jvs(281) - jvs(212) * jvs(282)
3907	jvs(283) = (jvs(283) + a) / jvs(236)
3908	a = 0.0; a = a - jvs(27) * jvs(270) - jvs(101) * jvs(274) - jvs(120) * jvs(276) - jvs(138) * jvs(277)&
3909	- jvs(144) * jvs(278) - jvs(156) * jvs(279) - jvs(167) * jvs(280) - jvs(188) * jvs(281) - jvs(213) * jvs(282)&
3910	- jvs(237) * jvs(283)
3911	jvs(284) = (jvs(284) + a) / jvs(251)
3912	a = 0.0; a = a - jvs(28) * jvs(270) - jvs(102) * jvs(274) - jvs(110) * jvs(275) - jvs(121) * jvs(276)&
3913	- jvs(139) * jvs(277) - jvs(145) * jvs(278) - jvs(157) * jvs(279) - jvs(168) * jvs(280) - jvs(189) * jvs(281)&
3914	- jvs(214) * jvs(282) - jvs(238) * jvs(283) - jvs(252) * jvs(284)
3915	jvs(285) = (jvs(285) + a) / jvs(267)
3916	jvs(286) = jvs(286) - jvs(29) * jvs(270) - jvs(33) * jvs(271) - jvs(42) * jvs(272) - jvs(103) * jvs(274)&
3917	- jvs(111) * jvs(275) - jvs(158) * jvs(279) - jvs(169) * jvs(280) - jvs(190) * jvs(281) - jvs(215) * jvs(282)&
3918	- jvs(239) * jvs(283) - jvs(253) * jvs(284) - jvs(268) * jvs(285)
3919	jvs(287) = jvs(287) - jvs(104) * jvs(274) - jvs(191) * jvs(281) - jvs(216) * jvs(282) - jvs(240) * jvs(283)&
3920	- jvs(254) * jvs(284) - jvs(269) * jvs(285)
3921	! i = 32
3922	jvs(288) = (jvs(288)) / jvs(6)
3923	jvs(289) = (jvs(289)) / jvs(63)
3924	a = 0.0; a = a - jvs(64) * jvs(289)
3925	jvs(290) = (jvs(290) + a) / jvs(107)
3926	a = 0.0; a = a - jvs(65) * jvs(289)
3927	jvs(291) = (jvs(291) + a) / jvs(135)
3928	jvs(292) = (jvs(292)) / jvs(152)
3929	a = 0.0; a = a - jvs(66) * jvs(289) - jvs(108) * jvs(290) - jvs(136) * jvs(291) - jvs(153) * jvs(292)
3930	jvs(293) = (jvs(293) + a) / jvs(163)
3931	a = 0.0; a = a - jvs(7) * jvs(288) - jvs(154) * jvs(292) - jvs(164) * jvs(293)
3932	jvs(294) = (jvs(294) + a) / jvs(185)
3933	a = 0.0; a = a - jvs(67) * jvs(289) - jvs(109) * jvs(290) - jvs(137) * jvs(291) - jvs(155) * jvs(292)&
3934	- jvs(165) * jvs(293) - jvs(186) * jvs(294)
3935	jvs(295) = (jvs(295) + a) / jvs(211)
3936	a = 0.0; a = a - jvs(166) * jvs(293) - jvs(187) * jvs(294) - jvs(212) * jvs(295)
3937	jvs(296) = (jvs(296) + a) / jvs(236)
3938	a = 0.0; a = a - jvs(138) * jvs(291) - jvs(156) * jvs(292) - jvs(167) * jvs(293) - jvs(188) * jvs(294)&
3939	- jvs(213) * jvs(295) - jvs(237) * jvs(296)
3940	jvs(297) = (jvs(297) + a) / jvs(251)
3941	a = 0.0; a = a - jvs(110) * jvs(290) - jvs(139) * jvs(291) - jvs(157) * jvs(292) - jvs(168) * jvs(293)&
3942	- jvs(189) * jvs(294) - jvs(214) * jvs(295) - jvs(238) * jvs(296) - jvs(252) * jvs(297)
3943	jvs(298) = (jvs(298) + a) / jvs(267)
3944	a = 0.0; a = a - jvs(111) * jvs(290) - jvs(158) * jvs(292) - jvs(169) * jvs(293) - jvs(190) * jvs(294)&
3945	- jvs(215) * jvs(295) - jvs(239) * jvs(296) - jvs(253) * jvs(297) - jvs(268) * jvs(298)
3946	jvs(299) = (jvs(299) + a) / jvs(286)
3947	jvs(300) = jvs(300) - jvs(8) * jvs(288) - jvs(191) * jvs(294) - jvs(216) * jvs(295) - jvs(240) * jvs(296)&
3948	- jvs(254) * jvs(297) - jvs(269) * jvs(298) - jvs(287) * jvs(299)
3949	RETURN
3950
3951	END SUBROUTINE kppdecomp
3952
3953	SUBROUTINE get_mechanismname
3954
3955	IMPLICIT NONE
3956
3957	! Set cs_mech for check with mechanism name from namelist
3958	cs_mech = 'cbm4'
3959
3960	RETURN
3961	END SUBROUTINE get_mechanismname
3962
3963
3964	SUBROUTINE chem_gasphase_integrate (time_step_len, conc, tempi, qvapi, fakti, photo, ierrf, xnacc, xnrej, istatus, l_debug, pe, &
3965	icntrl_i, rcntrl_i)
3966
3967	IMPLICIT NONE
3968
3969	REAL(dp), INTENT(IN) :: time_step_len
3970	REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: conc
3971	REAL(dp), DIMENSION(:, :), INTENT(IN) :: photo
3972	REAL(dp), DIMENSION(:), INTENT(IN) :: tempi
3973	REAL(dp), DIMENSION(:), INTENT(IN) :: qvapi
3974	REAL(dp), DIMENSION(:), INTENT(IN) :: fakti
3975	INTEGER, INTENT(OUT), OPTIONAL :: ierrf(:)
3976	INTEGER, INTENT(OUT), OPTIONAL :: xnacc(:)
3977	INTEGER, INTENT(OUT), OPTIONAL :: xnrej(:)
3978	INTEGER, INTENT(INOUT), OPTIONAL :: istatus(:)
3979	INTEGER, INTENT(IN), OPTIONAL :: pe
3980	LOGICAL, INTENT(IN), OPTIONAL :: l_debug
3981	INTEGER, DIMENSION(nkppctrl), INTENT(IN), OPTIONAL :: icntrl_i
3982	REAL(dp), DIMENSION(nkppctrl), INTENT(IN), OPTIONAL :: rcntrl_i
3983
3984	INTEGER :: k ! loop variable
3985	REAL(dp) :: dt
3986	INTEGER, DIMENSION(20) :: istatus_u
3987	INTEGER :: ierr_u
3988	INTEGER :: istatf
3989	INTEGER :: vl_dim_lo
3990
3991
3992	IF (PRESENT (istatus)) istatus = 0
3993	IF (PRESENT (icntrl_i)) icntrl = icntrl_i
3994	IF (PRESENT (rcntrl_i)) rcntrl = rcntrl_i
3995
3996	vl_glo = size(tempi, 1)
3997
3998	vl_dim_lo = vl_dim
3999	DO k=1, vl_glo, vl_dim_lo
4000	is = k
4001	ie = min(k+ vl_dim_lo-1, vl_glo)
4002	vl = ie-is+ 1
4003
4004	c(:) = conc(is, :)
4005
4006	temp = tempi(is)
4007
4008	qvap = qvapi(is)
4009
4010	fakt = fakti(is)
4011
4012	CALL initialize
4013
4014	phot(:) = photo(is, :)
4015
4016	CALL update_rconst
4017
4018	dt = time_step_len
4019
4020	! integrate from t=0 to t=dt
4021	CALL integrate(0._dp, dt, icntrl, rcntrl, istatus_u = istatus_u, ierr_u=ierr_u)
4022
4023
4024	IF (PRESENT(l_debug) .AND. PRESENT(pe)) THEN
4025	IF (l_debug) CALL error_output(conc(is, :), ierr_u, pe)
4026	ENDIF
4027
4028	conc(is, :) = c(:)
4029
4030	! RETURN diagnostic information
4031
4032	IF (PRESENT(ierrf)) ierrf(is) = ierr_u
4033	IF (PRESENT(xnacc)) xnacc(is) = istatus_u(4)
4034	IF (PRESENT(xnrej)) xnrej(is) = istatus_u(5)
4035
4036	IF (PRESENT (istatus)) THEN
4037	istatus(1:8) = istatus(1:8) + istatus_u(1:8)
4038	ENDIF
4039
4040	END DO
4041
4042
4043	! Deallocate input arrays
4044
4045
4046	data_loaded = .FALSE.
4047
4048	RETURN
4049	END SUBROUTINE chem_gasphase_integrate
4050
4051	END MODULE chem_gasphase_mod
4052

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |