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

Last change on this file since 4841 was 4841, checked in by forkel, 21 months 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 
1MODULE 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 
871SUBROUTINE 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     
904END SUBROUTINE initialize
905 
906SUBROUTINE 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
957END SUBROUTINE integrate
958 
959SUBROUTINE 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     
1113END SUBROUTINE fun
1114 
1115SUBROUTINE 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     
1211END SUBROUTINE kppsolve
1212 
1213SUBROUTINE 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     
2107END 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 
2122SUBROUTINE 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     
2214END SUBROUTINE update_rconst
2215 
2216!  END FUNCTION ARR2
2217REAL(kind=dp)FUNCTION arr2( a0, b0, temp)
2218    REAL(kind=dp):: temp
2219    REAL(kind=dp):: a0, b0
2220    arr2 = a0 * exp( - b0 / temp)
2221END FUNCTION arr2
2222 
2223SUBROUTINE 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
2271END SUBROUTINE initialize_kpp_ctrl
2272 
2273SUBROUTINE 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
2283END 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 
2386SUBROUTINE 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!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2674CONTAINS !  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 + 10*H = 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
2787timeloop: 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
2817untilaccepted: 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
2827stage: 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!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3532END SUBROUTINE rosenbrock
3533 
3534SUBROUTINE 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
3551END SUBROUTINE funtemplate
3552 
3553SUBROUTINE 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
3589END 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