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

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

renaming of get_mechanismname, do_emiss and do_depo, sorting in chem_modules

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