source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_simplep/chem_gasphase_mod.f90 @ 4841

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