source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_passive1/chem_gasphase_mod.f90 @ 3789

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

Removed unused variables from chem_gasphase_mod.f90

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