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

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

Removed unused variables from chem_gasphase_mod.f90

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