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

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

Reintegrated fixes/changes from branch chemistry

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