Home

Context Navigation

chem_gasphase_mod.f90 @ 3800

Last change on this file since 3800 was 3800, checked in by forkel, 4 years ago
added leading blanks to dummy statements
File size: 155.4 KB

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