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

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

removed read from unit 10 in chemistry_model_mod.f90, added get_mechanismname

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