source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_smog/chem_gasphase_mod.f90 @ 3298

Last change on this file since 3298 was 3298, checked in by kanani, 3 years ago

Merge chemistry branch at r3297 to trunk

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