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

Last change on this file since 3458 was 3458, checked in by kanani, 3 years ago

Reintegrated fixes/changes from branch chemistry

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