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

Last change on this file since 3582 was 3566, checked in by monakurppa, 6 years ago

Branch salsa: included chemical mechanism including salsa gases

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