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

Last change on this file since 3799 was 3799, checked in by forkel, 3 years ago

editing in kpp4palm: add statements for avoiding unused variables, remove $Id

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