Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |