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

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

Removed unused variables from chem_gasphase_mod.f90

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