source: palm/trunk/SOURCE/chem_gasphase_mod.f90 @ 3436

Last change on this file since 3436 was 3298, checked in by kanani, 6 years ago

Merge chemistry branch at r3297 to trunk

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