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