source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_salsa+simple/chem_gasphase_mod.f90 @ 3799

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

editing in kpp4palm: add statements for avoiding unused variables, remove $Id

File size: 88.0 KB
Line 
1MODULE chem_gasphase_mod
2 
3!   Mechanism: salsa+simple
4!
5!------------------------------------------------------------------------------!
6!
7! ******Module chem_gasphase_mod is automatically generated by kpp4palm ******
8!
9!   *********Please do NOT change this Code,it will be ovewritten *********
10!
11!------------------------------------------------------------------------------!
12! This file was created by KPP (http://people.cs.vt.edu/asandu/Software/Kpp/)
13! and kpp4palm (created by Klaus Ketelsen). kpp4palm is an adapted version
14! of KP4 (Jöckel,P.,Kerkweg,A.,Pozzer,A.,Sander,R.,Tost,H.,Riede,
15! H.,Baumgaertner,A.,Gromov,S.,and Kern,B.,2010: Development cycle 2 of
16! the Modular Earth Submodel System (MESSy2),Geosci. Model Dev.,3,717-752,
17! https://doi.org/10.5194/gmd-3-717-2010). KP4 is part of the Modular Earth
18! Submodel System (MESSy),which is is available under the  GNU General Public
19! License (GPL).
20!
21! KPP is free software; you can redistribute it and/or modify it under the terms
22! of the General Public Licence as published by the Free Software Foundation;
23! either version 2 of the License,or (at your option) any later version.
24! KPP is distributed in the hope that it will be useful,but WITHOUT ANY WARRANTY;
25! without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
26! PURPOSE. See the GNU General Public Licence for more details.
27!
28!------------------------------------------------------------------------------!
29! This file is part of the PALM model system.
30!
31! PALM is free software: you can redistribute it and/or modify it under the
32! terms of the GNU General Public License as published by the Free Software
33! Foundation,either version 3 of the License,or (at your option) any later
34! version.
35!
36! PALM is distributed in the hope that it will be useful,but WITHOUT ANY
37! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
38! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
39!
40! You should have received a copy of the GNU General Public License along with
41! PALM. If not,see <http://www.gnu.org/licenses/>.
42!
43! Copyright 1997-2019 Leibniz Universitaet Hannover
44!--------------------------------------------------------------------------------!
45!
46!
47! MODULE HEADER TEMPLATE
48!
49!  Initial version (Nov. 2016,ketelsen),for later modifications of module_header
50!  see comments in kpp4palm/src/create_kpp_module.C
51
52! Set kpp Double Precision to PALM Default Precision
53
54  USE kinds,           ONLY: dp=>wp
55
56  USE pegrid,          ONLY: myid, threads_per_task
57
58  IMPLICIT        NONE
59  PRIVATE
60  !SAVE  ! note: occurs again in automatically generated code ...
61
62!  PUBLIC :: IERR_NAMES
63 
64! PUBLIC :: SPC_NAMES,EQN_NAMES,EQN_TAGS,REQ_HET,REQ_AEROSOL,REQ_PHOTRAT &
65!         ,REQ_MCFCT,IP_MAX,jname
66
67  PUBLIC :: cs_mech
68  PUBLIC :: eqn_names, phot_names, spc_names
69  PUBLIC :: nmaxfixsteps
70  PUBLIC :: atol, rtol
71  PUBLIC :: nspec, nreact
72  PUBLIC :: temp
73  PUBLIC :: qvap
74  PUBLIC :: fakt
75  PUBLIC :: phot
76  PUBLIC :: rconst
77  PUBLIC :: nvar
78  PUBLIC :: nphot
79  PUBLIC :: vl_dim                     ! PUBLIC to ebable other MODULEs to distiguish between scalar and vec
80 
81  PUBLIC :: initialize, integrate, update_rconst
82  PUBLIC :: chem_gasphase_integrate
83  PUBLIC :: initialize_kpp_ctrl
84  PUBLIC :: get_mechanismname
85
86! END OF MODULE HEADER TEMPLATE
87                                                                 
88! Variables used for vector mode                                 
89                                                                 
90  LOGICAL, PARAMETER          :: l_vector = .FALSE.           
91  INTEGER, PARAMETER          :: i_lu_di = 2
92  INTEGER, PARAMETER          :: vl_dim = 1
93  INTEGER                     :: vl                             
94                                                                 
95  INTEGER                     :: vl_glo                         
96  INTEGER                     :: is, ie                           
97                                                                 
98                                                                 
99  LOGICAL                     :: data_loaded = .FALSE.             
100  REAL(dp), POINTER, DIMENSION(:), CONTIGUOUS    :: var             
101! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
102!
103! Parameter Module File
104!
105! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
106!       (http://www.cs.vt.edu/~asandu/Software/KPP)
107! KPP is distributed under GPL,the general public licence
108!       (http://www.gnu.org/copyleft/gpl.html)
109! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
110! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
111!     With important contributions from:
112!        M. Damian,Villanova University,USA
113!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
114!
115! File                 : chem_gasphase_mod_Parameters.f90
116! Time                 : Fri Mar 15 16:28:26 2019
117! Working directory    : /home/forkel-r/palmstuff/work/trunk20190315/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
118! Equation file        : chem_gasphase_mod.kpp
119! Output root filename : chem_gasphase_mod
120!
121! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
122
123
124
125
126
127
128! NSPEC - Number of chemical species
129  INTEGER, PARAMETER :: nspec = 14 
130! NVAR - Number of Variable species
131  INTEGER, PARAMETER :: nvar = 13 
132! NVARACT - Number of Active species
133  INTEGER, PARAMETER :: nvaract = 11 
134! NFIX - Number of Fixed species
135  INTEGER, PARAMETER :: nfix = 1 
136! NREACT - Number of reactions
137  INTEGER, PARAMETER :: nreact = 11 
138! NVARST - Starting of variables in conc. vect.
139  INTEGER, PARAMETER :: nvarst = 1 
140! NFIXST - Starting of fixed in conc. vect.
141  INTEGER, PARAMETER :: nfixst = 14 
142! NONZERO - Number of nonzero entries in Jacobian
143  INTEGER, PARAMETER :: nonzero = 39 
144! LU_NONZERO - Number of nonzero entries in LU factoriz. of Jacobian
145  INTEGER, PARAMETER :: lu_nonzero = 41 
146! CNVAR - (NVAR+1) Number of elements in compressed row format
147  INTEGER, PARAMETER :: cnvar = 14 
148! CNEQN - (NREACT+1) Number stoicm elements in compressed col format
149  INTEGER, PARAMETER :: cneqn = 12 
150! NHESS - Length of Sparse Hessian
151  INTEGER, PARAMETER :: nhess = 18 
152! NMASS - Number of atoms to check mass balance
153  INTEGER, PARAMETER :: nmass = 1 
154
155! Index declaration for variable species in C and VAR
156!   VAR(ind_spc) = C(ind_spc)
157
158  INTEGER, PARAMETER, PUBLIC :: ind_h2so4 = 1 
159  INTEGER, PARAMETER, PUBLIC :: ind_nh3 = 2 
160  INTEGER, PARAMETER, PUBLIC :: ind_ocnv = 3 
161  INTEGER, PARAMETER, PUBLIC :: ind_ocsv = 4 
162  INTEGER, PARAMETER, PUBLIC :: ind_hno3 = 5 
163  INTEGER, PARAMETER, PUBLIC :: ind_rcho = 6 
164  INTEGER, PARAMETER, PUBLIC :: ind_rh = 7 
165  INTEGER, PARAMETER, PUBLIC :: ind_o3 = 8 
166  INTEGER, PARAMETER, PUBLIC :: ind_oh = 9 
167  INTEGER, PARAMETER, PUBLIC :: ind_no2 = 10 
168  INTEGER, PARAMETER, PUBLIC :: ind_no = 11 
169  INTEGER, PARAMETER, PUBLIC :: ind_ho2 = 12 
170  INTEGER, PARAMETER, PUBLIC :: ind_ro2 = 13 
171
172! Index declaration for fixed species in C
173!   C(ind_spc)
174
175  INTEGER, PARAMETER, PUBLIC :: ind_h2o = 14 
176
177! Index declaration for fixed species in FIX
178!    FIX(indf_spc) = C(ind_spc) = C(NVAR+indf_spc)
179
180  INTEGER, PARAMETER :: indf_h2o = 1 
181
182! NJVRP - Length of sparse Jacobian JVRP
183  INTEGER, PARAMETER :: njvrp = 16 
184
185! NSTOICM - Length of Sparse Stoichiometric Matrix
186  INTEGER, PARAMETER :: nstoicm = 23 
187
188
189! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
190!
191! Global Data Module File
192!
193! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
194!       (http://www.cs.vt.edu/~asandu/Software/KPP)
195! KPP is distributed under GPL,the general public licence
196!       (http://www.gnu.org/copyleft/gpl.html)
197! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
198! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
199!     With important contributions from:
200!        M. Damian,Villanova University,USA
201!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
202!
203! File                 : chem_gasphase_mod_Global.f90
204! Time                 : Fri Mar 15 16:28:26 2019
205! Working directory    : /home/forkel-r/palmstuff/work/trunk20190315/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
206! Equation file        : chem_gasphase_mod.kpp
207! Output root filename : chem_gasphase_mod
208!
209! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
210
211
212
213
214
215
216! Declaration of global variables
217
218! C - Concentration of all species
219  REAL(kind=dp), TARGET    :: c(nspec)
220! VAR - Concentrations of variable species (global)
221! REAL(kind=dp):: var(nvar)  var is now POINTER
222! FIX - Concentrations of fixed species (global)
223  REAL(kind=dp):: fix(nfix)
224! VAR,FIX are chunks of array C
225! RCONST - Rate constants (global)
226  REAL(kind=dp):: rconst(nreact)
227! TIME - Current integration time
228  REAL(kind=dp):: time
229! TEMP - Temperature
230  REAL(kind=dp):: temp
231! ATOL - Absolute tolerance
232  REAL(kind=dp):: atol(nvar)
233! RTOL - Relative tolerance
234  REAL(kind=dp):: rtol(nvar)
235! STEPMIN - Lower bound for integration step
236  REAL(kind=dp):: stepmin
237! CFACTOR - Conversion factor for concentration units
238  REAL(kind=dp):: cfactor
239
240! INLINED global variable declarations
241
242! QVAP - Water vapor
243  REAL(dp)                                    :: qvap
244! FAKT - Conversion factor
245  REAL(dp)                                    :: fakt
246!   Declaration of global variable declarations for photolysis will come from IN
247
248! QVAP - Water vapor
249  REAL(kind=dp):: qvap
250! FAKT - Conversion factor
251  REAL(kind=dp):: fakt
252
253! CS_MECH for check of mechanism name with namelist
254  CHARACTER(len=30):: cs_mech
255
256! INLINED global variable declarations
257
258
259
260! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
261!
262! Sparse Jacobian Data Structures File
263!
264! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
265!       (http://www.cs.vt.edu/~asandu/Software/KPP)
266! KPP is distributed under GPL,the general public licence
267!       (http://www.gnu.org/copyleft/gpl.html)
268! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
269! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
270!     With important contributions from:
271!        M. Damian,Villanova University,USA
272!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
273!
274! File                 : chem_gasphase_mod_JacobianSP.f90
275! Time                 : Fri Mar 15 16:28:26 2019
276! Working directory    : /home/forkel-r/palmstuff/work/trunk20190315/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
277! Equation file        : chem_gasphase_mod.kpp
278! Output root filename : chem_gasphase_mod
279!
280! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
281
282
283
284
285
286
287! Sparse Jacobian Data
288
289
290  INTEGER, PARAMETER, DIMENSION(41):: lu_irow =  (/ &
291       1, 2, 3, 4, 5, 5, 5, 6, 6, 6, 7, 7, &
292       8, 8, 8, 9, 9, 9, 9, 9, 9, 10, 10, 10, &
293      10, 10, 10, 11, 11, 11, 11, 11, 12, 12, 12, 13, &
294      13, 13, 13, 13, 13 /) 
295
296  INTEGER, PARAMETER, DIMENSION(41):: lu_icol =  (/ &
297       1, 2, 3, 4, 5, 9, 10, 6, 11, 13, 7, 9, &
298       8, 10, 11, 7, 8, 9, 10, 11, 12, 8, 9, 10, &
299      11, 12, 13, 8, 10, 11, 12, 13, 11, 12, 13, 7, &
300       9, 10, 11, 12, 13 /) 
301
302  INTEGER, PARAMETER, DIMENSION(14):: lu_crow =  (/ &
303       1, 2, 3, 4, 5, 8, 11, 13, 16, 22, 28, 33, &
304      36, 42 /) 
305
306  INTEGER, PARAMETER, DIMENSION(14):: lu_diag =  (/ &
307       1, 2, 3, 4, 5, 8, 11, 13, 18, 24, 30, 34, &
308      41, 42 /) 
309
310
311
312! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
313!
314! Utility Data Module File
315!
316! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
317!       (http://www.cs.vt.edu/~asandu/Software/KPP)
318! KPP is distributed under GPL,the general public licence
319!       (http://www.gnu.org/copyleft/gpl.html)
320! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
321! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
322!     With important contributions from:
323!        M. Damian,Villanova University,USA
324!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
325!
326! File                 : chem_gasphase_mod_Monitor.f90
327! Time                 : Fri Mar 15 16:28:26 2019
328! Working directory    : /home/forkel-r/palmstuff/work/trunk20190315/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
329! Equation file        : chem_gasphase_mod.kpp
330! Output root filename : chem_gasphase_mod
331!
332! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
333
334
335
336
337
338  CHARACTER(len=15), PARAMETER, DIMENSION(14):: spc_names =  (/ &
339     'H2SO4          ','NH3            ','OCNV           ',&
340     'OCSV           ','HNO3           ','RCHO           ',&
341     'RH             ','O3             ','OH             ',&
342     'NO2            ','NO             ','HO2            ',&
343     'RO2            ','H2O            ' /)
344
345  CHARACTER(len=100), PARAMETER, DIMENSION(11):: eqn_names =  (/ &
346     '     NO2 --> O3 + NO                                                                                ',&
347     'O3 + H2O --> 2 OH                                                                                   ',&
348     ' O3 + NO --> NO2                                                                                    ',&
349     ' RH + OH --> RO2 + H2O                                                                              ',&
350     'NO + RO2 --> RCHO + NO2 + HO2                                                                       ',&
351     'NO + HO2 --> OH + NO2                                                                               ',&
352     'OH + NO2 --> HNO3                                                                                   ',&
353     '   H2SO4 --> H2SO4                                                                                  ',&
354     '     NH3 --> NH3                                                                                    ',&
355     '    OCNV --> OCNV                                                                                   ',&
356     '    OCSV --> OCSV                                                                                   ' /)
357
358! INLINED global variables
359
360  !   inline f90_data: declaration of global variables for photolysis
361  !   REAL(kind=dp):: phot(nphot)must eventually be moved to global later for
362  INTEGER, PARAMETER :: nphot = 2
363  !   phot photolysis frequencies
364  REAL(kind=dp):: phot(nphot)
365
366  INTEGER, PARAMETER, PUBLIC :: j_no2 = 1
367  INTEGER, PARAMETER, PUBLIC :: j_o31d = 2
368
369  CHARACTER(len=15), PARAMETER, DIMENSION(nphot):: phot_names =   (/ &
370     'J_NO2          ','J_O31D         '/)
371
372! End INLINED global variables
373
374
375! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
376 
377! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
378 
379! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
380 
381! Automatic generated PUBLIC Statements for ip_ and ihs_ variables
382 
383 
384!  variable definations from  individual module headers
385 
386! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
387!
388! Initialization File
389!
390! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
391!       (http://www.cs.vt.edu/~asandu/Software/KPP)
392! KPP is distributed under GPL,the general public licence
393!       (http://www.gnu.org/copyleft/gpl.html)
394! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
395! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
396!     With important contributions from:
397!        M. Damian,Villanova University,USA
398!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
399!
400! File                 : chem_gasphase_mod_Initialize.f90
401! Time                 : Fri Mar 15 16:28:26 2019
402! Working directory    : /home/forkel-r/palmstuff/work/trunk20190315/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
403! Equation file        : chem_gasphase_mod.kpp
404! Output root filename : chem_gasphase_mod
405!
406! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
407
408
409
410
411
412! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
413!
414! Numerical Integrator (Time-Stepping) File
415!
416! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
417!       (http://www.cs.vt.edu/~asandu/Software/KPP)
418! KPP is distributed under GPL,the general public licence
419!       (http://www.gnu.org/copyleft/gpl.html)
420! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
421! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
422!     With important contributions from:
423!        M. Damian,Villanova University,USA
424!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
425!
426! File                 : chem_gasphase_mod_Integrator.f90
427! Time                 : Fri Mar 15 16:28:26 2019
428! Working directory    : /home/forkel-r/palmstuff/work/trunk20190315/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
429! Equation file        : chem_gasphase_mod.kpp
430! Output root filename : chem_gasphase_mod
431!
432! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
433
434
435
436! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
437!
438! INTEGRATE - Integrator routine
439!   Arguments :
440!      TIN       - Start Time for Integration
441!      TOUT      - End Time for Integration
442!
443! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
444
445!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
446!  Rosenbrock - Implementation of several Rosenbrock methods:             !
447!               *Ros2                                                    !
448!               *Ros3                                                    !
449!               *Ros4                                                    !
450!               *Rodas3                                                  !
451!               *Rodas4                                                  !
452!  By default the code employs the KPP sparse linear algebra routines     !
453!  Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK)        !
454!                                                                         !
455!    (C)  Adrian Sandu,August 2004                                       !
456!    Virginia Polytechnic Institute and State University                  !
457!    Contact: sandu@cs.vt.edu                                             !
458!    Revised by Philipp Miehe and Adrian Sandu,May 2006                  !                               !
459!    This implementation is part of KPP - the Kinetic PreProcessor        !
460!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
461
462
463  SAVE
464 
465!~~~>  statistics on the work performed by the rosenbrock method
466  INTEGER, PARAMETER :: nfun=1, njac=2, nstp=3, nacc=4, &
467                        nrej=5, ndec=6, nsol=7, nsng=8, &
468                        ntexit=1, nhexit=2, nhnew = 3
469
470! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
471!
472! Linear Algebra Data and Routines File
473!
474! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
475!       (http://www.cs.vt.edu/~asandu/Software/KPP)
476! KPP is distributed under GPL,the general public licence
477!       (http://www.gnu.org/copyleft/gpl.html)
478! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
479! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
480!     With important contributions from:
481!        M. Damian,Villanova University,USA
482!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
483!
484! File                 : chem_gasphase_mod_LinearAlgebra.f90
485! Time                 : Fri Mar 15 16:28:26 2019
486! Working directory    : /home/forkel-r/palmstuff/work/trunk20190315/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
487! Equation file        : chem_gasphase_mod.kpp
488! Output root filename : chem_gasphase_mod
489!
490! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
491
492
493
494
495
496
497! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
498!
499! The ODE Jacobian of Chemical Model File
500!
501! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
502!       (http://www.cs.vt.edu/~asandu/Software/KPP)
503! KPP is distributed under GPL,the general public licence
504!       (http://www.gnu.org/copyleft/gpl.html)
505! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
506! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
507!     With important contributions from:
508!        M. Damian,Villanova University,USA
509!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
510!
511! File                 : chem_gasphase_mod_Jacobian.f90
512! Time                 : Fri Mar 15 16:28:26 2019
513! Working directory    : /home/forkel-r/palmstuff/work/trunk20190315/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
514! Equation file        : chem_gasphase_mod.kpp
515! Output root filename : chem_gasphase_mod
516!
517! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
518
519
520
521
522
523
524! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
525!
526! The ODE Function of Chemical Model File
527!
528! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
529!       (http://www.cs.vt.edu/~asandu/Software/KPP)
530! KPP is distributed under GPL,the general public licence
531!       (http://www.gnu.org/copyleft/gpl.html)
532! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
533! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
534!     With important contributions from:
535!        M. Damian,Villanova University,USA
536!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
537!
538! File                 : chem_gasphase_mod_Function.f90
539! Time                 : Fri Mar 15 16:28:26 2019
540! Working directory    : /home/forkel-r/palmstuff/work/trunk20190315/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
541! Equation file        : chem_gasphase_mod.kpp
542! Output root filename : chem_gasphase_mod
543!
544! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
545
546
547
548
549
550! A - Rate for each equation
551  REAL(kind=dp):: a(nreact)
552
553! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
554!
555! The Reaction Rates File
556!
557! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
558!       (http://www.cs.vt.edu/~asandu/Software/KPP)
559! KPP is distributed under GPL,the general public licence
560!       (http://www.gnu.org/copyleft/gpl.html)
561! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
562! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
563!     With important contributions from:
564!        M. Damian,Villanova University,USA
565!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
566!
567! File                 : chem_gasphase_mod_Rates.f90
568! Time                 : Fri Mar 15 16:28:26 2019
569! Working directory    : /home/forkel-r/palmstuff/work/trunk20190315/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
570! Equation file        : chem_gasphase_mod.kpp
571! Output root filename : chem_gasphase_mod
572!
573! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
574
575
576
577
578
579! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
580!
581! Auxiliary Routines File
582!
583! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
584!       (http://www.cs.vt.edu/~asandu/Software/KPP)
585! KPP is distributed under GPL,the general public licence
586!       (http://www.gnu.org/copyleft/gpl.html)
587! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
588! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
589!     With important contributions from:
590!        M. Damian,Villanova University,USA
591!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
592!
593! File                 : chem_gasphase_mod_Util.f90
594! Time                 : Fri Mar 15 16:28:26 2019
595! Working directory    : /home/forkel-r/palmstuff/work/trunk20190315/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
596! Equation file        : chem_gasphase_mod.kpp
597! Output root filename : chem_gasphase_mod
598!
599! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
600
601
602
603
604
605
606  ! header MODULE initialize_kpp_ctrl_template
607
608  ! notes:
609  ! - l_vector is automatically defined by kp4
610  ! - vl_dim is automatically defined by kp4
611  ! - i_lu_di is automatically defined by kp4
612  ! - wanted is automatically defined by xmecca
613  ! - icntrl rcntrl are automatically defined by kpp
614  ! - "USE messy_main_tools" is in MODULE_header of messy_mecca_kpp.f90
615  ! - SAVE will be automatically added by kp4
616
617  !SAVE
618
619  ! for fixed time step control
620  ! ... max. number of fixed time steps (sum must be 1)
621  INTEGER, PARAMETER         :: nmaxfixsteps = 50
622  ! ... switch for fixed time stepping
623  LOGICAL, PUBLIC            :: l_fixed_step = .FALSE.
624  INTEGER, PUBLIC            :: nfsteps = 1
625  ! ... number of kpp control PARAMETERs
626  INTEGER, PARAMETER, PUBLIC :: nkppctrl = 20
627  !
628  INTEGER, DIMENSION(nkppctrl), PUBLIC     :: icntrl = 0
629  REAL(dp), DIMENSION(nkppctrl), PUBLIC     :: rcntrl = 0.0_dp
630  REAL(dp), DIMENSION(nmaxfixsteps), PUBLIC :: t_steps = 0.0_dp
631
632  ! END header MODULE initialize_kpp_ctrl_template
633
634 
635! Interface Block
636 
637  INTERFACE            initialize
638    MODULE PROCEDURE   initialize
639  END INTERFACE        initialize
640 
641  INTERFACE            integrate
642    MODULE PROCEDURE   integrate
643  END INTERFACE        integrate
644 
645  INTERFACE            fun
646    MODULE PROCEDURE   fun
647  END INTERFACE        fun
648 
649  INTERFACE            kppsolve
650    MODULE PROCEDURE   kppsolve
651  END INTERFACE        kppsolve
652 
653  INTERFACE            jac_sp
654    MODULE PROCEDURE   jac_sp
655  END INTERFACE        jac_sp
656 
657  INTERFACE            k_arr
658    MODULE PROCEDURE   k_arr
659  END INTERFACE        k_arr
660 
661  INTERFACE            update_rconst
662    MODULE PROCEDURE   update_rconst
663  END INTERFACE        update_rconst
664 
665  INTERFACE            arr2
666    MODULE PROCEDURE   arr2
667  END INTERFACE        arr2
668 
669  INTERFACE            initialize_kpp_ctrl
670    MODULE PROCEDURE   initialize_kpp_ctrl
671  END INTERFACE        initialize_kpp_ctrl
672 
673  INTERFACE            error_output
674    MODULE PROCEDURE   error_output
675  END INTERFACE        error_output
676 
677  INTERFACE            wscal
678    MODULE PROCEDURE   wscal
679  END INTERFACE        wscal
680 
681!INTERFACE not working  INTERFACE            waxpy
682!INTERFACE not working    MODULE PROCEDURE   waxpy
683!INTERFACE not working  END INTERFACE        waxpy
684 
685  INTERFACE            rosenbrock
686    MODULE PROCEDURE   rosenbrock
687  END INTERFACE        rosenbrock
688 
689  INTERFACE            funtemplate
690    MODULE PROCEDURE   funtemplate
691  END INTERFACE        funtemplate
692 
693  INTERFACE            jactemplate
694    MODULE PROCEDURE   jactemplate
695  END INTERFACE        jactemplate
696 
697  INTERFACE            kppdecomp
698    MODULE PROCEDURE   kppdecomp
699  END INTERFACE        kppdecomp
700 
701  INTERFACE            get_mechanismname
702    MODULE PROCEDURE   get_mechanismname
703  END INTERFACE        get_mechanismname
704 
705  INTERFACE            chem_gasphase_integrate
706    MODULE PROCEDURE   chem_gasphase_integrate
707  END INTERFACE        chem_gasphase_integrate
708 
709 
710  ! openmp directives generated by kp4
711 
712  !$OMP THREADPRIVATE (vl, vl_glo, is, ie, data_loaded)
713  !$OMP THREADPRIVATE (c, var, fix, rconst, time, temp, stepmin, cfactor)
714  !$OMP THREADPRIVATE (qvap, fakt, cs_mech, a, icntrl, rcntrl)
715 
716 CONTAINS
717 
718SUBROUTINE initialize()
719
720
721 INTEGER         :: k
722
723  INTEGER :: i
724  REAL(kind=dp):: x
725  k = is
726  cfactor = 1.000000e+00_dp
727
728! Following line is just to avoid compiler message about unused variables
729IF ( lu_crow(1) == 1  .OR.  lu_icol(1) == 1  .OR.  lu_irow(1) == 1 )  CONTINUE 
730
731
732  x = (0.) * cfactor
733  DO i = 1 , nvar
734  ENDDO
735
736  x = (0.) * cfactor
737  DO i = 1 , nfix
738    fix(i) = x
739  ENDDO
740
741! constant rate coefficients
742! END constant rate coefficients
743
744! INLINED initializations
745
746! End INLINED initializations
747
748     
749END SUBROUTINE initialize
750 
751SUBROUTINE integrate( tin, tout, &
752  icntrl_u, rcntrl_u, istatus_u, rstatus_u, ierr_u)
753
754
755   REAL(kind=dp), INTENT(IN):: tin  ! start time
756   REAL(kind=dp), INTENT(IN):: tout ! END time
757   ! OPTIONAL input PARAMETERs and statistics
758   INTEGER,      INTENT(IN), OPTIONAL :: icntrl_u(20)
759   REAL(kind=dp), INTENT(IN), OPTIONAL :: rcntrl_u(20)
760   INTEGER,      INTENT(OUT), OPTIONAL :: istatus_u(20)
761   REAL(kind=dp), INTENT(OUT), OPTIONAL :: rstatus_u(20)
762   INTEGER,      INTENT(OUT), OPTIONAL :: ierr_u
763
764   REAL(kind=dp):: rcntrl(20), rstatus(20)
765   INTEGER       :: icntrl(20), istatus(20), ierr
766
767
768   icntrl(:) = 0
769   rcntrl(:) = 0.0_dp
770   istatus(:) = 0
771   rstatus(:) = 0.0_dp
772
773    !~~~> fine-tune the integrator:
774   icntrl(1) = 0      ! 0 - non- autonomous, 1 - autonomous
775   icntrl(2) = 0      ! 0 - vector tolerances, 1 - scalars
776
777   ! IF OPTIONAL PARAMETERs are given, and IF they are >0,
778   ! THEN they overwrite default settings.
779   IF (PRESENT(icntrl_u))THEN
780     WHERE(icntrl_u(:)> 0)icntrl(:) = icntrl_u(:)
781   ENDIF
782   IF (PRESENT(rcntrl_u))THEN
783     WHERE(rcntrl_u(:)> 0)rcntrl(:) = rcntrl_u(:)
784   ENDIF
785
786
787   CALL rosenbrock(nvar, var, tin, tout,  &
788         atol, rtol,               &
789         rcntrl, icntrl, rstatus, istatus, ierr)
790
791   !~~~> debug option: show no of steps
792   ! ntotal = ntotal + istatus(nstp)
793   ! PRINT*,'NSTEPS=',ISTATUS(Nstp),' (',Ntotal,')','  O3=',VAR(ind_O3)
794
795   stepmin = rstatus(nhexit)
796   ! IF OPTIONAL PARAMETERs are given for output they
797   ! are updated with the RETURN information
798   IF (PRESENT(istatus_u))istatus_u(:) = istatus(:)
799   IF (PRESENT(rstatus_u))rstatus_u(:) = rstatus(:)
800   IF (PRESENT(ierr_u))  ierr_u       = ierr
801
802END SUBROUTINE integrate
803 
804SUBROUTINE fun(v, f, rct, vdot)
805
806! V - Concentrations of variable species (local)
807  REAL(kind=dp):: v(nvar)
808! F - Concentrations of fixed species (local)
809  REAL(kind=dp):: f(nfix)
810! RCT - Rate constants (local)
811  REAL(kind=dp):: rct(nreact)
812! Vdot - Time derivative of variable species concentrations
813  REAL(kind=dp):: vdot(nvar)
814
815
816! Following line is just to avoid compiler message about unused variables
817IF ( f(nfix) > 0.0_dp )  CONTINUE
818
819! Computation of equation rates
820  a(1) = rct(1) * v(10)
821  a(2) = rct(2) * v(8) * f(1)
822  a(3) = rct(3) * v(8) * v(11)
823  a(4) = rct(4) * v(7) * v(9)
824  a(5) = rct(5) * v(11) * v(13)
825  a(6) = rct(6) * v(11) * v(12)
826  a(7) = rct(7) * v(9) * v(10)
827
828! Aggregate function
829  vdot(1) = 0
830  vdot(2) = 0
831  vdot(3) = 0
832  vdot(4) = 0
833  vdot(5) = a(7)
834  vdot(6) = a(5)
835  vdot(7) = - a(4)
836  vdot(8) = a(1) - a(2) - a(3)
837  vdot(9) = 2* a(2) - a(4) + a(6) - a(7)
838  vdot(10) = - a(1) + a(3) + a(5) + a(6) - a(7)
839  vdot(11) = a(1) - a(3) - a(5) - a(6)
840  vdot(12) = a(5) - a(6)
841  vdot(13) = a(4) - a(5)
842     
843END SUBROUTINE fun
844 
845SUBROUTINE kppsolve(jvs, x)
846
847! JVS - sparse Jacobian of variables
848  REAL(kind=dp):: jvs(lu_nonzero)
849! X - Vector for variables
850  REAL(kind=dp):: x(nvar)
851
852  x(9) = x(9) - jvs(16) * x(7) - jvs(17) * x(8)
853  x(10) = x(10) - jvs(22) * x(8) - jvs(23) * x(9)
854  x(11) = x(11) - jvs(28) * x(8) - jvs(29) * x(10)
855  x(12) = x(12) - jvs(33) * x(11)
856  x(13) = x(13) - jvs(36) * x(7) - jvs(37) * x(9) - jvs(38) * x(10) - jvs(39) * x(11) - jvs(40) * x(12)
857  x(13) = x(13) / jvs(41)
858  x(12) = (x(12) - jvs(35) * x(13)) /(jvs(34))
859  x(11) = (x(11) - jvs(31) * x(12) - jvs(32) * x(13)) /(jvs(30))
860  x(10) = (x(10) - jvs(25) * x(11) - jvs(26) * x(12) - jvs(27) * x(13)) /(jvs(24))
861  x(9) = (x(9) - jvs(19) * x(10) - jvs(20) * x(11) - jvs(21) * x(12)) /(jvs(18))
862  x(8) = (x(8) - jvs(14) * x(10) - jvs(15) * x(11)) /(jvs(13))
863  x(7) = (x(7) - jvs(12) * x(9)) /(jvs(11))
864  x(6) = (x(6) - jvs(9) * x(11) - jvs(10) * x(13)) /(jvs(8))
865  x(5) = (x(5) - jvs(6) * x(9) - jvs(7) * x(10)) /(jvs(5))
866  x(4) = x(4) / jvs(4)
867  x(3) = x(3) / jvs(3)
868  x(2) = x(2) / jvs(2)
869  x(1) = x(1) / jvs(1)
870     
871END SUBROUTINE kppsolve
872 
873SUBROUTINE jac_sp(v, f, rct, jvs)
874
875! V - Concentrations of variable species (local)
876  REAL(kind=dp):: v(nvar)
877! F - Concentrations of fixed species (local)
878  REAL(kind=dp):: f(nfix)
879! RCT - Rate constants (local)
880  REAL(kind=dp):: rct(nreact)
881! JVS - sparse Jacobian of variables
882  REAL(kind=dp):: jvs(lu_nonzero)
883
884
885! Local variables
886! B - Temporary array
887  REAL(kind=dp):: b(17)
888!
889! Following line is just to avoid compiler message about unused variables
890IF ( f(nfix) > 0.0_dp )  CONTINUE
891
892! B(1) = dA(1)/dV(10)
893  b(1) = rct(1)
894! B(2) = dA(2)/dV(8)
895  b(2) = rct(2) * f(1)
896! B(4) = dA(3)/dV(8)
897  b(4) = rct(3) * v(11)
898! B(5) = dA(3)/dV(11)
899  b(5) = rct(3) * v(8)
900! B(6) = dA(4)/dV(7)
901  b(6) = rct(4) * v(9)
902! B(7) = dA(4)/dV(9)
903  b(7) = rct(4) * v(7)
904! B(8) = dA(5)/dV(11)
905  b(8) = rct(5) * v(13)
906! B(9) = dA(5)/dV(13)
907  b(9) = rct(5) * v(11)
908! B(10) = dA(6)/dV(11)
909  b(10) = rct(6) * v(12)
910! B(11) = dA(6)/dV(12)
911  b(11) = rct(6) * v(11)
912! B(12) = dA(7)/dV(9)
913  b(12) = rct(7) * v(10)
914! B(13) = dA(7)/dV(10)
915  b(13) = rct(7) * v(9)
916! B(14) = dA(8)/dV(1)
917  b(14) = rct(8)
918! B(15) = dA(9)/dV(2)
919  b(15) = rct(9)
920! B(16) = dA(10)/dV(3)
921  b(16) = rct(10)
922! B(17) = dA(11)/dV(4)
923  b(17) = rct(11)
924
925! Construct the Jacobian terms from B's
926! JVS(1) = Jac_FULL(1,1)
927  jvs(1) = 0
928! JVS(2) = Jac_FULL(2,2)
929  jvs(2) = 0
930! JVS(3) = Jac_FULL(3,3)
931  jvs(3) = 0
932! JVS(4) = Jac_FULL(4,4)
933  jvs(4) = 0
934! JVS(5) = Jac_FULL(5,5)
935  jvs(5) = 0
936! JVS(6) = Jac_FULL(5,9)
937  jvs(6) = b(12)
938! JVS(7) = Jac_FULL(5,10)
939  jvs(7) = b(13)
940! JVS(8) = Jac_FULL(6,6)
941  jvs(8) = 0
942! JVS(9) = Jac_FULL(6,11)
943  jvs(9) = b(8)
944! JVS(10) = Jac_FULL(6,13)
945  jvs(10) = b(9)
946! JVS(11) = Jac_FULL(7,7)
947  jvs(11) = - b(6)
948! JVS(12) = Jac_FULL(7,9)
949  jvs(12) = - b(7)
950! JVS(13) = Jac_FULL(8,8)
951  jvs(13) = - b(2) - b(4)
952! JVS(14) = Jac_FULL(8,10)
953  jvs(14) = b(1)
954! JVS(15) = Jac_FULL(8,11)
955  jvs(15) = - b(5)
956! JVS(16) = Jac_FULL(9,7)
957  jvs(16) = - b(6)
958! JVS(17) = Jac_FULL(9,8)
959  jvs(17) = 2* b(2)
960! JVS(18) = Jac_FULL(9,9)
961  jvs(18) = - b(7) - b(12)
962! JVS(19) = Jac_FULL(9,10)
963  jvs(19) = - b(13)
964! JVS(20) = Jac_FULL(9,11)
965  jvs(20) = b(10)
966! JVS(21) = Jac_FULL(9,12)
967  jvs(21) = b(11)
968! JVS(22) = Jac_FULL(10,8)
969  jvs(22) = b(4)
970! JVS(23) = Jac_FULL(10,9)
971  jvs(23) = - b(12)
972! JVS(24) = Jac_FULL(10,10)
973  jvs(24) = - b(1) - b(13)
974! JVS(25) = Jac_FULL(10,11)
975  jvs(25) = b(5) + b(8) + b(10)
976! JVS(26) = Jac_FULL(10,12)
977  jvs(26) = b(11)
978! JVS(27) = Jac_FULL(10,13)
979  jvs(27) = b(9)
980! JVS(28) = Jac_FULL(11,8)
981  jvs(28) = - b(4)
982! JVS(29) = Jac_FULL(11,10)
983  jvs(29) = b(1)
984! JVS(30) = Jac_FULL(11,11)
985  jvs(30) = - b(5) - b(8) - b(10)
986! JVS(31) = Jac_FULL(11,12)
987  jvs(31) = - b(11)
988! JVS(32) = Jac_FULL(11,13)
989  jvs(32) = - b(9)
990! JVS(33) = Jac_FULL(12,11)
991  jvs(33) = b(8) - b(10)
992! JVS(34) = Jac_FULL(12,12)
993  jvs(34) = - b(11)
994! JVS(35) = Jac_FULL(12,13)
995  jvs(35) = b(9)
996! JVS(36) = Jac_FULL(13,7)
997  jvs(36) = b(6)
998! JVS(37) = Jac_FULL(13,9)
999  jvs(37) = b(7)
1000! JVS(38) = Jac_FULL(13,10)
1001  jvs(38) = 0
1002! JVS(39) = Jac_FULL(13,11)
1003  jvs(39) = - b(8)
1004! JVS(40) = Jac_FULL(13,12)
1005  jvs(40) = 0
1006! JVS(41) = Jac_FULL(13,13)
1007  jvs(41) = - b(9)
1008     
1009END SUBROUTINE jac_sp
1010 
1011  elemental REAL(kind=dp)FUNCTION k_arr (k_298, tdep, temp)
1012    ! arrhenius FUNCTION
1013
1014    REAL,    INTENT(IN):: k_298 ! k at t = 298.15k
1015    REAL,    INTENT(IN):: tdep  ! temperature dependence
1016    REAL(kind=dp), INTENT(IN):: temp  ! temperature
1017
1018    intrinsic exp
1019
1020    k_arr = k_298 * exp(tdep* (1._dp/temp- 3.3540e-3_dp))! 1/298.15=3.3540e-3
1021
1022  END FUNCTION k_arr
1023 
1024SUBROUTINE update_rconst()
1025 INTEGER         :: k
1026
1027  k = is
1028
1029! Begin INLINED RCONST
1030
1031
1032! End INLINED RCONST
1033
1034  rconst(1) = (phot(j_no2))
1035  rconst(2) = (2.0_dp * 2.2e-10_dp * phot(j_o31d) / (arr2(1.9e+8_dp , -390.0_dp , temp)))
1036  rconst(3) = (arr2(1.80e-12_dp , 1370.0_dp , temp))
1037  rconst(4) = (arr2(2.00e-11_dp , 500.0_dp , temp))
1038  rconst(5) = (arr2(4.20e-12_dp , -180.0_dp , temp))
1039  rconst(6) = (arr2(3.70e-12_dp , -240.0_dp , temp))
1040  rconst(7) = (arr2(1.15e-11_dp , 0.0_dp , temp))
1041  rconst(8) = (1.0_dp)
1042  rconst(9) = (1.0_dp)
1043  rconst(10) = (1.0_dp)
1044  rconst(11) = (1.0_dp)
1045     
1046END SUBROUTINE update_rconst
1047 
1048!  END FUNCTION ARR2
1049REAL(kind=dp)FUNCTION arr2( a0, b0, temp)
1050    REAL(kind=dp):: temp
1051    REAL(kind=dp):: a0, b0
1052    arr2 = a0 * exp( - b0 / temp)
1053END FUNCTION arr2
1054 
1055SUBROUTINE initialize_kpp_ctrl(status)
1056
1057
1058  ! i/o
1059  INTEGER,         INTENT(OUT):: status
1060
1061  ! local
1062  REAL(dp):: tsum
1063  INTEGER  :: i
1064
1065  ! check fixed time steps
1066  tsum = 0.0_dp
1067  DO i=1, nmaxfixsteps
1068     IF (t_steps(i)< tiny(0.0_dp))exit
1069     tsum = tsum + t_steps(i)
1070  ENDDO
1071
1072  nfsteps = i- 1
1073
1074  l_fixed_step = (nfsteps > 0).and.((tsum - 1.0)< tiny(0.0_dp))
1075
1076  IF (l_vector)THEN
1077     WRITE(*,*) ' MODE             : VECTOR (LENGTH=',VL_DIM,')'
1078  ELSE
1079     WRITE(*,*) ' MODE             : SCALAR'
1080  ENDIF
1081  !
1082  WRITE(*,*) ' DE-INDEXING MODE :',I_LU_DI
1083  !
1084  WRITE(*,*) ' ICNTRL           : ',icntrl
1085  WRITE(*,*) ' RCNTRL           : ',rcntrl
1086  !
1087  ! note: this is ONLY meaningful for vectorized (kp4)rosenbrock- methods
1088  IF (l_vector)THEN
1089     IF (l_fixed_step)THEN
1090        WRITE(*,*) ' TIME STEPS       : FIXED (',t_steps(1:nfsteps),')'
1091     ELSE
1092        WRITE(*,*) ' TIME STEPS       : AUTOMATIC'
1093     ENDIF
1094  ELSE
1095     WRITE(*,*) ' TIME STEPS       : AUTOMATIC '//&
1096          &'(t_steps (CTRL_KPP) ignored in SCALAR MODE)'
1097  ENDIF
1098  ! mz_pj_20070531-
1099
1100  status = 0
1101
1102
1103END SUBROUTINE initialize_kpp_ctrl
1104 
1105SUBROUTINE error_output(c, ierr, pe)
1106
1107
1108  INTEGER, INTENT(IN):: ierr
1109  INTEGER, INTENT(IN):: pe
1110  REAL(dp), DIMENSION(:), INTENT(IN):: c
1111
1112  write(6,*) 'ERROR in chem_gasphase_mod ',ierr,C(1),PE
1113
1114
1115END SUBROUTINE error_output
1116 
1117      SUBROUTINE wscal(n, alpha, x, incx)
1118!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1119!     constant times a vector: x(1:N) <- Alpha*x(1:N)
1120!     only for incX=incY=1
1121!     after BLAS
1122!     replace this by the function from the optimized BLAS implementation:
1123!         CALL SSCAL(N,Alpha,X,1) or  CALL DSCAL(N,Alpha,X,1)
1124!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1125
1126      INTEGER  :: i, incx, m, mp1, n
1127      REAL(kind=dp) :: x(n), alpha
1128      REAL(kind=dp), PARAMETER  :: zero=0.0_dp, one=1.0_dp
1129
1130! Following line is just to avoid compiler message about unused variables
1131IF ( incx == 0 )  CONTINUE
1132
1133      IF (alpha .eq. one)RETURN
1134      IF (n .le. 0)RETURN
1135
1136      m = mod(n, 5)
1137      IF ( m .ne. 0)THEN
1138        IF (alpha .eq. (- one))THEN
1139          DO i = 1, m
1140            x(i) = - x(i)
1141          ENDDO
1142        ELSEIF (alpha .eq. zero)THEN
1143          DO i = 1, m
1144            x(i) = zero
1145          ENDDO
1146        ELSE
1147          DO i = 1, m
1148            x(i) = alpha* x(i)
1149          ENDDO
1150        ENDIF
1151        IF ( n .lt. 5)RETURN
1152      ENDIF
1153      mp1 = m + 1
1154      IF (alpha .eq. (- one))THEN
1155        DO i = mp1, n, 5
1156          x(i)   = - x(i)
1157          x(i + 1) = - x(i + 1)
1158          x(i + 2) = - x(i + 2)
1159          x(i + 3) = - x(i + 3)
1160          x(i + 4) = - x(i + 4)
1161        ENDDO
1162      ELSEIF (alpha .eq. zero)THEN
1163        DO i = mp1, n, 5
1164          x(i)   = zero
1165          x(i + 1) = zero
1166          x(i + 2) = zero
1167          x(i + 3) = zero
1168          x(i + 4) = zero
1169        ENDDO
1170      ELSE
1171        DO i = mp1, n, 5
1172          x(i)   = alpha* x(i)
1173          x(i + 1) = alpha* x(i + 1)
1174          x(i + 2) = alpha* x(i + 2)
1175          x(i + 3) = alpha* x(i + 3)
1176          x(i + 4) = alpha* x(i + 4)
1177        ENDDO
1178      ENDIF
1179
1180      END SUBROUTINE wscal
1181 
1182      SUBROUTINE waxpy(n, alpha, x, incx, y, incy)
1183!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1184!     constant times a vector plus a vector: y <- y + Alpha*x
1185!     only for incX=incY=1
1186!     after BLAS
1187!     replace this by the function from the optimized BLAS implementation:
1188!         CALL SAXPY(N,Alpha,X,1,Y,1) or  CALL DAXPY(N,Alpha,X,1,Y,1)
1189!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1190
1191      INTEGER  :: i, incx, incy, m, mp1, n
1192      REAL(kind=dp):: x(n), y(n), alpha
1193      REAL(kind=dp), PARAMETER :: zero = 0.0_dp
1194
1195
1196! Following line is just to avoid compiler message about unused variables
1197IF ( incx == 0  .OR.  incy == 0 )  CONTINUE
1198      IF (alpha .eq. zero)RETURN
1199      IF (n .le. 0)RETURN
1200
1201      m = mod(n, 4)
1202      IF ( m .ne. 0)THEN
1203        DO i = 1, m
1204          y(i) = y(i) + alpha* x(i)
1205        ENDDO
1206        IF ( n .lt. 4)RETURN
1207      ENDIF
1208      mp1 = m + 1
1209      DO i = mp1, n, 4
1210        y(i) = y(i) + alpha* x(i)
1211        y(i + 1) = y(i + 1) + alpha* x(i + 1)
1212        y(i + 2) = y(i + 2) + alpha* x(i + 2)
1213        y(i + 3) = y(i + 3) + alpha* x(i + 3)
1214      ENDDO
1215     
1216      END SUBROUTINE waxpy
1217 
1218SUBROUTINE rosenbrock(n, y, tstart, tend, &
1219           abstol, reltol,             &
1220           rcntrl, icntrl, rstatus, istatus, ierr)
1221!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1222!
1223!    Solves the system y'=F(t,y) using a Rosenbrock method defined by:
1224!
1225!     G = 1/(H*gamma(1)) - Jac(t0,Y0)
1226!     T_i = t0 + Alpha(i)*H
1227!     Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
1228!     G *K_i = Fun( T_i,Y_i)+ \sum_{j=1}^S C(i,j)/H *K_j +
1229!         gamma(i)*dF/dT(t0,Y0)
1230!     Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
1231!
1232!    For details on Rosenbrock methods and their implementation consult:
1233!      E. Hairer and G. Wanner
1234!      "Solving ODEs II. Stiff and differential-algebraic problems".
1235!      Springer series in computational mathematics,Springer-Verlag,1996.
1236!    The codes contained in the book inspired this implementation.
1237!
1238!    (C)  Adrian Sandu,August 2004
1239!    Virginia Polytechnic Institute and State University
1240!    Contact: sandu@cs.vt.edu
1241!    Revised by Philipp Miehe and Adrian Sandu,May 2006                 
1242!    This implementation is part of KPP - the Kinetic PreProcessor
1243!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1244!
1245!~~~>   input arguments:
1246!
1247!-     y(n)  = vector of initial conditions (at t=tstart)
1248!-    [tstart, tend]  = time range of integration
1249!     (if Tstart>Tend the integration is performed backwards in time)
1250!-    reltol, abstol = user precribed accuracy
1251!- SUBROUTINE  fun( t, y, ydot) = ode FUNCTION,
1252!                       returns Ydot = Y' = F(T,Y)
1253!- SUBROUTINE  jac( t, y, jcb) = jacobian of the ode FUNCTION,
1254!                       returns Jcb = dFun/dY
1255!-    icntrl(1:20)  = INTEGER inputs PARAMETERs
1256!-    rcntrl(1:20)  = REAL inputs PARAMETERs
1257!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1258!
1259!~~~>     output arguments:
1260!
1261!-    y(n)  - > vector of final states (at t- >tend)
1262!-    istatus(1:20) - > INTEGER output PARAMETERs
1263!-    rstatus(1:20) - > REAL output PARAMETERs
1264!-    ierr            - > job status upon RETURN
1265!                        success (positive value) or
1266!                        failure (negative value)
1267!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1268!
1269!~~~>     input PARAMETERs:
1270!
1271!    Note: For input parameters equal to zero the default values of the
1272!       corresponding variables are used.
1273!
1274!    ICNTRL(1) = 1: F = F(y)   Independent of T (AUTONOMOUS)
1275!              = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
1276!
1277!    ICNTRL(2) = 0: AbsTol,RelTol are N-dimensional vectors
1278!              = 1: AbsTol,RelTol are scalars
1279!
1280!    ICNTRL(3)  -> selection of a particular Rosenbrock method
1281!        = 0 :    Rodas3 (default)
1282!        = 1 :    Ros2
1283!        = 2 :    Ros3
1284!        = 3 :    Ros4
1285!        = 4 :    Rodas3
1286!        = 5 :    Rodas4
1287!
1288!    ICNTRL(4)  -> maximum number of integration steps
1289!        For ICNTRL(4) =0) the default value of 100000 is used
1290!
1291!    RCNTRL(1)  -> Hmin,lower bound for the integration step size
1292!          It is strongly recommended to keep Hmin = ZERO
1293!    RCNTRL(2)  -> Hmax,upper bound for the integration step size
1294!    RCNTRL(3)  -> Hstart,starting value for the integration step size
1295!
1296!    RCNTRL(4)  -> FacMin,lower bound on step decrease factor (default=0.2)
1297!    RCNTRL(5)  -> FacMax,upper bound on step increase factor (default=6)
1298!    RCNTRL(6)  -> FacRej,step decrease factor after multiple rejections
1299!                          (default=0.1)
1300!    RCNTRL(7)  -> FacSafe,by which the new step is slightly smaller
1301!         than the predicted value  (default=0.9)
1302!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1303!
1304!
1305!    OUTPUT ARGUMENTS:
1306!    -----------------
1307!
1308!    T           -> T value for which the solution has been computed
1309!                     (after successful return T=Tend).
1310!
1311!    Y(N)        -> Numerical solution at T
1312!
1313!    IDID        -> Reports on successfulness upon return:
1314!                    = 1 for success
1315!                    < 0 for error (value equals error code)
1316!
1317!    ISTATUS(1)  -> No. of function calls
1318!    ISTATUS(2)  -> No. of jacobian calls
1319!    ISTATUS(3)  -> No. of steps
1320!    ISTATUS(4)  -> No. of accepted steps
1321!    ISTATUS(5)  -> No. of rejected steps (except at very beginning)
1322!    ISTATUS(6)  -> No. of LU decompositions
1323!    ISTATUS(7)  -> No. of forward/backward substitutions
1324!    ISTATUS(8)  -> No. of singular matrix decompositions
1325!
1326!    RSTATUS(1)  -> Texit,the time corresponding to the
1327!                     computed Y upon return
1328!    RSTATUS(2)  -> Hexit,last accepted step before exit
1329!    RSTATUS(3)  -> Hnew,last predicted step (not yet taken)
1330!                   For multiple restarts,use Hnew as Hstart
1331!                     in the subsequent run
1332!
1333!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1334
1335
1336!~~~>  arguments
1337   INTEGER,      INTENT(IN)  :: n
1338   REAL(kind=dp), INTENT(INOUT):: y(n)
1339   REAL(kind=dp), INTENT(IN)  :: tstart, tend
1340   REAL(kind=dp), INTENT(IN)  :: abstol(n), reltol(n)
1341   INTEGER,      INTENT(IN)  :: icntrl(20)
1342   REAL(kind=dp), INTENT(IN)  :: rcntrl(20)
1343   INTEGER,      INTENT(INOUT):: istatus(20)
1344   REAL(kind=dp), INTENT(INOUT):: rstatus(20)
1345   INTEGER, INTENT(OUT) :: ierr
1346!~~~>  PARAMETERs of the rosenbrock method, up to 6 stages
1347   INTEGER ::  ros_s, rosmethod
1348   INTEGER, PARAMETER :: rs2=1, rs3=2, rs4=3, rd3=4, rd4=5, rg3=6
1349   REAL(kind=dp):: ros_a(15), ros_c(15), ros_m(6), ros_e(6), &
1350                    ros_alpha(6), ros_gamma(6), ros_elo
1351   LOGICAL :: ros_newf(6)
1352   CHARACTER(len=12):: ros_name
1353!~~~>  local variables
1354   REAL(kind=dp):: roundoff, facmin, facmax, facrej, facsafe
1355   REAL(kind=dp):: hmin, hmax, hstart
1356   REAL(kind=dp):: texit
1357   INTEGER       :: i, uplimtol, max_no_steps
1358   LOGICAL       :: autonomous, vectortol
1359!~~~>   PARAMETERs
1360   REAL(kind=dp), PARAMETER :: zero = 0.0_dp, one  = 1.0_dp
1361   REAL(kind=dp), PARAMETER :: deltamin = 1.0e-5_dp
1362
1363!~~~>  initialize statistics
1364   istatus(1:8) = 0
1365   rstatus(1:3) = zero
1366
1367!~~~>  autonomous or time dependent ode. default is time dependent.
1368   autonomous = .not.(icntrl(1) == 0)
1369
1370!~~~>  for scalar tolerances (icntrl(2).ne.0) the code uses abstol(1)and reltol(1)
1371!   For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:N) and RelTol(1:N)
1372   IF (icntrl(2) == 0)THEN
1373      vectortol = .TRUE.
1374      uplimtol  = n
1375   ELSE
1376      vectortol = .FALSE.
1377      uplimtol  = 1
1378   ENDIF
1379
1380!~~~>   initialize the particular rosenbrock method selected
1381   select CASE (icntrl(3))
1382     CASE (1)
1383       CALL ros2
1384     CASE (2)
1385       CALL ros3
1386     CASE (3)
1387       CALL ros4
1388     CASE (0, 4)
1389       CALL rodas3
1390     CASE (5)
1391       CALL rodas4
1392     CASE (6)
1393       CALL rang3
1394     CASE default
1395       PRINT *,'Unknown Rosenbrock method: ICNTRL(3) =',ICNTRL(3) 
1396       CALL ros_errormsg(- 2, tstart, zero, ierr)
1397       RETURN
1398   END select
1399
1400!~~~>   the maximum number of steps admitted
1401   IF (icntrl(4) == 0)THEN
1402      max_no_steps = 200000
1403   ELSEIF (icntrl(4)> 0)THEN
1404      max_no_steps=icntrl(4)
1405   ELSE
1406      PRINT *,'User-selected max no. of steps: ICNTRL(4) =',ICNTRL(4)
1407      CALL ros_errormsg(- 1, tstart, zero, ierr)
1408      RETURN
1409   ENDIF
1410
1411!~~~>  unit roundoff (1+ roundoff>1)
1412   roundoff = epsilon(one)
1413
1414!~~~>  lower bound on the step size: (positive value)
1415   IF (rcntrl(1) == zero)THEN
1416      hmin = zero
1417   ELSEIF (rcntrl(1)> zero)THEN
1418      hmin = rcntrl(1)
1419   ELSE
1420      PRINT *,'User-selected Hmin: RCNTRL(1) =',RCNTRL(1)
1421      CALL ros_errormsg(- 3, tstart, zero, ierr)
1422      RETURN
1423   ENDIF
1424!~~~>  upper bound on the step size: (positive value)
1425   IF (rcntrl(2) == zero)THEN
1426      hmax = abs(tend-tstart)
1427   ELSEIF (rcntrl(2)> zero)THEN
1428      hmax = min(abs(rcntrl(2)), abs(tend-tstart))
1429   ELSE
1430      PRINT *,'User-selected Hmax: RCNTRL(2) =',RCNTRL(2)
1431      CALL ros_errormsg(- 3, tstart, zero, ierr)
1432      RETURN
1433   ENDIF
1434!~~~>  starting step size: (positive value)
1435   IF (rcntrl(3) == zero)THEN
1436      hstart = max(hmin, deltamin)
1437   ELSEIF (rcntrl(3)> zero)THEN
1438      hstart = min(abs(rcntrl(3)), abs(tend-tstart))
1439   ELSE
1440      PRINT *,'User-selected Hstart: RCNTRL(3) =',RCNTRL(3)
1441      CALL ros_errormsg(- 3, tstart, zero, ierr)
1442      RETURN
1443   ENDIF
1444!~~~>  step size can be changed s.t.  facmin < hnew/hold < facmax
1445   IF (rcntrl(4) == zero)THEN
1446      facmin = 0.2_dp
1447   ELSEIF (rcntrl(4)> zero)THEN
1448      facmin = rcntrl(4)
1449   ELSE
1450      PRINT *,'User-selected FacMin: RCNTRL(4) =',RCNTRL(4)
1451      CALL ros_errormsg(- 4, tstart, zero, ierr)
1452      RETURN
1453   ENDIF
1454   IF (rcntrl(5) == zero)THEN
1455      facmax = 6.0_dp
1456   ELSEIF (rcntrl(5)> zero)THEN
1457      facmax = rcntrl(5)
1458   ELSE
1459      PRINT *,'User-selected FacMax: RCNTRL(5) =',RCNTRL(5)
1460      CALL ros_errormsg(- 4, tstart, zero, ierr)
1461      RETURN
1462   ENDIF
1463!~~~>   facrej: factor to decrease step after 2 succesive rejections
1464   IF (rcntrl(6) == zero)THEN
1465      facrej = 0.1_dp
1466   ELSEIF (rcntrl(6)> zero)THEN
1467      facrej = rcntrl(6)
1468   ELSE
1469      PRINT *,'User-selected FacRej: RCNTRL(6) =',RCNTRL(6)
1470      CALL ros_errormsg(- 4, tstart, zero, ierr)
1471      RETURN
1472   ENDIF
1473!~~~>   facsafe: safety factor in the computation of new step size
1474   IF (rcntrl(7) == zero)THEN
1475      facsafe = 0.9_dp
1476   ELSEIF (rcntrl(7)> zero)THEN
1477      facsafe = rcntrl(7)
1478   ELSE
1479      PRINT *,'User-selected FacSafe: RCNTRL(7) =',RCNTRL(7)
1480      CALL ros_errormsg(- 4, tstart, zero, ierr)
1481      RETURN
1482   ENDIF
1483!~~~>  check IF tolerances are reasonable
1484    DO i=1, uplimtol
1485      IF ((abstol(i)<= zero).or. (reltol(i)<= 10.0_dp* roundoff)&
1486         .or. (reltol(i)>= 1.0_dp))THEN
1487        PRINT *,' AbsTol(',i,') = ',AbsTol(i)
1488        PRINT *,' RelTol(',i,') = ',RelTol(i)
1489        CALL ros_errormsg(- 5, tstart, zero, ierr)
1490        RETURN
1491      ENDIF
1492    ENDDO
1493
1494
1495!~~~>  CALL rosenbrock method
1496   CALL ros_integrator(y, tstart, tend, texit,  &
1497        abstol, reltol,                         &
1498!  Integration parameters
1499        autonomous, vectortol, max_no_steps,    &
1500        roundoff, hmin, hmax, hstart,           &
1501        facmin, facmax, facrej, facsafe,        &
1502!  Error indicator
1503        ierr)
1504
1505!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1506CONTAINS !  SUBROUTINEs internal to rosenbrock
1507!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1508   
1509!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
1510 SUBROUTINE ros_errormsg(code, t, h, ierr)
1511!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
1512!    Handles all error messages
1513!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
1514 
1515   REAL(kind=dp), INTENT(IN):: t, h
1516   INTEGER, INTENT(IN) :: code
1517   INTEGER, INTENT(OUT):: ierr
1518   
1519   ierr = code
1520   print * , &
1521     'Forced exit from Rosenbrock due to the following error:' 
1522     
1523   select CASE (code)
1524    CASE (- 1) 
1525      PRINT *,'--> Improper value for maximal no of steps'
1526    CASE (- 2) 
1527      PRINT *,'--> Selected Rosenbrock method not implemented'
1528    CASE (- 3) 
1529      PRINT *,'--> Hmin/Hmax/Hstart must be positive'
1530    CASE (- 4) 
1531      PRINT *,'--> FacMin/FacMax/FacRej must be positive'
1532    CASE (- 5)
1533      PRINT *,'--> Improper tolerance values'
1534    CASE (- 6)
1535      PRINT *,'--> No of steps exceeds maximum bound'
1536    CASE (- 7)
1537      PRINT *,'--> Step size too small: T + 10*H = T',&
1538            ' or H < Roundoff'
1539    CASE (- 8) 
1540      PRINT *,'--> Matrix is repeatedly singular'
1541    CASE default
1542      PRINT *,'Unknown Error code: ',Code
1543   END select
1544   
1545   print * , "t=", t, "and h=", h
1546     
1547 END SUBROUTINE ros_errormsg
1548
1549!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1550 SUBROUTINE ros_integrator (y, tstart, tend, t, &
1551        abstol, reltol,                         &
1552!~~~> integration PARAMETERs
1553        autonomous, vectortol, max_no_steps,    &
1554        roundoff, hmin, hmax, hstart,           &
1555        facmin, facmax, facrej, facsafe,        &
1556!~~~> error indicator
1557        ierr)
1558!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1559!   Template for the implementation of a generic Rosenbrock method
1560!      defined by ros_S (no of stages)
1561!      and its coefficients ros_{A,C,M,E,Alpha,Gamma}
1562!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1563
1564
1565!~~~> input: the initial condition at tstart; output: the solution at t
1566   REAL(kind=dp), INTENT(INOUT):: y(n)
1567!~~~> input: integration interval
1568   REAL(kind=dp), INTENT(IN):: tstart, tend
1569!~~~> output: time at which the solution is RETURNed (t=tendIF success)
1570   REAL(kind=dp), INTENT(OUT)::  t
1571!~~~> input: tolerances
1572   REAL(kind=dp), INTENT(IN)::  abstol(n), reltol(n)
1573!~~~> input: integration PARAMETERs
1574   LOGICAL, INTENT(IN):: autonomous, vectortol
1575   REAL(kind=dp), INTENT(IN):: hstart, hmin, hmax
1576   INTEGER, INTENT(IN):: max_no_steps
1577   REAL(kind=dp), INTENT(IN):: roundoff, facmin, facmax, facrej, facsafe
1578!~~~> output: error indicator
1579   INTEGER, INTENT(OUT):: ierr
1580! ~~~~ Local variables
1581   REAL(kind=dp):: ynew(n), fcn0(n), fcn(n)
1582   REAL(kind=dp):: k(n* ros_s), dfdt(n)
1583#ifdef full_algebra   
1584   REAL(kind=dp):: jac0(n, n), ghimj(n, n)
1585#else
1586   REAL(kind=dp):: jac0(lu_nonzero), ghimj(lu_nonzero)
1587#endif
1588   REAL(kind=dp):: h, hnew, hc, hg, fac, tau
1589   REAL(kind=dp):: err, yerr(n)
1590   INTEGER :: pivot(n), direction, ioffset, j, istage
1591   LOGICAL :: rejectlasth, rejectmoreh, singular
1592!~~~>  local PARAMETERs
1593   REAL(kind=dp), PARAMETER :: zero = 0.0_dp, one  = 1.0_dp
1594   REAL(kind=dp), PARAMETER :: deltamin = 1.0e-5_dp
1595!~~~>  locally called FUNCTIONs
1596!    REAL(kind=dp) WLAMCH
1597!    EXTERNAL WLAMCH
1598!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1599
1600
1601!~~~>  initial preparations
1602   t = tstart
1603   rstatus(nhexit) = zero
1604   h = min( max(abs(hmin), abs(hstart)), abs(hmax))
1605   IF (abs(h)<= 10.0_dp* roundoff)h = deltamin
1606
1607   IF (tend  >=  tstart)THEN
1608     direction = + 1
1609   ELSE
1610     direction = - 1
1611   ENDIF
1612   h = direction* h
1613
1614   rejectlasth=.FALSE.
1615   rejectmoreh=.FALSE.
1616
1617!~~~> time loop begins below
1618
1619timeloop: DO WHILE((direction > 0).and.((t- tend) + roundoff <= zero)&
1620       .or. (direction < 0).and.((tend-t) + roundoff <= zero))
1621
1622   IF (istatus(nstp)> max_no_steps)THEN  ! too many steps
1623      CALL ros_errormsg(- 6, t, h, ierr)
1624      RETURN
1625   ENDIF
1626   IF (((t+ 0.1_dp* h) == t).or.(h <= roundoff))THEN  ! step size too small
1627      CALL ros_errormsg(- 7, t, h, ierr)
1628      RETURN
1629   ENDIF
1630
1631!~~~>  limit h IF necessary to avoid going beyond tend
1632   h = min(h, abs(tend-t))
1633
1634!~~~>   compute the FUNCTION at current time
1635   CALL funtemplate(t, y, fcn0)
1636   istatus(nfun) = istatus(nfun) + 1
1637
1638!~~~>  compute the FUNCTION derivative with respect to t
1639   IF (.not.autonomous)THEN
1640      CALL ros_funtimederivative(t, roundoff, y, &
1641                fcn0, dfdt)
1642   ENDIF
1643
1644!~~~>   compute the jacobian at current time
1645   CALL jactemplate(t, y, jac0)
1646   istatus(njac) = istatus(njac) + 1
1647
1648!~~~>  repeat step calculation until current step accepted
1649untilaccepted: do
1650
1651   CALL ros_preparematrix(h, direction, ros_gamma(1), &
1652          jac0, ghimj, pivot, singular)
1653   IF (singular)THEN ! more than 5 consecutive failed decompositions
1654       CALL ros_errormsg(- 8, t, h, ierr)
1655       RETURN
1656   ENDIF
1657
1658!~~~>   compute the stages
1659stage: DO istage = 1, ros_s
1660
1661      ! current istage offset. current istage vector is k(ioffset+ 1:ioffset+ n)
1662       ioffset = n* (istage-1)
1663
1664      ! for the 1st istage the FUNCTION has been computed previously
1665       IF (istage == 1)THEN
1666         !slim: CALL wcopy(n, fcn0, 1, fcn, 1)
1667       fcn(1:n) = fcn0(1:n)
1668      ! istage>1 and a new FUNCTION evaluation is needed at the current istage
1669       ELSEIF(ros_newf(istage))THEN
1670         !slim: CALL wcopy(n, y, 1, ynew, 1)
1671       ynew(1:n) = y(1:n)
1672         DO j = 1, istage-1
1673           CALL waxpy(n, ros_a((istage-1) * (istage-2) /2+ j), &
1674            k(n* (j- 1) + 1), 1, ynew, 1)
1675         ENDDO
1676         tau = t + ros_alpha(istage) * direction* h
1677         CALL funtemplate(tau, ynew, fcn)
1678         istatus(nfun) = istatus(nfun) + 1
1679       ENDIF ! IF istage == 1 ELSEIF ros_newf(istage)
1680       !slim: CALL wcopy(n, fcn, 1, k(ioffset+ 1), 1)
1681       k(ioffset+ 1:ioffset+ n) = fcn(1:n)
1682       DO j = 1, istage-1
1683         hc = ros_c((istage-1) * (istage-2) /2+ j) /(direction* h)
1684         CALL waxpy(n, hc, k(n* (j- 1) + 1), 1, k(ioffset+ 1), 1)
1685       ENDDO
1686       IF ((.not. autonomous).and.(ros_gamma(istage).ne.zero))THEN
1687         hg = direction* h* ros_gamma(istage)
1688         CALL waxpy(n, hg, dfdt, 1, k(ioffset+ 1), 1)
1689       ENDIF
1690       CALL ros_solve(ghimj, pivot, k(ioffset+ 1))
1691
1692   END DO stage
1693
1694
1695!~~~>  compute the new solution
1696   !slim: CALL wcopy(n, y, 1, ynew, 1)
1697   ynew(1:n) = y(1:n)
1698   DO j=1, ros_s
1699         CALL waxpy(n, ros_m(j), k(n* (j- 1) + 1), 1, ynew, 1)
1700   ENDDO
1701
1702!~~~>  compute the error estimation
1703   !slim: CALL wscal(n, zero, yerr, 1)
1704   yerr(1:n) = zero
1705   DO j=1, ros_s
1706        CALL waxpy(n, ros_e(j), k(n* (j- 1) + 1), 1, yerr, 1)
1707   ENDDO
1708   err = ros_errornorm(y, ynew, yerr, abstol, reltol, vectortol)
1709
1710!~~~> new step size is bounded by facmin <= hnew/h <= facmax
1711   fac  = min(facmax, max(facmin, facsafe/err** (one/ros_elo)))
1712   hnew = h* fac
1713
1714!~~~>  check the error magnitude and adjust step size
1715   istatus(nstp) = istatus(nstp) + 1
1716   IF ((err <= one).or.(h <= hmin))THEN  !~~~> accept step
1717      istatus(nacc) = istatus(nacc) + 1
1718      !slim: CALL wcopy(n, ynew, 1, y, 1)
1719      y(1:n) = ynew(1:n)
1720      t = t + direction* h
1721      hnew = max(hmin, min(hnew, hmax))
1722      IF (rejectlasth)THEN  ! no step size increase after a rejected step
1723         hnew = min(hnew, h)
1724      ENDIF
1725      rstatus(nhexit) = h
1726      rstatus(nhnew) = hnew
1727      rstatus(ntexit) = t
1728      rejectlasth = .FALSE.
1729      rejectmoreh = .FALSE.
1730      h = hnew
1731      exit untilaccepted ! exit the loop: WHILE step not accepted
1732   ELSE           !~~~> reject step
1733      IF (rejectmoreh)THEN
1734         hnew = h* facrej
1735      ENDIF
1736      rejectmoreh = rejectlasth
1737      rejectlasth = .TRUE.
1738      h = hnew
1739      IF (istatus(nacc)>= 1) istatus(nrej) = istatus(nrej) + 1
1740   ENDIF ! err <= 1
1741
1742   END DO untilaccepted
1743
1744   END DO timeloop
1745
1746!~~~> succesful exit
1747   ierr = 1  !~~~> the integration was successful
1748
1749  END SUBROUTINE ros_integrator
1750
1751
1752!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1753  REAL(kind=dp)FUNCTION ros_errornorm(y, ynew, yerr, &
1754                               abstol, reltol, vectortol)
1755!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1756!~~~> computes the "scaled norm" of the error vector yerr
1757!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1758
1759! Input arguments
1760   REAL(kind=dp), INTENT(IN):: y(n), ynew(n), &
1761          yerr(n), abstol(n), reltol(n)
1762   LOGICAL, INTENT(IN)::  vectortol
1763! Local variables
1764   REAL(kind=dp):: err, scale, ymax
1765   INTEGER  :: i
1766   REAL(kind=dp), PARAMETER :: zero = 0.0_dp
1767
1768   err = zero
1769   DO i=1, n
1770     ymax = max(abs(y(i)), abs(ynew(i)))
1771     IF (vectortol)THEN
1772       scale = abstol(i) + reltol(i) * ymax
1773     ELSE
1774       scale = abstol(1) + reltol(1) * ymax
1775     ENDIF
1776     err = err+ (yerr(i) /scale) ** 2
1777   ENDDO
1778   err  = sqrt(err/n)
1779
1780   ros_errornorm = max(err, 1.0d-10)
1781
1782  END FUNCTION ros_errornorm
1783
1784
1785!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1786  SUBROUTINE ros_funtimederivative(t, roundoff, y, &
1787                fcn0, dfdt)
1788!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1789!~~~> the time partial derivative of the FUNCTION by finite differences
1790!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1791
1792!~~~> input arguments
1793   REAL(kind=dp), INTENT(IN):: t, roundoff, y(n), fcn0(n)
1794!~~~> output arguments
1795   REAL(kind=dp), INTENT(OUT):: dfdt(n)
1796!~~~> local variables
1797   REAL(kind=dp):: delta
1798   REAL(kind=dp), PARAMETER :: one = 1.0_dp, deltamin = 1.0e-6_dp
1799
1800   delta = sqrt(roundoff) * max(deltamin, abs(t))
1801   CALL funtemplate(t+ delta, y, dfdt)
1802   istatus(nfun) = istatus(nfun) + 1
1803   CALL waxpy(n, (- one), fcn0, 1, dfdt, 1)
1804   CALL wscal(n, (one/delta), dfdt, 1)
1805
1806  END SUBROUTINE ros_funtimederivative
1807
1808
1809!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1810  SUBROUTINE ros_preparematrix(h, direction, gam, &
1811             jac0, ghimj, pivot, singular)
1812! --- --- --- --- --- --- --- --- --- --- --- --- ---
1813!  Prepares the LHS matrix for stage calculations
1814!  1.  Construct Ghimj = 1/(H*ham) - Jac0
1815!      "(Gamma H) Inverse Minus Jacobian"
1816!  2.  Repeat LU decomposition of Ghimj until successful.
1817!       -half the step size if LU decomposition fails and retry
1818!       -exit after 5 consecutive fails
1819! --- --- --- --- --- --- --- --- --- --- --- --- ---
1820
1821!~~~> input arguments
1822#ifdef full_algebra   
1823   REAL(kind=dp), INTENT(IN)::  jac0(n, n)
1824#else
1825   REAL(kind=dp), INTENT(IN)::  jac0(lu_nonzero)
1826#endif   
1827   REAL(kind=dp), INTENT(IN)::  gam
1828   INTEGER, INTENT(IN)::  direction
1829!~~~> output arguments
1830#ifdef full_algebra   
1831   REAL(kind=dp), INTENT(OUT):: ghimj(n, n)
1832#else
1833   REAL(kind=dp), INTENT(OUT):: ghimj(lu_nonzero)
1834#endif   
1835   LOGICAL, INTENT(OUT)::  singular
1836   INTEGER, INTENT(OUT)::  pivot(n)
1837!~~~> inout arguments
1838   REAL(kind=dp), INTENT(INOUT):: h   ! step size is decreased when lu fails
1839!~~~> local variables
1840   INTEGER  :: i, ising, nconsecutive
1841   REAL(kind=dp):: ghinv
1842   REAL(kind=dp), PARAMETER :: one  = 1.0_dp, half = 0.5_dp
1843
1844   nconsecutive = 0
1845   singular = .TRUE.
1846
1847   DO WHILE (singular)
1848
1849!~~~>    construct ghimj = 1/(h* gam) - jac0
1850#ifdef full_algebra   
1851     !slim: CALL wcopy(n* n, jac0, 1, ghimj, 1)
1852     !slim: CALL wscal(n* n, (- one), ghimj, 1)
1853     ghimj = - jac0
1854     ghinv = one/(direction* h* gam)
1855     DO i=1, n
1856       ghimj(i, i) = ghimj(i, i) + ghinv
1857     ENDDO
1858#else
1859     !slim: CALL wcopy(lu_nonzero, jac0, 1, ghimj, 1)
1860     !slim: CALL wscal(lu_nonzero, (- one), ghimj, 1)
1861     ghimj(1:lu_nonzero) = - jac0(1:lu_nonzero)
1862     ghinv = one/(direction* h* gam)
1863     DO i=1, n
1864       ghimj(lu_diag(i)) = ghimj(lu_diag(i)) + ghinv
1865     ENDDO
1866#endif   
1867!~~~>    compute lu decomposition
1868     CALL ros_decomp( ghimj, pivot, ising)
1869     IF (ising == 0)THEN
1870!~~~>    IF successful done
1871        singular = .FALSE.
1872     ELSE ! ising .ne. 0
1873!~~~>    IF unsuccessful half the step size; IF 5 consecutive fails THEN RETURN
1874        istatus(nsng) = istatus(nsng) + 1
1875        nconsecutive = nconsecutive+1
1876        singular = .TRUE.
1877        PRINT*,'Warning: LU Decomposition returned ISING = ',ISING
1878        IF (nconsecutive <= 5)THEN ! less than 5 consecutive failed decompositions
1879           h = h* half
1880        ELSE  ! more than 5 consecutive failed decompositions
1881           RETURN
1882        ENDIF  ! nconsecutive
1883      ENDIF    ! ising
1884
1885   END DO ! WHILE singular
1886
1887  END SUBROUTINE ros_preparematrix
1888
1889
1890!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1891  SUBROUTINE ros_decomp( a, pivot, ising)
1892!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1893!  Template for the LU decomposition
1894!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1895!~~~> inout variables
1896#ifdef full_algebra   
1897   REAL(kind=dp), INTENT(INOUT):: a(n, n)
1898#else   
1899   REAL(kind=dp), INTENT(INOUT):: a(lu_nonzero)
1900#endif
1901!~~~> output variables
1902   INTEGER, INTENT(OUT):: pivot(n), ising
1903
1904#ifdef full_algebra   
1905   CALL  dgetrf( n, n, a, n, pivot, ising)
1906#else   
1907   CALL kppdecomp(a, ising)
1908   pivot(1) = 1
1909#endif
1910   istatus(ndec) = istatus(ndec) + 1
1911
1912  END SUBROUTINE ros_decomp
1913
1914
1915!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1916  SUBROUTINE ros_solve( a, pivot, b)
1917!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1918!  Template for the forward/backward substitution (using pre-computed LU decomposition)
1919!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1920!~~~> input variables
1921#ifdef full_algebra   
1922   REAL(kind=dp), INTENT(IN):: a(n, n)
1923   INTEGER :: ising
1924#else   
1925   REAL(kind=dp), INTENT(IN):: a(lu_nonzero)
1926#endif
1927   INTEGER, INTENT(IN):: pivot(n)
1928!~~~> inout variables
1929   REAL(kind=dp), INTENT(INOUT):: b(n)
1930
1931! Following line is just to avoid compiler message about unused variables
1932IF ( pivot(1) == 0 )  CONTINUE
1933
1934#ifdef full_algebra   
1935   CALL  DGETRS( 'N',N ,1,A,N,Pivot,b,N,ISING)
1936   IF (info < 0)THEN
1937      print* , "error in dgetrs. ising=", ising
1938   ENDIF 
1939#else   
1940   CALL kppsolve( a, b)
1941#endif
1942
1943   istatus(nsol) = istatus(nsol) + 1
1944
1945  END SUBROUTINE ros_solve
1946
1947
1948
1949!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1950  SUBROUTINE ros2
1951!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1952! --- AN L-STABLE METHOD,2 stages,order 2
1953!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1954
1955   double precision g
1956
1957    g = 1.0_dp + 1.0_dp/sqrt(2.0_dp)
1958    rosmethod = rs2
1959!~~~> name of the method
1960    ros_Name = 'ROS-2'
1961!~~~> number of stages
1962    ros_s = 2
1963
1964!~~~> the coefficient matrices a and c are strictly lower triangular.
1965!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1966!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
1967!   The general mapping formula is:
1968!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
1969!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
1970
1971    ros_a(1) = (1.0_dp) /g
1972    ros_c(1) = (- 2.0_dp) /g
1973!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
1974!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
1975    ros_newf(1) = .TRUE.
1976    ros_newf(2) = .TRUE.
1977!~~~> m_i = coefficients for new step solution
1978    ros_m(1) = (3.0_dp) /(2.0_dp* g)
1979    ros_m(2) = (1.0_dp) /(2.0_dp* g)
1980! E_i = Coefficients for error estimator
1981    ros_e(1) = 1.0_dp/(2.0_dp* g)
1982    ros_e(2) = 1.0_dp/(2.0_dp* g)
1983!~~~> ros_elo = estimator of local order - the minimum between the
1984!    main and the embedded scheme orders plus one
1985    ros_elo = 2.0_dp
1986!~~~> y_stage_i ~ y( t + h* alpha_i)
1987    ros_alpha(1) = 0.0_dp
1988    ros_alpha(2) = 1.0_dp
1989!~~~> gamma_i = \sum_j  gamma_{i, j}
1990    ros_gamma(1) = g
1991    ros_gamma(2) = -g
1992
1993 END SUBROUTINE ros2
1994
1995
1996!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1997  SUBROUTINE ros3
1998!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1999! --- AN L-STABLE METHOD,3 stages,order 3,2 function evaluations
2000!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2001
2002   rosmethod = rs3
2003!~~~> name of the method
2004   ros_Name = 'ROS-3'
2005!~~~> number of stages
2006   ros_s = 3
2007
2008!~~~> the coefficient matrices a and c are strictly lower triangular.
2009!   The lower triangular (subdiagonal) elements are stored in row-wise order:
2010!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
2011!   The general mapping formula is:
2012!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
2013!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
2014
2015   ros_a(1) = 1.0_dp
2016   ros_a(2) = 1.0_dp
2017   ros_a(3) = 0.0_dp
2018
2019   ros_c(1) = - 0.10156171083877702091975600115545e+01_dp
2020   ros_c(2) =  0.40759956452537699824805835358067e+01_dp
2021   ros_c(3) =  0.92076794298330791242156818474003e+01_dp
2022!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
2023!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
2024   ros_newf(1) = .TRUE.
2025   ros_newf(2) = .TRUE.
2026   ros_newf(3) = .FALSE.
2027!~~~> m_i = coefficients for new step solution
2028   ros_m(1) =  0.1e+01_dp
2029   ros_m(2) =  0.61697947043828245592553615689730e+01_dp
2030   ros_m(3) = - 0.42772256543218573326238373806514_dp
2031! E_i = Coefficients for error estimator
2032   ros_e(1) =  0.5_dp
2033   ros_e(2) = - 0.29079558716805469821718236208017e+01_dp
2034   ros_e(3) =  0.22354069897811569627360909276199_dp
2035!~~~> ros_elo = estimator of local order - the minimum between the
2036!    main and the embedded scheme orders plus 1
2037   ros_elo = 3.0_dp
2038!~~~> y_stage_i ~ y( t + h* alpha_i)
2039   ros_alpha(1) = 0.0_dp
2040   ros_alpha(2) = 0.43586652150845899941601945119356_dp
2041   ros_alpha(3) = 0.43586652150845899941601945119356_dp
2042!~~~> gamma_i = \sum_j  gamma_{i, j}
2043   ros_gamma(1) = 0.43586652150845899941601945119356_dp
2044   ros_gamma(2) = 0.24291996454816804366592249683314_dp
2045   ros_gamma(3) = 0.21851380027664058511513169485832e+01_dp
2046
2047  END SUBROUTINE ros3
2048
2049!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2050
2051
2052!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2053  SUBROUTINE ros4
2054!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2055!     L-STABLE ROSENBROCK METHOD OF ORDER 4,WITH 4 STAGES
2056!     L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
2057!
2058!      E. HAIRER AND G. WANNER,SOLVING ORDINARY DIFFERENTIAL
2059!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
2060!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
2061!      SPRINGER-VERLAG (1990)
2062!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2063
2064
2065   rosmethod = rs4
2066!~~~> name of the method
2067   ros_Name = 'ROS-4'
2068!~~~> number of stages
2069   ros_s = 4
2070
2071!~~~> the coefficient matrices a and c are strictly lower triangular.
2072!   The lower triangular (subdiagonal) elements are stored in row-wise order:
2073!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
2074!   The general mapping formula is:
2075!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
2076!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
2077
2078   ros_a(1) = 0.2000000000000000e+01_dp
2079   ros_a(2) = 0.1867943637803922e+01_dp
2080   ros_a(3) = 0.2344449711399156_dp
2081   ros_a(4) = ros_a(2)
2082   ros_a(5) = ros_a(3)
2083   ros_a(6) = 0.0_dp
2084
2085   ros_c(1) = -0.7137615036412310e+01_dp
2086   ros_c(2) = 0.2580708087951457e+01_dp
2087   ros_c(3) = 0.6515950076447975_dp
2088   ros_c(4) = -0.2137148994382534e+01_dp
2089   ros_c(5) = -0.3214669691237626_dp
2090   ros_c(6) = -0.6949742501781779_dp
2091!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
2092!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
2093   ros_newf(1) = .TRUE.
2094   ros_newf(2) = .TRUE.
2095   ros_newf(3) = .TRUE.
2096   ros_newf(4) = .FALSE.
2097!~~~> m_i = coefficients for new step solution
2098   ros_m(1) = 0.2255570073418735e+01_dp
2099   ros_m(2) = 0.2870493262186792_dp
2100   ros_m(3) = 0.4353179431840180_dp
2101   ros_m(4) = 0.1093502252409163e+01_dp
2102!~~~> e_i  = coefficients for error estimator
2103   ros_e(1) = -0.2815431932141155_dp
2104   ros_e(2) = -0.7276199124938920e-01_dp
2105   ros_e(3) = -0.1082196201495311_dp
2106   ros_e(4) = -0.1093502252409163e+01_dp
2107!~~~> ros_elo  = estimator of local order - the minimum between the
2108!    main and the embedded scheme orders plus 1
2109   ros_elo  = 4.0_dp
2110!~~~> y_stage_i ~ y( t + h* alpha_i)
2111   ros_alpha(1) = 0.0_dp
2112   ros_alpha(2) = 0.1145640000000000e+01_dp
2113   ros_alpha(3) = 0.6552168638155900_dp
2114   ros_alpha(4) = ros_alpha(3)
2115!~~~> gamma_i = \sum_j  gamma_{i, j}
2116   ros_gamma(1) = 0.5728200000000000_dp
2117   ros_gamma(2) = -0.1769193891319233e+01_dp
2118   ros_gamma(3) = 0.7592633437920482_dp
2119   ros_gamma(4) = -0.1049021087100450_dp
2120
2121  END SUBROUTINE ros4
2122
2123!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2124  SUBROUTINE rodas3
2125!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2126! --- A STIFFLY-STABLE METHOD,4 stages,order 3
2127!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2128
2129
2130   rosmethod = rd3
2131!~~~> name of the method
2132   ros_Name = 'RODAS-3'
2133!~~~> number of stages
2134   ros_s = 4
2135
2136!~~~> the coefficient matrices a and c are strictly lower triangular.
2137!   The lower triangular (subdiagonal) elements are stored in row-wise order:
2138!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
2139!   The general mapping formula is:
2140!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
2141!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
2142
2143   ros_a(1) = 0.0_dp
2144   ros_a(2) = 2.0_dp
2145   ros_a(3) = 0.0_dp
2146   ros_a(4) = 2.0_dp
2147   ros_a(5) = 0.0_dp
2148   ros_a(6) = 1.0_dp
2149
2150   ros_c(1) = 4.0_dp
2151   ros_c(2) = 1.0_dp
2152   ros_c(3) = -1.0_dp
2153   ros_c(4) = 1.0_dp
2154   ros_c(5) = -1.0_dp
2155   ros_c(6) = -(8.0_dp/3.0_dp)
2156
2157!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
2158!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
2159   ros_newf(1) = .TRUE.
2160   ros_newf(2) = .FALSE.
2161   ros_newf(3) = .TRUE.
2162   ros_newf(4) = .TRUE.
2163!~~~> m_i = coefficients for new step solution
2164   ros_m(1) = 2.0_dp
2165   ros_m(2) = 0.0_dp
2166   ros_m(3) = 1.0_dp
2167   ros_m(4) = 1.0_dp
2168!~~~> e_i  = coefficients for error estimator
2169   ros_e(1) = 0.0_dp
2170   ros_e(2) = 0.0_dp
2171   ros_e(3) = 0.0_dp
2172   ros_e(4) = 1.0_dp
2173!~~~> ros_elo  = estimator of local order - the minimum between the
2174!    main and the embedded scheme orders plus 1
2175   ros_elo  = 3.0_dp
2176!~~~> y_stage_i ~ y( t + h* alpha_i)
2177   ros_alpha(1) = 0.0_dp
2178   ros_alpha(2) = 0.0_dp
2179   ros_alpha(3) = 1.0_dp
2180   ros_alpha(4) = 1.0_dp
2181!~~~> gamma_i = \sum_j  gamma_{i, j}
2182   ros_gamma(1) = 0.5_dp
2183   ros_gamma(2) = 1.5_dp
2184   ros_gamma(3) = 0.0_dp
2185   ros_gamma(4) = 0.0_dp
2186
2187  END SUBROUTINE rodas3
2188
2189!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2190  SUBROUTINE rodas4
2191!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2192!     STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4,WITH 6 STAGES
2193!
2194!      E. HAIRER AND G. WANNER,SOLVING ORDINARY DIFFERENTIAL
2195!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
2196!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
2197!      SPRINGER-VERLAG (1996)
2198!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2199
2200
2201    rosmethod = rd4
2202!~~~> name of the method
2203    ros_Name = 'RODAS-4'
2204!~~~> number of stages
2205    ros_s = 6
2206
2207!~~~> y_stage_i ~ y( t + h* alpha_i)
2208    ros_alpha(1) = 0.000_dp
2209    ros_alpha(2) = 0.386_dp
2210    ros_alpha(3) = 0.210_dp
2211    ros_alpha(4) = 0.630_dp
2212    ros_alpha(5) = 1.000_dp
2213    ros_alpha(6) = 1.000_dp
2214
2215!~~~> gamma_i = \sum_j  gamma_{i, j}
2216    ros_gamma(1) = 0.2500000000000000_dp
2217    ros_gamma(2) = -0.1043000000000000_dp
2218    ros_gamma(3) = 0.1035000000000000_dp
2219    ros_gamma(4) = -0.3620000000000023e-01_dp
2220    ros_gamma(5) = 0.0_dp
2221    ros_gamma(6) = 0.0_dp
2222
2223!~~~> the coefficient matrices a and c are strictly lower triangular.
2224!   The lower triangular (subdiagonal) elements are stored in row-wise order:
2225!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
2226!   The general mapping formula is:  A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
2227!                  C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
2228
2229    ros_a(1) = 0.1544000000000000e+01_dp
2230    ros_a(2) = 0.9466785280815826_dp
2231    ros_a(3) = 0.2557011698983284_dp
2232    ros_a(4) = 0.3314825187068521e+01_dp
2233    ros_a(5) = 0.2896124015972201e+01_dp
2234    ros_a(6) = 0.9986419139977817_dp
2235    ros_a(7) = 0.1221224509226641e+01_dp
2236    ros_a(8) = 0.6019134481288629e+01_dp
2237    ros_a(9) = 0.1253708332932087e+02_dp
2238    ros_a(10) = -0.6878860361058950_dp
2239    ros_a(11) = ros_a(7)
2240    ros_a(12) = ros_a(8)
2241    ros_a(13) = ros_a(9)
2242    ros_a(14) = ros_a(10)
2243    ros_a(15) = 1.0_dp
2244
2245    ros_c(1) = -0.5668800000000000e+01_dp
2246    ros_c(2) = -0.2430093356833875e+01_dp
2247    ros_c(3) = -0.2063599157091915_dp
2248    ros_c(4) = -0.1073529058151375_dp
2249    ros_c(5) = -0.9594562251023355e+01_dp
2250    ros_c(6) = -0.2047028614809616e+02_dp
2251    ros_c(7) = 0.7496443313967647e+01_dp
2252    ros_c(8) = -0.1024680431464352e+02_dp
2253    ros_c(9) = -0.3399990352819905e+02_dp
2254    ros_c(10) = 0.1170890893206160e+02_dp
2255    ros_c(11) = 0.8083246795921522e+01_dp
2256    ros_c(12) = -0.7981132988064893e+01_dp
2257    ros_c(13) = -0.3152159432874371e+02_dp
2258    ros_c(14) = 0.1631930543123136e+02_dp
2259    ros_c(15) = -0.6058818238834054e+01_dp
2260
2261!~~~> m_i = coefficients for new step solution
2262    ros_m(1) = ros_a(7)
2263    ros_m(2) = ros_a(8)
2264    ros_m(3) = ros_a(9)
2265    ros_m(4) = ros_a(10)
2266    ros_m(5) = 1.0_dp
2267    ros_m(6) = 1.0_dp
2268
2269!~~~> e_i  = coefficients for error estimator
2270    ros_e(1) = 0.0_dp
2271    ros_e(2) = 0.0_dp
2272    ros_e(3) = 0.0_dp
2273    ros_e(4) = 0.0_dp
2274    ros_e(5) = 0.0_dp
2275    ros_e(6) = 1.0_dp
2276
2277!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
2278!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
2279    ros_newf(1) = .TRUE.
2280    ros_newf(2) = .TRUE.
2281    ros_newf(3) = .TRUE.
2282    ros_newf(4) = .TRUE.
2283    ros_newf(5) = .TRUE.
2284    ros_newf(6) = .TRUE.
2285
2286!~~~> ros_elo  = estimator of local order - the minimum between the
2287!        main and the embedded scheme orders plus 1
2288    ros_elo = 4.0_dp
2289
2290  END SUBROUTINE rodas4
2291!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2292  SUBROUTINE rang3
2293!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2294! STIFFLY-STABLE W METHOD OF ORDER 3,WITH 4 STAGES
2295!
2296! J. RANG and L. ANGERMANN
2297! NEW ROSENBROCK W-METHODS OF ORDER 3
2298! FOR PARTIAL DIFFERENTIAL ALGEBRAIC
2299!        EQUATIONS OF INDEX 1                                             
2300! BIT Numerical Mathematics (2005) 45: 761-787
2301!  DOI: 10.1007/s10543-005-0035-y
2302! Table 4.1-4.2
2303!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2304
2305
2306    rosmethod = rg3
2307!~~~> name of the method
2308    ros_Name = 'RANG-3'
2309!~~~> number of stages
2310    ros_s = 4
2311
2312    ros_a(1) = 5.09052051067020d+00;
2313    ros_a(2) = 5.09052051067020d+00;
2314    ros_a(3) = 0.0d0;
2315    ros_a(4) = 4.97628111010787d+00;
2316    ros_a(5) = 2.77268164715849d-02;
2317    ros_a(6) = 2.29428036027904d-01;
2318
2319    ros_c(1) = - 1.16790812312283d+01;
2320    ros_c(2) = - 1.64057326467367d+01;
2321    ros_c(3) = - 2.77268164715850d-01;
2322    ros_c(4) = - 8.38103960500476d+00;
2323    ros_c(5) = - 8.48328409199343d-01;
2324    ros_c(6) =  2.87009860433106d-01;
2325
2326    ros_m(1) =  5.22582761233094d+00;
2327    ros_m(2) = - 5.56971148154165d-01;
2328    ros_m(3) =  3.57979469353645d-01;
2329    ros_m(4) =  1.72337398521064d+00;
2330
2331    ros_e(1) = - 5.16845212784040d+00;
2332    ros_e(2) = - 1.26351942603842d+00;
2333    ros_e(3) = - 1.11022302462516d-16;
2334    ros_e(4) =  2.22044604925031d-16;
2335
2336    ros_alpha(1) = 0.0d00;
2337    ros_alpha(2) = 2.21878746765329d+00;
2338    ros_alpha(3) = 2.21878746765329d+00;
2339    ros_alpha(4) = 1.55392337535788d+00;
2340
2341    ros_gamma(1) =  4.35866521508459d-01;
2342    ros_gamma(2) = - 1.78292094614483d+00;
2343    ros_gamma(3) = - 2.46541900496934d+00;
2344    ros_gamma(4) = - 8.05529997906370d-01;
2345
2346
2347!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
2348!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
2349    ros_newf(1) = .TRUE.
2350    ros_newf(2) = .TRUE.
2351    ros_newf(3) = .TRUE.
2352    ros_newf(4) = .TRUE.
2353
2354!~~~> ros_elo  = estimator of local order - the minimum between the
2355!        main and the embedded scheme orders plus 1
2356    ros_elo = 3.0_dp
2357
2358  END SUBROUTINE rang3
2359!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2360
2361!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2362!   End of the set of internal Rosenbrock subroutines
2363!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2364END SUBROUTINE rosenbrock
2365 
2366SUBROUTINE funtemplate( t, y, ydot)
2367!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2368!  Template for the ODE function call.
2369!  Updates the rate coefficients (and possibly the fixed species) at each call
2370!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2371!~~~> input variables
2372   REAL(kind=dp):: t, y(nvar)
2373!~~~> output variables
2374   REAL(kind=dp):: ydot(nvar)
2375!~~~> local variables
2376   REAL(kind=dp):: told
2377
2378   told = time
2379   time = t
2380   CALL fun( y, fix, rconst, ydot)
2381   time = told
2382
2383END SUBROUTINE funtemplate
2384 
2385SUBROUTINE jactemplate( t, y, jcb)
2386!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2387!  Template for the ODE Jacobian call.
2388!  Updates the rate coefficients (and possibly the fixed species) at each call
2389!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2390!~~~> input variables
2391    REAL(kind=dp):: t, y(nvar)
2392!~~~> output variables
2393#ifdef full_algebra   
2394    REAL(kind=dp):: jv(lu_nonzero), jcb(nvar, nvar)
2395#else
2396    REAL(kind=dp):: jcb(lu_nonzero)
2397#endif   
2398!~~~> local variables
2399    REAL(kind=dp):: told
2400#ifdef full_algebra   
2401    INTEGER :: i, j
2402#endif   
2403
2404    told = time
2405    time = t
2406#ifdef full_algebra   
2407    CALL jac_sp(y, fix, rconst, jv)
2408    DO j=1, nvar
2409      DO i=1, nvar
2410         jcb(i, j) = 0.0_dp
2411      ENDDO
2412    ENDDO
2413    DO i=1, lu_nonzero
2414       jcb(lu_irow(i), lu_icol(i)) = jv(i)
2415    ENDDO
2416#else
2417    CALL jac_sp( y, fix, rconst, jcb)
2418#endif   
2419    time = told
2420
2421END SUBROUTINE jactemplate
2422 
2423  SUBROUTINE kppdecomp( jvs, ier)                                   
2424  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2425  !        sparse lu factorization                                   
2426  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2427  ! loop expansion generated by kp4                                   
2428                                                                     
2429    INTEGER  :: ier                                                   
2430    REAL(kind=dp):: jvs(lu_nonzero), a                         
2431                                                                     
2432    a = 0.                                                           
2433    ier = 0                                                           
2434                                                                     
2435!   i = 1
2436!   i = 2
2437!   i = 3
2438!   i = 4
2439!   i = 5
2440!   i = 6
2441!   i = 7
2442!   i = 8
2443!   i = 9
2444    jvs(16) = (jvs(16)) / jvs(11) 
2445    jvs(17) = (jvs(17)) / jvs(13) 
2446    jvs(18) = jvs(18) - jvs(12) * jvs(16)
2447    jvs(19) = jvs(19) - jvs(14) * jvs(17)
2448    jvs(20) = jvs(20) - jvs(15) * jvs(17)
2449!   i = 10
2450    jvs(22) = (jvs(22)) / jvs(13) 
2451    jvs(23) = (jvs(23)) / jvs(18) 
2452    jvs(24) = jvs(24) - jvs(14) * jvs(22) - jvs(19) * jvs(23)
2453    jvs(25) = jvs(25) - jvs(15) * jvs(22) - jvs(20) * jvs(23)
2454    jvs(26) = jvs(26) - jvs(21) * jvs(23)
2455!   i = 11
2456    jvs(28) = (jvs(28)) / jvs(13) 
2457    a = 0.0; a = a  - jvs(14) * jvs(28)
2458    jvs(29) = (jvs(29) + a) / jvs(24) 
2459    jvs(30) = jvs(30) - jvs(15) * jvs(28) - jvs(25) * jvs(29)
2460    jvs(31) = jvs(31) - jvs(26) * jvs(29)
2461    jvs(32) = jvs(32) - jvs(27) * jvs(29)
2462!   i = 12
2463    jvs(33) = (jvs(33)) / jvs(30) 
2464    jvs(34) = jvs(34) - jvs(31) * jvs(33)
2465    jvs(35) = jvs(35) - jvs(32) * jvs(33)
2466!   i = 13
2467    jvs(36) = (jvs(36)) / jvs(11) 
2468    a = 0.0; a = a  - jvs(12) * jvs(36)
2469    jvs(37) = (jvs(37) + a) / jvs(18) 
2470    a = 0.0; a = a  - jvs(19) * jvs(37)
2471    jvs(38) = (jvs(38) + a) / jvs(24) 
2472    a = 0.0; a = a  - jvs(20) * jvs(37) - jvs(25) * jvs(38)
2473    jvs(39) = (jvs(39) + a) / jvs(30) 
2474    a = 0.0; a = a  - jvs(21) * jvs(37) - jvs(26) * jvs(38) - jvs(31) * jvs(39)
2475    jvs(40) = (jvs(40) + a) / jvs(34) 
2476    jvs(41) = jvs(41) - jvs(27) * jvs(38) - jvs(32) * jvs(39) - jvs(35) * jvs(40)
2477    RETURN                                                           
2478                                                                     
2479  END SUBROUTINE kppdecomp                                           
2480 
2481SUBROUTINE get_mechanismname                                       
2482                                                                   
2483  IMPLICIT NONE                                                     
2484
2485! Set cs_mech for check with mechanism name from namelist
2486    cs_mech = 'salsa+simple'
2487                                                                   
2488  RETURN                                                           
2489END SUBROUTINE get_mechanismname                                   
2490                                                                   
2491 
2492SUBROUTINE chem_gasphase_integrate (time_step_len, conc, tempi, qvapi, fakti, photo, ierrf, xnacc, xnrej, istatus, l_debug, pe, &
2493                     icntrl_i, rcntrl_i) 
2494                                                                   
2495  IMPLICIT NONE                                                     
2496                                                                   
2497  REAL(dp), INTENT(IN)                  :: time_step_len           
2498  REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: conc                   
2499  REAL(dp), DIMENSION(:, :), INTENT(IN)   :: photo                   
2500  REAL(dp), DIMENSION(:), INTENT(IN)     :: tempi                   
2501  REAL(dp), DIMENSION(:), INTENT(IN)     :: qvapi                   
2502  REAL(dp), DIMENSION(:), INTENT(IN)     :: fakti                   
2503  INTEGER, INTENT(OUT), OPTIONAL        :: ierrf(:)               
2504  INTEGER, INTENT(OUT), OPTIONAL        :: xnacc(:)               
2505  INTEGER, INTENT(OUT), OPTIONAL        :: xnrej(:)               
2506  INTEGER, INTENT(INOUT), OPTIONAL      :: istatus(:)             
2507  INTEGER, INTENT(IN), OPTIONAL         :: pe                     
2508  LOGICAL, INTENT(IN), OPTIONAL         :: l_debug                 
2509  INTEGER, DIMENSION(nkppctrl), INTENT(IN), OPTIONAL  :: icntrl_i         
2510  REAL(dp), DIMENSION(nkppctrl), INTENT(IN), OPTIONAL  :: rcntrl_i         
2511                                                                   
2512  INTEGER                                 :: k   ! loop variable     
2513  REAL(dp)                               :: dt                     
2514  INTEGER, DIMENSION(20)                :: istatus_u               
2515  INTEGER                                :: ierr_u                 
2516  INTEGER                                :: vl_dim_lo               
2517                                                                   
2518                                                                   
2519  IF (PRESENT (istatus)) istatus = 0                             
2520  IF (PRESENT (icntrl_i)) icntrl  = icntrl_i                     
2521  IF (PRESENT (rcntrl_i)) rcntrl  = rcntrl_i                     
2522                                                                   
2523  var => c(1:nvar)                                                 
2524                                                                   
2525  vl_glo = size(tempi, 1)                                           
2526                                                                   
2527  vl_dim_lo = vl_dim                                               
2528  DO k=1, vl_glo, vl_dim_lo                                           
2529    is = k                                                         
2530    ie = min(k+ vl_dim_lo-1, vl_glo)                                 
2531    vl = ie-is+ 1                                                   
2532                                                                   
2533    c(:) = conc(is, :)                                           
2534                                                                   
2535    temp = tempi(is)                                             
2536                                                                   
2537    qvap = qvapi(is)                                             
2538                                                                   
2539    fakt = fakti(is)                                             
2540                                                                   
2541    CALL initialize                                                 
2542                                                                   
2543    phot(:) = photo(is, :)                                           
2544                                                                   
2545    CALL update_rconst                                             
2546                                                                   
2547    dt = time_step_len                                             
2548                                                                   
2549    ! integrate from t=0 to t=dt                                   
2550    CALL integrate(0._dp, dt, icntrl, rcntrl, istatus_u = istatus_u, ierr_u=ierr_u)
2551                                                                   
2552                                                                   
2553   IF (PRESENT(l_debug) .AND. PRESENT(pe)) THEN                       
2554      IF (l_debug) CALL error_output(conc(is, :), ierr_u, pe)           
2555   ENDIF                                                             
2556                                                                     
2557    conc(is, :) = c(:)                                               
2558                                                                   
2559    ! RETURN diagnostic information                                 
2560                                                                   
2561    IF (PRESENT(ierrf))   ierrf(is) = ierr_u                     
2562    IF (PRESENT(xnacc))   xnacc(is) = istatus_u(4)               
2563    IF (PRESENT(xnrej))   xnrej(is) = istatus_u(5)               
2564                                                                   
2565    IF (PRESENT (istatus)) THEN                                   
2566      istatus(1:8) = istatus(1:8) + istatus_u(1:8)               
2567    ENDIF                                                         
2568                                                                   
2569  END DO                                                           
2570 
2571                                                                   
2572! Deallocate input arrays                                           
2573                                                                   
2574                                                                   
2575  data_loaded = .FALSE.                                             
2576                                                                   
2577  RETURN                                                           
2578END SUBROUTINE chem_gasphase_integrate                             
2579
2580END MODULE chem_gasphase_mod
2581 
Note: See TracBrowser for help on using the repository browser.