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

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

removed read from unit 10 in chemistry_model_mod.f90, added get_mechanismname

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