Home

Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

chem_gasphase_mod.f90

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