Home

Context Navigation

chem_gasphase_mod.f90 @ 3678

Last change on this file since 3678 was 3623, checked in by forkel, 6 years ago
removed inlined declaration of qvap and fakt - was forgotten for CBM4
File size: 153.9 KB

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