Home

Context Navigation

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

Last change on this file since 3944 was 3944, checked in by maronga, 6 years ago
fix for last commit
File size: 155.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |