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

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

Merge chemistry branch at r3297 to trunk

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