source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_salsa+simple/chem_gasphase_mod.f90 @ 3789

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

Removed unused variables from chem_gasphase_mod.f90

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