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

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

fix for last commit

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