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

Last change on this file since 3789 was 3789, checked in by forkel, 5 years ago

Removed unused variables from chem_gasphase_mod.f90

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