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

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

added leading blanks to dummy statements

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