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

Last change on this file since 4901 was 4841, checked in by forkel, 19 months ago

updated copyright statements for chem_gasphase_mod.f90. This must be done in UTIL/chemistry/gasphase_preproc/kpp4palm/templates/module_header and NOT in SOURCE/chem_gasphase_mod.f90.

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