source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_simplep/chem_gasphase_mod.f90 @ 3789

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

Removed unused variables from chem_gasphase_mod.f90

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