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

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

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

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