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

Last change on this file was 4841, checked in by forkel, 19 months ago

updated copyright statements for chem_gasphase_mod.f90. This must be done in UTIL/chemistry/gasphase_preproc/kpp4palm/templates/module_header and NOT in SOURCE/chem_gasphase_mod.f90.

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