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

Last change on this file since 3944 was 3944, checked in by maronga, 6 years ago

fix for last commit

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