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

Last change on this file since 4370 was 4370, checked in by raasch, 23 months ago

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

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