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