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

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

Branch salsa: included chemical mechanism including salsa gases

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