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

Last change on this file since 3949 was 3949, checked in by forkel, 6 years ago

fix in def_salsa+simple/chem_gasphase_mod.kpp and finaling def_salsa+phstat

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