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

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

Merge chemistry branch at r3297 to trunk

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