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

Last change on this file since 4370 was 4370, checked in by raasch, 2 years ago

bugfixes for previous commit: unused variables removed, Temperton-fft usage on GPU, openacc porting of vector version of Obukhov length calculation, collective read switched off on NEC to avoid hanging; some vector directives added in prognostic equations to force vectorization on Intel19 compiler, configuration files for NEC Aurora added

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