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

Last change on this file since 3698 was 3698, checked in by suehring, 3 years ago

Fix bad commit in gasphase_preproc

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