source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_simple/chem_gasphase_mod.f90 @ 3800

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

added leading blanks to dummy statements

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