Home

Context Navigation

chem_gasphase_mod.f90 @ 3298

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