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

Last change on this file was 4841, checked in by forkel, 19 months ago

updated copyright statements for chem_gasphase_mod.f90. This must be done in UTIL/chemistry/gasphase_preproc/kpp4palm/templates/module_header and NOT in SOURCE/chem_gasphase_mod.f90.

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