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

Last change on this file since 3833 was 3833, checked in by forkel, 2 years ago

removed USE chem_gasphase_mod from chem_modules, apply USE chem_gasphase for nvar, nspec, cs_mech and spc_names instead

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