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

Last change on this file since 4370 was 4370, checked in by raasch, 2 years ago

bugfixes for previous commit: unused variables removed, Temperton-fft usage on GPU, openacc porting of vector version of Obukhov length calculation, collective read switched off on NEC to avoid hanging; some vector directives added in prognostic equations to force vectorization on Intel19 compiler, configuration files for NEC Aurora added

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