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

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

removed read from unit 10 in chemistry_model_mod.f90, added get_mechanismname

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