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

Last change on this file since 3585 was 3585, checked in by forkel, 4 years ago

Comments in all chem_gasphase_mod.kpp, add def_salsagas/*.f90, removed executables

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