source: palm/trunk/TUTORIALS/cases/chem_smog/USER_CODE/chem_gasphase_mod.f90 @ 4782

Last change on this file since 4782 was 4409, checked in by Giersch, 5 years ago

Tutorials updated, all tested with revision 4402

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