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

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

removed inlined declaration of qvap and fakt - was forgotten for CBM4

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