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

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

Reintegrated fixes/changes from branch chemistry

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