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

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