source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_salsagas/chem_gasphase_mod.f90 @ 3780

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

removed read from unit 10 in chemistry_model_mod.f90, added get_mechanismname

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