source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_phstat/chem_gasphase_mod.f90 @ 3298

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

Merge chemistry branch at r3297 to trunk

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