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

Last change on this file since 3395 was 3298, checked in by kanani, 5 years ago

Merge chemistry branch at r3297 to trunk

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