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

Last change on this file since 3698 was 3698, checked in by suehring, 3 years ago

Fix bad commit in gasphase_preproc

File size: 85.0 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 :: eqn_names, phot_names, spc_names
68  PUBLIC :: nmaxfixsteps
69  PUBLIC :: atol, rtol
70  PUBLIC :: nspec, nreact
71  PUBLIC :: temp
72  PUBLIC :: qvap
73  PUBLIC :: fakt
74  PUBLIC :: phot
75  PUBLIC :: rconst
76  PUBLIC :: nvar
77  PUBLIC :: nphot
78  PUBLIC :: vl_dim                     ! PUBLIC to ebable other MODULEs to distiguish between scalar and vec
79 
80  PUBLIC :: initialize, integrate, update_rconst
81  PUBLIC :: chem_gasphase_integrate
82  PUBLIC :: initialize_kpp_ctrl
83
84! END OF MODULE HEADER TEMPLATE
85                                                                 
86! Variables used for vector mode                                 
87                                                                 
88  LOGICAL, PARAMETER          :: l_vector = .FALSE.           
89  INTEGER, PARAMETER          :: i_lu_di = 2
90  INTEGER, PARAMETER          :: vl_dim = 1
91  INTEGER                     :: vl                             
92                                                                 
93  INTEGER                     :: vl_glo                         
94  INTEGER                     :: is, ie                           
95                                                                 
96                                                                 
97  INTEGER, DIMENSION(vl_dim)  :: kacc, krej                       
98  INTEGER, DIMENSION(vl_dim)  :: ierrv                           
99  LOGICAL                     :: data_loaded = .FALSE.             
100! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
101!
102! Parameter Module File
103!
104! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
105!       (http://www.cs.vt.edu/~asandu/Software/KPP)
106! KPP is distributed under GPL,the general public licence
107!       (http://www.gnu.org/copyleft/gpl.html)
108! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
109! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
110!     With important contributions from:
111!        M. Damian,Villanova University,USA
112!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
113!
114! File                 : chem_gasphase_mod_Parameters.f90
115! Time                 : Thu Dec 20 14:57:53 2018
116! Working directory    : /home/forkel-r/palmstuff/work/trunk20181220/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                 : Thu Dec 20 14:57:53 2018
201! Working directory    : /home/forkel-r/palmstuff/work/trunk20181220/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! TSTART - Integration start time
230  REAL(kind=dp):: tstart
231! ATOL - Absolute tolerance
232  REAL(kind=dp):: atol(nvar)
233! RTOL - Relative tolerance
234  REAL(kind=dp):: rtol(nvar)
235! STEPMIN - Lower bound for integration step
236  REAL(kind=dp):: stepmin
237! CFACTOR - Conversion factor for concentration units
238  REAL(kind=dp):: cfactor
239
240! INLINED global variable declarations
241
242! QVAP - Water vapor
243  REAL(kind=dp):: qvap
244! FAKT - Conversion factor
245  REAL(kind=dp):: fakt
246
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                 : Thu Dec 20 14:57:53 2018
268! Working directory    : /home/forkel-r/palmstuff/work/trunk20181220/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                 : Thu Dec 20 14:57:53 2018
318! Working directory    : /home/forkel-r/palmstuff/work/trunk20181220/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                 : Thu Dec 20 14:57:53 2018
388! Working directory    : /home/forkel-r/palmstuff/work/trunk20181220/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                 : Thu Dec 20 14:57:53 2018
414! Working directory    : /home/forkel-r/palmstuff/work/trunk20181220/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                 : Thu Dec 20 14:57:53 2018
472! Working directory    : /home/forkel-r/palmstuff/work/trunk20181220/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                 : Thu Dec 20 14:57:53 2018
499! Working directory    : /home/forkel-r/palmstuff/work/trunk20181220/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                 : Thu Dec 20 14:57:53 2018
526! Working directory    : /home/forkel-r/palmstuff/work/trunk20181220/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                 : Thu Dec 20 14:57:53 2018
555! Working directory    : /home/forkel-r/palmstuff/work/trunk20181220/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                 : Thu Dec 20 14:57:53 2018
581! Working directory    : /home/forkel-r/palmstuff/work/trunk20181220/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            chem_gasphase_integrate
688    MODULE PROCEDURE   chem_gasphase_integrate
689  END INTERFACE        chem_gasphase_integrate
690 
691 
692 CONTAINS
693 
694SUBROUTINE initialize()
695
696
697  INTEGER         :: j, k
698
699  INTEGER :: i
700  REAL(kind=dp):: x
701  k = is
702  cfactor = 1.000000e+00_dp
703
704  x = (0.) * cfactor
705  DO i = 1 , nvar
706  ENDDO
707
708  x = (0.) * cfactor
709  DO i = 1 , nfix
710    fix(i) = x
711  ENDDO
712
713! constant rate coefficients
714! END constant rate coefficients
715
716! INLINED initializations
717
718! End INLINED initializations
719
720     
721END SUBROUTINE initialize
722 
723SUBROUTINE integrate( tin, tout, &
724  icntrl_u, rcntrl_u, istatus_u, rstatus_u, ierr_u)
725
726
727   REAL(kind=dp), INTENT(IN):: tin  ! start time
728   REAL(kind=dp), INTENT(IN):: tout ! END time
729   ! OPTIONAL input PARAMETERs and statistics
730   INTEGER,      INTENT(IN), OPTIONAL :: icntrl_u(20)
731   REAL(kind=dp), INTENT(IN), OPTIONAL :: rcntrl_u(20)
732   INTEGER,      INTENT(OUT), OPTIONAL :: istatus_u(20)
733   REAL(kind=dp), INTENT(OUT), OPTIONAL :: rstatus_u(20)
734   INTEGER,      INTENT(OUT), OPTIONAL :: ierr_u
735
736   REAL(kind=dp):: rcntrl(20), rstatus(20)
737   INTEGER       :: icntrl(20), istatus(20), ierr
738
739   INTEGER, SAVE :: ntotal = 0
740
741   icntrl(:) = 0
742   rcntrl(:) = 0.0_dp
743   istatus(:) = 0
744   rstatus(:) = 0.0_dp
745
746    !~~~> fine-tune the integrator:
747   icntrl(1) = 0      ! 0 - non- autonomous, 1 - autonomous
748   icntrl(2) = 0      ! 0 - vector tolerances, 1 - scalars
749
750   ! IF OPTIONAL PARAMETERs are given, and IF they are >0,
751   ! THEN they overwrite default settings.
752   IF (PRESENT(icntrl_u))THEN
753     WHERE(icntrl_u(:)> 0)icntrl(:) = icntrl_u(:)
754   ENDIF
755   IF (PRESENT(rcntrl_u))THEN
756     WHERE(rcntrl_u(:)> 0)rcntrl(:) = rcntrl_u(:)
757   ENDIF
758
759
760   CALL rosenbrock(nvar, var, tin, tout,  &
761         atol, rtol,               &
762         rcntrl, icntrl, rstatus, istatus, ierr)
763
764   !~~~> debug option: show no of steps
765   ! ntotal = ntotal + istatus(nstp)
766   ! PRINT*,'NSTEPS=',ISTATUS(Nstp),' (',Ntotal,')','  O3=',VAR(ind_O3)
767
768   stepmin = rstatus(nhexit)
769   ! IF OPTIONAL PARAMETERs are given for output they
770   ! are updated with the RETURN information
771   IF (PRESENT(istatus_u))istatus_u(:) = istatus(:)
772   IF (PRESENT(rstatus_u))rstatus_u(:) = rstatus(:)
773   IF (PRESENT(ierr_u))  ierr_u       = ierr
774
775END SUBROUTINE integrate
776 
777SUBROUTINE fun(v, f, rct, vdot)
778
779! V - Concentrations of variable species (local)
780  REAL(kind=dp):: v(nvar)
781! F - Concentrations of fixed species (local)
782  REAL(kind=dp):: f(nfix)
783! RCT - Rate constants (local)
784  REAL(kind=dp):: rct(nreact)
785! Vdot - Time derivative of variable species concentrations
786  REAL(kind=dp):: vdot(nvar)
787
788
789! Computation of equation rates
790  a(1) = rct(1) * v(10)
791  a(2) = rct(2) * v(6) * f(1)
792  a(3) = rct(3) * v(6) * v(9)
793  a(4) = rct(4) * v(4) * v(7)
794  a(5) = rct(5) * v(8) * v(9)
795  a(6) = rct(6) * v(5) * v(9)
796  a(7) = rct(7) * v(7) * v(10)
797
798! Aggregate function
799  vdot(1) = 0
800  vdot(2) = a(7)
801  vdot(3) = a(5)
802  vdot(4) = - a(4)
803  vdot(5) = a(5) - a(6)
804  vdot(6) = a(1) - a(2) - a(3)
805  vdot(7) = 2* a(2) - a(4) + a(6) - a(7)
806  vdot(8) = a(4) - a(5)
807  vdot(9) = a(1) - a(3) - a(5) - a(6)
808  vdot(10) = - a(1) + a(3) + a(5) + a(6) - a(7)
809     
810END SUBROUTINE fun
811 
812SUBROUTINE kppsolve(jvs, x)
813
814! JVS - sparse Jacobian of variables
815  REAL(kind=dp):: jvs(lu_nonzero)
816! X - Vector for variables
817  REAL(kind=dp):: x(nvar)
818
819  x(7) = x(7) - jvs(16) * x(4) - jvs(17) * x(5) - jvs(18) * x(6)
820  x(8) = x(8) - jvs(23) * x(4) - jvs(24) * x(7)
821  x(9) = x(9) - jvs(28) * x(5) - jvs(29) * x(6) - jvs(30) * x(8)
822  x(10) = x(10) - jvs(33) * x(5) - jvs(34) * x(6) - jvs(35) * x(7) - jvs(36) * x(8) - jvs(37) * x(9)
823  x(10) = x(10) / jvs(38)
824  x(9) = (x(9) - jvs(32) * x(10)) /(jvs(31))
825  x(8) = (x(8) - jvs(26) * x(9) - jvs(27) * x(10)) /(jvs(25))
826  x(7) = (x(7) - jvs(20) * x(8) - jvs(21) * x(9) - jvs(22) * x(10)) /(jvs(19))
827  x(6) = (x(6) - jvs(14) * x(9) - jvs(15) * x(10)) /(jvs(13))
828  x(5) = (x(5) - jvs(11) * x(8) - jvs(12) * x(9)) /(jvs(10))
829  x(4) = (x(4) - jvs(9) * x(7)) /(jvs(8))
830  x(3) = (x(3) - jvs(6) * x(8) - jvs(7) * x(9)) /(jvs(5))
831  x(2) = (x(2) - jvs(3) * x(7) - jvs(4) * x(10)) /(jvs(2))
832  x(1) = x(1) / jvs(1)
833     
834END SUBROUTINE kppsolve
835 
836SUBROUTINE jac_sp(v, f, rct, jvs)
837
838! V - Concentrations of variable species (local)
839  REAL(kind=dp):: v(nvar)
840! F - Concentrations of fixed species (local)
841  REAL(kind=dp):: f(nfix)
842! RCT - Rate constants (local)
843  REAL(kind=dp):: rct(nreact)
844! JVS - sparse Jacobian of variables
845  REAL(kind=dp):: jvs(lu_nonzero)
846
847
848! Local variables
849! B - Temporary array
850  REAL(kind=dp):: b(14)
851
852! B(1) = dA(1)/dV(10)
853  b(1) = rct(1)
854! B(2) = dA(2)/dV(6)
855  b(2) = rct(2) * f(1)
856! B(4) = dA(3)/dV(6)
857  b(4) = rct(3) * v(9)
858! B(5) = dA(3)/dV(9)
859  b(5) = rct(3) * v(6)
860! B(6) = dA(4)/dV(4)
861  b(6) = rct(4) * v(7)
862! B(7) = dA(4)/dV(7)
863  b(7) = rct(4) * v(4)
864! B(8) = dA(5)/dV(8)
865  b(8) = rct(5) * v(9)
866! B(9) = dA(5)/dV(9)
867  b(9) = rct(5) * v(8)
868! B(10) = dA(6)/dV(5)
869  b(10) = rct(6) * v(9)
870! B(11) = dA(6)/dV(9)
871  b(11) = rct(6) * v(5)
872! B(12) = dA(7)/dV(7)
873  b(12) = rct(7) * v(10)
874! B(13) = dA(7)/dV(10)
875  b(13) = rct(7) * v(7)
876! B(14) = dA(8)/dV(1)
877  b(14) = rct(8)
878
879! Construct the Jacobian terms from B's
880! JVS(1) = Jac_FULL(1,1)
881  jvs(1) = 0
882! JVS(2) = Jac_FULL(2,2)
883  jvs(2) = 0
884! JVS(3) = Jac_FULL(2,7)
885  jvs(3) = b(12)
886! JVS(4) = Jac_FULL(2,10)
887  jvs(4) = b(13)
888! JVS(5) = Jac_FULL(3,3)
889  jvs(5) = 0
890! JVS(6) = Jac_FULL(3,8)
891  jvs(6) = b(8)
892! JVS(7) = Jac_FULL(3,9)
893  jvs(7) = b(9)
894! JVS(8) = Jac_FULL(4,4)
895  jvs(8) = - b(6)
896! JVS(9) = Jac_FULL(4,7)
897  jvs(9) = - b(7)
898! JVS(10) = Jac_FULL(5,5)
899  jvs(10) = - b(10)
900! JVS(11) = Jac_FULL(5,8)
901  jvs(11) = b(8)
902! JVS(12) = Jac_FULL(5,9)
903  jvs(12) = b(9) - b(11)
904! JVS(13) = Jac_FULL(6,6)
905  jvs(13) = - b(2) - b(4)
906! JVS(14) = Jac_FULL(6,9)
907  jvs(14) = - b(5)
908! JVS(15) = Jac_FULL(6,10)
909  jvs(15) = b(1)
910! JVS(16) = Jac_FULL(7,4)
911  jvs(16) = - b(6)
912! JVS(17) = Jac_FULL(7,5)
913  jvs(17) = b(10)
914! JVS(18) = Jac_FULL(7,6)
915  jvs(18) = 2* b(2)
916! JVS(19) = Jac_FULL(7,7)
917  jvs(19) = - b(7) - b(12)
918! JVS(20) = Jac_FULL(7,8)
919  jvs(20) = 0
920! JVS(21) = Jac_FULL(7,9)
921  jvs(21) = b(11)
922! JVS(22) = Jac_FULL(7,10)
923  jvs(22) = - b(13)
924! JVS(23) = Jac_FULL(8,4)
925  jvs(23) = b(6)
926! JVS(24) = Jac_FULL(8,7)
927  jvs(24) = b(7)
928! JVS(25) = Jac_FULL(8,8)
929  jvs(25) = - b(8)
930! JVS(26) = Jac_FULL(8,9)
931  jvs(26) = - b(9)
932! JVS(27) = Jac_FULL(8,10)
933  jvs(27) = 0
934! JVS(28) = Jac_FULL(9,5)
935  jvs(28) = - b(10)
936! JVS(29) = Jac_FULL(9,6)
937  jvs(29) = - b(4)
938! JVS(30) = Jac_FULL(9,8)
939  jvs(30) = - b(8)
940! JVS(31) = Jac_FULL(9,9)
941  jvs(31) = - b(5) - b(9) - b(11)
942! JVS(32) = Jac_FULL(9,10)
943  jvs(32) = b(1)
944! JVS(33) = Jac_FULL(10,5)
945  jvs(33) = b(10)
946! JVS(34) = Jac_FULL(10,6)
947  jvs(34) = b(4)
948! JVS(35) = Jac_FULL(10,7)
949  jvs(35) = - b(12)
950! JVS(36) = Jac_FULL(10,8)
951  jvs(36) = b(8)
952! JVS(37) = Jac_FULL(10,9)
953  jvs(37) = b(5) + b(9) + b(11)
954! JVS(38) = Jac_FULL(10,10)
955  jvs(38) = - b(1) - b(13)
956     
957END SUBROUTINE jac_sp
958 
959  elemental REAL(kind=dp)FUNCTION k_arr (k_298, tdep, temp)
960    ! arrhenius FUNCTION
961
962    REAL,    INTENT(IN):: k_298 ! k at t = 298.15k
963    REAL,    INTENT(IN):: tdep  ! temperature dependence
964    REAL(kind=dp), INTENT(IN):: temp  ! temperature
965
966    intrinsic exp
967
968    k_arr = k_298 * exp(tdep* (1._dp/temp- 3.3540e-3_dp))! 1/298.15=3.3540e-3
969
970  END FUNCTION k_arr
971 
972SUBROUTINE update_rconst()
973 INTEGER         :: k
974
975  k = is
976
977! Begin INLINED RCONST
978
979
980! End INLINED RCONST
981
982  rconst(1) = (phot(j_no2))
983  rconst(2) = (2.0_dp * 2.2e-10_dp * phot(j_o31d) / (arr2(1.9e+8_dp , -390.0_dp , temp)))
984  rconst(3) = (arr2(1.8e-12_dp , 1370.0_dp , temp))
985  rconst(4) = (arr2(2.e-11_dp , 500.0_dp , temp))
986  rconst(5) = (arr2(4.2e-12_dp , -180.0_dp , temp))
987  rconst(6) = (arr2(3.7e-12_dp , -240.0_dp , temp))
988  rconst(7) = (arr2(1.15e-11_dp , 0.0_dp , temp))
989  rconst(8) = (1.0_dp)
990     
991END SUBROUTINE update_rconst
992 
993!  END FUNCTION ARR2
994REAL(kind=dp)FUNCTION arr2( a0, b0, temp)
995    REAL(kind=dp):: temp
996    REAL(kind=dp):: a0, b0
997    arr2 = a0 * exp( - b0 / temp)
998END FUNCTION arr2
999 
1000SUBROUTINE initialize_kpp_ctrl(status)
1001
1002
1003  ! i/o
1004  INTEGER,         INTENT(OUT):: status
1005
1006  ! local
1007  REAL(dp):: tsum
1008  INTEGER  :: i
1009
1010  ! check fixed time steps
1011  tsum = 0.0_dp
1012  DO i=1, nmaxfixsteps
1013     IF (t_steps(i)< tiny(0.0_dp))exit
1014     tsum = tsum + t_steps(i)
1015  ENDDO
1016
1017  nfsteps = i- 1
1018
1019  l_fixed_step = (nfsteps > 0).and.((tsum - 1.0)< tiny(0.0_dp))
1020
1021  IF (l_vector)THEN
1022     WRITE(*,*) ' MODE             : VECTOR (LENGTH=',VL_DIM,')'
1023  ELSE
1024     WRITE(*,*) ' MODE             : SCALAR'
1025  ENDIF
1026  !
1027  WRITE(*,*) ' DE-INDEXING MODE :',I_LU_DI
1028  !
1029  WRITE(*,*) ' ICNTRL           : ',icntrl
1030  WRITE(*,*) ' RCNTRL           : ',rcntrl
1031  !
1032  ! note: this is ONLY meaningful for vectorized (kp4)rosenbrock- methods
1033  IF (l_vector)THEN
1034     IF (l_fixed_step)THEN
1035        WRITE(*,*) ' TIME STEPS       : FIXED (',t_steps(1:nfsteps),')'
1036     ELSE
1037        WRITE(*,*) ' TIME STEPS       : AUTOMATIC'
1038     ENDIF
1039  ELSE
1040     WRITE(*,*) ' TIME STEPS       : AUTOMATIC '//&
1041          &'(t_steps (CTRL_KPP) ignored in SCALAR MODE)'
1042  ENDIF
1043  ! mz_pj_20070531-
1044
1045  status = 0
1046
1047
1048END SUBROUTINE initialize_kpp_ctrl
1049 
1050SUBROUTINE error_output(c, ierr, pe)
1051
1052
1053  INTEGER, INTENT(IN):: ierr
1054  INTEGER, INTENT(IN):: pe
1055  REAL(dp), DIMENSION(:), INTENT(IN):: c
1056
1057  write(6,*) 'ERROR in chem_gasphase_mod ',ierr,C(1)
1058
1059
1060END SUBROUTINE error_output
1061 
1062      SUBROUTINE wscal(n, alpha, x, incx)
1063!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1064!     constant times a vector: x(1:N) <- Alpha*x(1:N)
1065!     only for incX=incY=1
1066!     after BLAS
1067!     replace this by the function from the optimized BLAS implementation:
1068!         CALL SSCAL(N,Alpha,X,1) or  CALL DSCAL(N,Alpha,X,1)
1069!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1070
1071      INTEGER  :: i, incx, m, mp1, n
1072      REAL(kind=dp) :: x(n), alpha
1073      REAL(kind=dp), PARAMETER  :: zero=0.0_dp, one=1.0_dp
1074
1075      IF (alpha .eq. one)RETURN
1076      IF (n .le. 0)RETURN
1077
1078      m = mod(n, 5)
1079      IF ( m .ne. 0)THEN
1080        IF (alpha .eq. (- one))THEN
1081          DO i = 1, m
1082            x(i) = - x(i)
1083          ENDDO
1084        ELSEIF (alpha .eq. zero)THEN
1085          DO i = 1, m
1086            x(i) = zero
1087          ENDDO
1088        ELSE
1089          DO i = 1, m
1090            x(i) = alpha* x(i)
1091          ENDDO
1092        ENDIF
1093        IF ( n .lt. 5)RETURN
1094      ENDIF
1095      mp1 = m + 1
1096      IF (alpha .eq. (- one))THEN
1097        DO i = mp1, n, 5
1098          x(i)   = - x(i)
1099          x(i + 1) = - x(i + 1)
1100          x(i + 2) = - x(i + 2)
1101          x(i + 3) = - x(i + 3)
1102          x(i + 4) = - x(i + 4)
1103        ENDDO
1104      ELSEIF (alpha .eq. zero)THEN
1105        DO i = mp1, n, 5
1106          x(i)   = zero
1107          x(i + 1) = zero
1108          x(i + 2) = zero
1109          x(i + 3) = zero
1110          x(i + 4) = zero
1111        ENDDO
1112      ELSE
1113        DO i = mp1, n, 5
1114          x(i)   = alpha* x(i)
1115          x(i + 1) = alpha* x(i + 1)
1116          x(i + 2) = alpha* x(i + 2)
1117          x(i + 3) = alpha* x(i + 3)
1118          x(i + 4) = alpha* x(i + 4)
1119        ENDDO
1120      ENDIF
1121
1122      END SUBROUTINE wscal
1123 
1124      SUBROUTINE waxpy(n, alpha, x, incx, y, incy)
1125!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1126!     constant times a vector plus a vector: y <- y + Alpha*x
1127!     only for incX=incY=1
1128!     after BLAS
1129!     replace this by the function from the optimized BLAS implementation:
1130!         CALL SAXPY(N,Alpha,X,1,Y,1) or  CALL DAXPY(N,Alpha,X,1,Y,1)
1131!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1132
1133      INTEGER  :: i, incx, incy, m, mp1, n
1134      REAL(kind=dp):: x(n), y(n), alpha
1135      REAL(kind=dp), PARAMETER :: zero = 0.0_dp
1136
1137      IF (alpha .eq. zero)RETURN
1138      IF (n .le. 0)RETURN
1139
1140      m = mod(n, 4)
1141      IF ( m .ne. 0)THEN
1142        DO i = 1, m
1143          y(i) = y(i) + alpha* x(i)
1144        ENDDO
1145        IF ( n .lt. 4)RETURN
1146      ENDIF
1147      mp1 = m + 1
1148      DO i = mp1, n, 4
1149        y(i) = y(i) + alpha* x(i)
1150        y(i + 1) = y(i + 1) + alpha* x(i + 1)
1151        y(i + 2) = y(i + 2) + alpha* x(i + 2)
1152        y(i + 3) = y(i + 3) + alpha* x(i + 3)
1153      ENDDO
1154     
1155      END SUBROUTINE waxpy
1156 
1157SUBROUTINE rosenbrock(n, y, tstart, tend, &
1158           abstol, reltol,             &
1159           rcntrl, icntrl, rstatus, istatus, ierr)
1160!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1161!
1162!    Solves the system y'=F(t,y) using a Rosenbrock method defined by:
1163!
1164!     G = 1/(H*gamma(1)) - Jac(t0,Y0)
1165!     T_i = t0 + Alpha(i)*H
1166!     Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
1167!     G *K_i = Fun( T_i,Y_i)+ \sum_{j=1}^S C(i,j)/H *K_j +
1168!         gamma(i)*dF/dT(t0,Y0)
1169!     Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
1170!
1171!    For details on Rosenbrock methods and their implementation consult:
1172!      E. Hairer and G. Wanner
1173!      "Solving ODEs II. Stiff and differential-algebraic problems".
1174!      Springer series in computational mathematics,Springer-Verlag,1996.
1175!    The codes contained in the book inspired this implementation.
1176!
1177!    (C)  Adrian Sandu,August 2004
1178!    Virginia Polytechnic Institute and State University
1179!    Contact: sandu@cs.vt.edu
1180!    Revised by Philipp Miehe and Adrian Sandu,May 2006                 
1181!    This implementation is part of KPP - the Kinetic PreProcessor
1182!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1183!
1184!~~~>   input arguments:
1185!
1186!-     y(n)  = vector of initial conditions (at t=tstart)
1187!-    [tstart, tend]  = time range of integration
1188!     (if Tstart>Tend the integration is performed backwards in time)
1189!-    reltol, abstol = user precribed accuracy
1190!- SUBROUTINE  fun( t, y, ydot) = ode FUNCTION,
1191!                       returns Ydot = Y' = F(T,Y)
1192!- SUBROUTINE  jac( t, y, jcb) = jacobian of the ode FUNCTION,
1193!                       returns Jcb = dFun/dY
1194!-    icntrl(1:20)  = INTEGER inputs PARAMETERs
1195!-    rcntrl(1:20)  = REAL inputs PARAMETERs
1196!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1197!
1198!~~~>     output arguments:
1199!
1200!-    y(n)  - > vector of final states (at t- >tend)
1201!-    istatus(1:20) - > INTEGER output PARAMETERs
1202!-    rstatus(1:20) - > REAL output PARAMETERs
1203!-    ierr            - > job status upon RETURN
1204!                        success (positive value) or
1205!                        failure (negative value)
1206!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1207!
1208!~~~>     input PARAMETERs:
1209!
1210!    Note: For input parameters equal to zero the default values of the
1211!       corresponding variables are used.
1212!
1213!    ICNTRL(1) = 1: F = F(y)   Independent of T (AUTONOMOUS)
1214!              = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
1215!
1216!    ICNTRL(2) = 0: AbsTol,RelTol are N-dimensional vectors
1217!              = 1: AbsTol,RelTol are scalars
1218!
1219!    ICNTRL(3)  -> selection of a particular Rosenbrock method
1220!        = 0 :    Rodas3 (default)
1221!        = 1 :    Ros2
1222!        = 2 :    Ros3
1223!        = 3 :    Ros4
1224!        = 4 :    Rodas3
1225!        = 5 :    Rodas4
1226!
1227!    ICNTRL(4)  -> maximum number of integration steps
1228!        For ICNTRL(4) =0) the default value of 100000 is used
1229!
1230!    RCNTRL(1)  -> Hmin,lower bound for the integration step size
1231!          It is strongly recommended to keep Hmin = ZERO
1232!    RCNTRL(2)  -> Hmax,upper bound for the integration step size
1233!    RCNTRL(3)  -> Hstart,starting value for the integration step size
1234!
1235!    RCNTRL(4)  -> FacMin,lower bound on step decrease factor (default=0.2)
1236!    RCNTRL(5)  -> FacMax,upper bound on step increase factor (default=6)
1237!    RCNTRL(6)  -> FacRej,step decrease factor after multiple rejections
1238!                          (default=0.1)
1239!    RCNTRL(7)  -> FacSafe,by which the new step is slightly smaller
1240!         than the predicted value  (default=0.9)
1241!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1242!
1243!
1244!    OUTPUT ARGUMENTS:
1245!    -----------------
1246!
1247!    T           -> T value for which the solution has been computed
1248!                     (after successful return T=Tend).
1249!
1250!    Y(N)        -> Numerical solution at T
1251!
1252!    IDID        -> Reports on successfulness upon return:
1253!                    = 1 for success
1254!                    < 0 for error (value equals error code)
1255!
1256!    ISTATUS(1)  -> No. of function calls
1257!    ISTATUS(2)  -> No. of jacobian calls
1258!    ISTATUS(3)  -> No. of steps
1259!    ISTATUS(4)  -> No. of accepted steps
1260!    ISTATUS(5)  -> No. of rejected steps (except at very beginning)
1261!    ISTATUS(6)  -> No. of LU decompositions
1262!    ISTATUS(7)  -> No. of forward/backward substitutions
1263!    ISTATUS(8)  -> No. of singular matrix decompositions
1264!
1265!    RSTATUS(1)  -> Texit,the time corresponding to the
1266!                     computed Y upon return
1267!    RSTATUS(2)  -> Hexit,last accepted step before exit
1268!    RSTATUS(3)  -> Hnew,last predicted step (not yet taken)
1269!                   For multiple restarts,use Hnew as Hstart
1270!                     in the subsequent run
1271!
1272!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1273
1274
1275!~~~>  arguments
1276   INTEGER,      INTENT(IN)  :: n
1277   REAL(kind=dp), INTENT(INOUT):: y(n)
1278   REAL(kind=dp), INTENT(IN)  :: tstart, tend
1279   REAL(kind=dp), INTENT(IN)  :: abstol(n), reltol(n)
1280   INTEGER,      INTENT(IN)  :: icntrl(20)
1281   REAL(kind=dp), INTENT(IN)  :: rcntrl(20)
1282   INTEGER,      INTENT(INOUT):: istatus(20)
1283   REAL(kind=dp), INTENT(INOUT):: rstatus(20)
1284   INTEGER, INTENT(OUT) :: ierr
1285!~~~>  PARAMETERs of the rosenbrock method, up to 6 stages
1286   INTEGER ::  ros_s, rosmethod
1287   INTEGER, PARAMETER :: rs2=1, rs3=2, rs4=3, rd3=4, rd4=5, rg3=6
1288   REAL(kind=dp):: ros_a(15), ros_c(15), ros_m(6), ros_e(6), &
1289                    ros_alpha(6), ros_gamma(6), ros_elo
1290   LOGICAL :: ros_newf(6)
1291   CHARACTER(len=12):: ros_name
1292!~~~>  local variables
1293   REAL(kind=dp):: roundoff, facmin, facmax, facrej, facsafe
1294   REAL(kind=dp):: hmin, hmax, hstart
1295   REAL(kind=dp):: texit
1296   INTEGER       :: i, uplimtol, max_no_steps
1297   LOGICAL       :: autonomous, vectortol
1298!~~~>   PARAMETERs
1299   REAL(kind=dp), PARAMETER :: zero = 0.0_dp, one  = 1.0_dp
1300   REAL(kind=dp), PARAMETER :: deltamin = 1.0e-5_dp
1301
1302!~~~>  initialize statistics
1303   istatus(1:8) = 0
1304   rstatus(1:3) = zero
1305
1306!~~~>  autonomous or time dependent ode. default is time dependent.
1307   autonomous = .not.(icntrl(1) == 0)
1308
1309!~~~>  for scalar tolerances (icntrl(2).ne.0) the code uses abstol(1)and reltol(1)
1310!   For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:N) and RelTol(1:N)
1311   IF (icntrl(2) == 0)THEN
1312      vectortol = .TRUE.
1313      uplimtol  = n
1314   ELSE
1315      vectortol = .FALSE.
1316      uplimtol  = 1
1317   ENDIF
1318
1319!~~~>   initialize the particular rosenbrock method selected
1320   select CASE (icntrl(3))
1321     CASE (1)
1322       CALL ros2
1323     CASE (2)
1324       CALL ros3
1325     CASE (3)
1326       CALL ros4
1327     CASE (0, 4)
1328       CALL rodas3
1329     CASE (5)
1330       CALL rodas4
1331     CASE (6)
1332       CALL rang3
1333     CASE default
1334       PRINT *,'Unknown Rosenbrock method: ICNTRL(3) =',ICNTRL(3) 
1335       CALL ros_errormsg(- 2, tstart, zero, ierr)
1336       RETURN
1337   END select
1338
1339!~~~>   the maximum number of steps admitted
1340   IF (icntrl(4) == 0)THEN
1341      max_no_steps = 200000
1342   ELSEIF (icntrl(4)> 0)THEN
1343      max_no_steps=icntrl(4)
1344   ELSE
1345      PRINT *,'User-selected max no. of steps: ICNTRL(4) =',ICNTRL(4)
1346      CALL ros_errormsg(- 1, tstart, zero, ierr)
1347      RETURN
1348   ENDIF
1349
1350!~~~>  unit roundoff (1+ roundoff>1)
1351   roundoff = epsilon(one)
1352
1353!~~~>  lower bound on the step size: (positive value)
1354   IF (rcntrl(1) == zero)THEN
1355      hmin = zero
1356   ELSEIF (rcntrl(1)> zero)THEN
1357      hmin = rcntrl(1)
1358   ELSE
1359      PRINT *,'User-selected Hmin: RCNTRL(1) =',RCNTRL(1)
1360      CALL ros_errormsg(- 3, tstart, zero, ierr)
1361      RETURN
1362   ENDIF
1363!~~~>  upper bound on the step size: (positive value)
1364   IF (rcntrl(2) == zero)THEN
1365      hmax = abs(tend-tstart)
1366   ELSEIF (rcntrl(2)> zero)THEN
1367      hmax = min(abs(rcntrl(2)), abs(tend-tstart))
1368   ELSE
1369      PRINT *,'User-selected Hmax: RCNTRL(2) =',RCNTRL(2)
1370      CALL ros_errormsg(- 3, tstart, zero, ierr)
1371      RETURN
1372   ENDIF
1373!~~~>  starting step size: (positive value)
1374   IF (rcntrl(3) == zero)THEN
1375      hstart = max(hmin, deltamin)
1376   ELSEIF (rcntrl(3)> zero)THEN
1377      hstart = min(abs(rcntrl(3)), abs(tend-tstart))
1378   ELSE
1379      PRINT *,'User-selected Hstart: RCNTRL(3) =',RCNTRL(3)
1380      CALL ros_errormsg(- 3, tstart, zero, ierr)
1381      RETURN
1382   ENDIF
1383!~~~>  step size can be changed s.t.  facmin < hnew/hold < facmax
1384   IF (rcntrl(4) == zero)THEN
1385      facmin = 0.2_dp
1386   ELSEIF (rcntrl(4)> zero)THEN
1387      facmin = rcntrl(4)
1388   ELSE
1389      PRINT *,'User-selected FacMin: RCNTRL(4) =',RCNTRL(4)
1390      CALL ros_errormsg(- 4, tstart, zero, ierr)
1391      RETURN
1392   ENDIF
1393   IF (rcntrl(5) == zero)THEN
1394      facmax = 6.0_dp
1395   ELSEIF (rcntrl(5)> zero)THEN
1396      facmax = rcntrl(5)
1397   ELSE
1398      PRINT *,'User-selected FacMax: RCNTRL(5) =',RCNTRL(5)
1399      CALL ros_errormsg(- 4, tstart, zero, ierr)
1400      RETURN
1401   ENDIF
1402!~~~>   facrej: factor to decrease step after 2 succesive rejections
1403   IF (rcntrl(6) == zero)THEN
1404      facrej = 0.1_dp
1405   ELSEIF (rcntrl(6)> zero)THEN
1406      facrej = rcntrl(6)
1407   ELSE
1408      PRINT *,'User-selected FacRej: RCNTRL(6) =',RCNTRL(6)
1409      CALL ros_errormsg(- 4, tstart, zero, ierr)
1410      RETURN
1411   ENDIF
1412!~~~>   facsafe: safety factor in the computation of new step size
1413   IF (rcntrl(7) == zero)THEN
1414      facsafe = 0.9_dp
1415   ELSEIF (rcntrl(7)> zero)THEN
1416      facsafe = rcntrl(7)
1417   ELSE
1418      PRINT *,'User-selected FacSafe: RCNTRL(7) =',RCNTRL(7)
1419      CALL ros_errormsg(- 4, tstart, zero, ierr)
1420      RETURN
1421   ENDIF
1422!~~~>  check IF tolerances are reasonable
1423    DO i=1, uplimtol
1424      IF ((abstol(i)<= zero).or. (reltol(i)<= 10.0_dp* roundoff)&
1425         .or. (reltol(i)>= 1.0_dp))THEN
1426        PRINT *,' AbsTol(',i,') = ',AbsTol(i)
1427        PRINT *,' RelTol(',i,') = ',RelTol(i)
1428        CALL ros_errormsg(- 5, tstart, zero, ierr)
1429        RETURN
1430      ENDIF
1431    ENDDO
1432
1433
1434!~~~>  CALL rosenbrock method
1435   CALL ros_integrator(y, tstart, tend, texit,  &
1436        abstol, reltol,                         &
1437!  Integration parameters
1438        autonomous, vectortol, max_no_steps,    &
1439        roundoff, hmin, hmax, hstart,           &
1440        facmin, facmax, facrej, facsafe,        &
1441!  Error indicator
1442        ierr)
1443
1444!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1445CONTAINS !  SUBROUTINEs internal to rosenbrock
1446!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1447   
1448!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
1449 SUBROUTINE ros_errormsg(code, t, h, ierr)
1450!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
1451!    Handles all error messages
1452!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
1453 
1454   REAL(kind=dp), INTENT(IN):: t, h
1455   INTEGER, INTENT(IN) :: code
1456   INTEGER, INTENT(OUT):: ierr
1457   
1458   ierr = code
1459   print * , &
1460     'Forced exit from Rosenbrock due to the following error:' 
1461     
1462   select CASE (code)
1463    CASE (- 1) 
1464      PRINT *,'--> Improper value for maximal no of steps'
1465    CASE (- 2) 
1466      PRINT *,'--> Selected Rosenbrock method not implemented'
1467    CASE (- 3) 
1468      PRINT *,'--> Hmin/Hmax/Hstart must be positive'
1469    CASE (- 4) 
1470      PRINT *,'--> FacMin/FacMax/FacRej must be positive'
1471    CASE (- 5)
1472      PRINT *,'--> Improper tolerance values'
1473    CASE (- 6)
1474      PRINT *,'--> No of steps exceeds maximum bound'
1475    CASE (- 7)
1476      PRINT *,'--> Step size too small: T + 10*H = T',&
1477            ' or H < Roundoff'
1478    CASE (- 8) 
1479      PRINT *,'--> Matrix is repeatedly singular'
1480    CASE default
1481      PRINT *,'Unknown Error code: ',Code
1482   END select
1483   
1484   print * , "t=", t, "and h=", h
1485     
1486 END SUBROUTINE ros_errormsg
1487
1488!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1489 SUBROUTINE ros_integrator (y, tstart, tend, t, &
1490        abstol, reltol,                         &
1491!~~~> integration PARAMETERs
1492        autonomous, vectortol, max_no_steps,    &
1493        roundoff, hmin, hmax, hstart,           &
1494        facmin, facmax, facrej, facsafe,        &
1495!~~~> error indicator
1496        ierr)
1497!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1498!   Template for the implementation of a generic Rosenbrock method
1499!      defined by ros_S (no of stages)
1500!      and its coefficients ros_{A,C,M,E,Alpha,Gamma}
1501!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1502
1503
1504!~~~> input: the initial condition at tstart; output: the solution at t
1505   REAL(kind=dp), INTENT(INOUT):: y(n)
1506!~~~> input: integration interval
1507   REAL(kind=dp), INTENT(IN):: tstart, tend
1508!~~~> output: time at which the solution is RETURNed (t=tendIF success)
1509   REAL(kind=dp), INTENT(OUT)::  t
1510!~~~> input: tolerances
1511   REAL(kind=dp), INTENT(IN)::  abstol(n), reltol(n)
1512!~~~> input: integration PARAMETERs
1513   LOGICAL, INTENT(IN):: autonomous, vectortol
1514   REAL(kind=dp), INTENT(IN):: hstart, hmin, hmax
1515   INTEGER, INTENT(IN):: max_no_steps
1516   REAL(kind=dp), INTENT(IN):: roundoff, facmin, facmax, facrej, facsafe
1517!~~~> output: error indicator
1518   INTEGER, INTENT(OUT):: ierr
1519! ~~~~ Local variables
1520   REAL(kind=dp):: ynew(n), fcn0(n), fcn(n)
1521   REAL(kind=dp):: k(n* ros_s), dfdt(n)
1522#ifdef full_algebra   
1523   REAL(kind=dp):: jac0(n, n), ghimj(n, n)
1524#else
1525   REAL(kind=dp):: jac0(lu_nonzero), ghimj(lu_nonzero)
1526#endif
1527   REAL(kind=dp):: h, hnew, hc, hg, fac, tau
1528   REAL(kind=dp):: err, yerr(n)
1529   INTEGER :: pivot(n), direction, ioffset, j, istage
1530   LOGICAL :: rejectlasth, rejectmoreh, singular
1531!~~~>  local PARAMETERs
1532   REAL(kind=dp), PARAMETER :: zero = 0.0_dp, one  = 1.0_dp
1533   REAL(kind=dp), PARAMETER :: deltamin = 1.0e-5_dp
1534!~~~>  locally called FUNCTIONs
1535!    REAL(kind=dp) WLAMCH
1536!    EXTERNAL WLAMCH
1537!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1538
1539
1540!~~~>  initial preparations
1541   t = tstart
1542   rstatus(nhexit) = zero
1543   h = min( max(abs(hmin), abs(hstart)), abs(hmax))
1544   IF (abs(h)<= 10.0_dp* roundoff)h = deltamin
1545
1546   IF (tend  >=  tstart)THEN
1547     direction = + 1
1548   ELSE
1549     direction = - 1
1550   ENDIF
1551   h = direction* h
1552
1553   rejectlasth=.FALSE.
1554   rejectmoreh=.FALSE.
1555
1556!~~~> time loop begins below
1557
1558timeloop: DO WHILE((direction > 0).and.((t- tend) + roundoff <= zero)&
1559       .or. (direction < 0).and.((tend-t) + roundoff <= zero))
1560
1561   IF (istatus(nstp)> max_no_steps)THEN  ! too many steps
1562      CALL ros_errormsg(- 6, t, h, ierr)
1563      RETURN
1564   ENDIF
1565   IF (((t+ 0.1_dp* h) == t).or.(h <= roundoff))THEN  ! step size too small
1566      CALL ros_errormsg(- 7, t, h, ierr)
1567      RETURN
1568   ENDIF
1569
1570!~~~>  limit h IF necessary to avoid going beyond tend
1571   h = min(h, abs(tend-t))
1572
1573!~~~>   compute the FUNCTION at current time
1574   CALL funtemplate(t, y, fcn0)
1575   istatus(nfun) = istatus(nfun) + 1
1576
1577!~~~>  compute the FUNCTION derivative with respect to t
1578   IF (.not.autonomous)THEN
1579      CALL ros_funtimederivative(t, roundoff, y, &
1580                fcn0, dfdt)
1581   ENDIF
1582
1583!~~~>   compute the jacobian at current time
1584   CALL jactemplate(t, y, jac0)
1585   istatus(njac) = istatus(njac) + 1
1586
1587!~~~>  repeat step calculation until current step accepted
1588untilaccepted: do
1589
1590   CALL ros_preparematrix(h, direction, ros_gamma(1), &
1591          jac0, ghimj, pivot, singular)
1592   IF (singular)THEN ! more than 5 consecutive failed decompositions
1593       CALL ros_errormsg(- 8, t, h, ierr)
1594       RETURN
1595   ENDIF
1596
1597!~~~>   compute the stages
1598stage: DO istage = 1, ros_s
1599
1600      ! current istage offset. current istage vector is k(ioffset+ 1:ioffset+ n)
1601       ioffset = n* (istage-1)
1602
1603      ! for the 1st istage the FUNCTION has been computed previously
1604       IF (istage == 1)THEN
1605         !slim: CALL wcopy(n, fcn0, 1, fcn, 1)
1606       fcn(1:n) = fcn0(1:n)
1607      ! istage>1 and a new FUNCTION evaluation is needed at the current istage
1608       ELSEIF(ros_newf(istage))THEN
1609         !slim: CALL wcopy(n, y, 1, ynew, 1)
1610       ynew(1:n) = y(1:n)
1611         DO j = 1, istage-1
1612           CALL waxpy(n, ros_a((istage-1) * (istage-2) /2+ j), &
1613            k(n* (j- 1) + 1), 1, ynew, 1)
1614         ENDDO
1615         tau = t + ros_alpha(istage) * direction* h
1616         CALL funtemplate(tau, ynew, fcn)
1617         istatus(nfun) = istatus(nfun) + 1
1618       ENDIF ! IF istage == 1 ELSEIF ros_newf(istage)
1619       !slim: CALL wcopy(n, fcn, 1, k(ioffset+ 1), 1)
1620       k(ioffset+ 1:ioffset+ n) = fcn(1:n)
1621       DO j = 1, istage-1
1622         hc = ros_c((istage-1) * (istage-2) /2+ j) /(direction* h)
1623         CALL waxpy(n, hc, k(n* (j- 1) + 1), 1, k(ioffset+ 1), 1)
1624       ENDDO
1625       IF ((.not. autonomous).and.(ros_gamma(istage).ne.zero))THEN
1626         hg = direction* h* ros_gamma(istage)
1627         CALL waxpy(n, hg, dfdt, 1, k(ioffset+ 1), 1)
1628       ENDIF
1629       CALL ros_solve(ghimj, pivot, k(ioffset+ 1))
1630
1631   END DO stage
1632
1633
1634!~~~>  compute the new solution
1635   !slim: CALL wcopy(n, y, 1, ynew, 1)
1636   ynew(1:n) = y(1:n)
1637   DO j=1, ros_s
1638         CALL waxpy(n, ros_m(j), k(n* (j- 1) + 1), 1, ynew, 1)
1639   ENDDO
1640
1641!~~~>  compute the error estimation
1642   !slim: CALL wscal(n, zero, yerr, 1)
1643   yerr(1:n) = zero
1644   DO j=1, ros_s
1645        CALL waxpy(n, ros_e(j), k(n* (j- 1) + 1), 1, yerr, 1)
1646   ENDDO
1647   err = ros_errornorm(y, ynew, yerr, abstol, reltol, vectortol)
1648
1649!~~~> new step size is bounded by facmin <= hnew/h <= facmax
1650   fac  = min(facmax, max(facmin, facsafe/err** (one/ros_elo)))
1651   hnew = h* fac
1652
1653!~~~>  check the error magnitude and adjust step size
1654   istatus(nstp) = istatus(nstp) + 1
1655   IF ((err <= one).or.(h <= hmin))THEN  !~~~> accept step
1656      istatus(nacc) = istatus(nacc) + 1
1657      !slim: CALL wcopy(n, ynew, 1, y, 1)
1658      y(1:n) = ynew(1:n)
1659      t = t + direction* h
1660      hnew = max(hmin, min(hnew, hmax))
1661      IF (rejectlasth)THEN  ! no step size increase after a rejected step
1662         hnew = min(hnew, h)
1663      ENDIF
1664      rstatus(nhexit) = h
1665      rstatus(nhnew) = hnew
1666      rstatus(ntexit) = t
1667      rejectlasth = .FALSE.
1668      rejectmoreh = .FALSE.
1669      h = hnew
1670      exit untilaccepted ! exit the loop: WHILE step not accepted
1671   ELSE           !~~~> reject step
1672      IF (rejectmoreh)THEN
1673         hnew = h* facrej
1674      ENDIF
1675      rejectmoreh = rejectlasth
1676      rejectlasth = .TRUE.
1677      h = hnew
1678      IF (istatus(nacc)>= 1) istatus(nrej) = istatus(nrej) + 1
1679   ENDIF ! err <= 1
1680
1681   END DO untilaccepted
1682
1683   END DO timeloop
1684
1685!~~~> succesful exit
1686   ierr = 1  !~~~> the integration was successful
1687
1688  END SUBROUTINE ros_integrator
1689
1690
1691!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1692  REAL(kind=dp)FUNCTION ros_errornorm(y, ynew, yerr, &
1693                               abstol, reltol, vectortol)
1694!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1695!~~~> computes the "scaled norm" of the error vector yerr
1696!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1697
1698! Input arguments
1699   REAL(kind=dp), INTENT(IN):: y(n), ynew(n), &
1700          yerr(n), abstol(n), reltol(n)
1701   LOGICAL, INTENT(IN)::  vectortol
1702! Local variables
1703   REAL(kind=dp):: err, scale, ymax
1704   INTEGER  :: i
1705   REAL(kind=dp), PARAMETER :: zero = 0.0_dp
1706
1707   err = zero
1708   DO i=1, n
1709     ymax = max(abs(y(i)), abs(ynew(i)))
1710     IF (vectortol)THEN
1711       scale = abstol(i) + reltol(i) * ymax
1712     ELSE
1713       scale = abstol(1) + reltol(1) * ymax
1714     ENDIF
1715     err = err+ (yerr(i) /scale) ** 2
1716   ENDDO
1717   err  = sqrt(err/n)
1718
1719   ros_errornorm = max(err, 1.0d-10)
1720
1721  END FUNCTION ros_errornorm
1722
1723
1724!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1725  SUBROUTINE ros_funtimederivative(t, roundoff, y, &
1726                fcn0, dfdt)
1727!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1728!~~~> the time partial derivative of the FUNCTION by finite differences
1729!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1730
1731!~~~> input arguments
1732   REAL(kind=dp), INTENT(IN):: t, roundoff, y(n), fcn0(n)
1733!~~~> output arguments
1734   REAL(kind=dp), INTENT(OUT):: dfdt(n)
1735!~~~> local variables
1736   REAL(kind=dp):: delta
1737   REAL(kind=dp), PARAMETER :: one = 1.0_dp, deltamin = 1.0e-6_dp
1738
1739   delta = sqrt(roundoff) * max(deltamin, abs(t))
1740   CALL funtemplate(t+ delta, y, dfdt)
1741   istatus(nfun) = istatus(nfun) + 1
1742   CALL waxpy(n, (- one), fcn0, 1, dfdt, 1)
1743   CALL wscal(n, (one/delta), dfdt, 1)
1744
1745  END SUBROUTINE ros_funtimederivative
1746
1747
1748!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1749  SUBROUTINE ros_preparematrix(h, direction, gam, &
1750             jac0, ghimj, pivot, singular)
1751! --- --- --- --- --- --- --- --- --- --- --- --- ---
1752!  Prepares the LHS matrix for stage calculations
1753!  1.  Construct Ghimj = 1/(H*ham) - Jac0
1754!      "(Gamma H) Inverse Minus Jacobian"
1755!  2.  Repeat LU decomposition of Ghimj until successful.
1756!       -half the step size if LU decomposition fails and retry
1757!       -exit after 5 consecutive fails
1758! --- --- --- --- --- --- --- --- --- --- --- --- ---
1759
1760!~~~> input arguments
1761#ifdef full_algebra   
1762   REAL(kind=dp), INTENT(IN)::  jac0(n, n)
1763#else
1764   REAL(kind=dp), INTENT(IN)::  jac0(lu_nonzero)
1765#endif   
1766   REAL(kind=dp), INTENT(IN)::  gam
1767   INTEGER, INTENT(IN)::  direction
1768!~~~> output arguments
1769#ifdef full_algebra   
1770   REAL(kind=dp), INTENT(OUT):: ghimj(n, n)
1771#else
1772   REAL(kind=dp), INTENT(OUT):: ghimj(lu_nonzero)
1773#endif   
1774   LOGICAL, INTENT(OUT)::  singular
1775   INTEGER, INTENT(OUT)::  pivot(n)
1776!~~~> inout arguments
1777   REAL(kind=dp), INTENT(INOUT):: h   ! step size is decreased when lu fails
1778!~~~> local variables
1779   INTEGER  :: i, ising, nconsecutive
1780   REAL(kind=dp):: ghinv
1781   REAL(kind=dp), PARAMETER :: one  = 1.0_dp, half = 0.5_dp
1782
1783   nconsecutive = 0
1784   singular = .TRUE.
1785
1786   DO WHILE (singular)
1787
1788!~~~>    construct ghimj = 1/(h* gam) - jac0
1789#ifdef full_algebra   
1790     !slim: CALL wcopy(n* n, jac0, 1, ghimj, 1)
1791     !slim: CALL wscal(n* n, (- one), ghimj, 1)
1792     ghimj = - jac0
1793     ghinv = one/(direction* h* gam)
1794     DO i=1, n
1795       ghimj(i, i) = ghimj(i, i) + ghinv
1796     ENDDO
1797#else
1798     !slim: CALL wcopy(lu_nonzero, jac0, 1, ghimj, 1)
1799     !slim: CALL wscal(lu_nonzero, (- one), ghimj, 1)
1800     ghimj(1:lu_nonzero) = - jac0(1:lu_nonzero)
1801     ghinv = one/(direction* h* gam)
1802     DO i=1, n
1803       ghimj(lu_diag(i)) = ghimj(lu_diag(i)) + ghinv
1804     ENDDO
1805#endif   
1806!~~~>    compute lu decomposition
1807     CALL ros_decomp( ghimj, pivot, ising)
1808     IF (ising == 0)THEN
1809!~~~>    IF successful done
1810        singular = .FALSE.
1811     ELSE ! ising .ne. 0
1812!~~~>    IF unsuccessful half the step size; IF 5 consecutive fails THEN RETURN
1813        istatus(nsng) = istatus(nsng) + 1
1814        nconsecutive = nconsecutive+1
1815        singular = .TRUE.
1816        PRINT*,'Warning: LU Decomposition returned ISING = ',ISING
1817        IF (nconsecutive <= 5)THEN ! less than 5 consecutive failed decompositions
1818           h = h* half
1819        ELSE  ! more than 5 consecutive failed decompositions
1820           RETURN
1821        ENDIF  ! nconsecutive
1822      ENDIF    ! ising
1823
1824   END DO ! WHILE singular
1825
1826  END SUBROUTINE ros_preparematrix
1827
1828
1829!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1830  SUBROUTINE ros_decomp( a, pivot, ising)
1831!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1832!  Template for the LU decomposition
1833!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1834!~~~> inout variables
1835#ifdef full_algebra   
1836   REAL(kind=dp), INTENT(INOUT):: a(n, n)
1837#else   
1838   REAL(kind=dp), INTENT(INOUT):: a(lu_nonzero)
1839#endif
1840!~~~> output variables
1841   INTEGER, INTENT(OUT):: pivot(n), ising
1842
1843#ifdef full_algebra   
1844   CALL  dgetrf( n, n, a, n, pivot, ising)
1845#else   
1846   CALL kppdecomp(a, ising)
1847   pivot(1) = 1
1848#endif
1849   istatus(ndec) = istatus(ndec) + 1
1850
1851  END SUBROUTINE ros_decomp
1852
1853
1854!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1855  SUBROUTINE ros_solve( a, pivot, b)
1856!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1857!  Template for the forward/backward substitution (using pre-computed LU decomposition)
1858!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1859!~~~> input variables
1860#ifdef full_algebra   
1861   REAL(kind=dp), INTENT(IN):: a(n, n)
1862   INTEGER :: ising
1863#else   
1864   REAL(kind=dp), INTENT(IN):: a(lu_nonzero)
1865#endif
1866   INTEGER, INTENT(IN):: pivot(n)
1867!~~~> inout variables
1868   REAL(kind=dp), INTENT(INOUT):: b(n)
1869
1870#ifdef full_algebra   
1871   CALL  DGETRS( 'N',N ,1,A,N,Pivot,b,N,ISING)
1872   IF (info < 0)THEN
1873      print* , "error in dgetrs. ising=", ising
1874   ENDIF 
1875#else   
1876   CALL kppsolve( a, b)
1877#endif
1878
1879   istatus(nsol) = istatus(nsol) + 1
1880
1881  END SUBROUTINE ros_solve
1882
1883
1884
1885!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1886  SUBROUTINE ros2
1887!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1888! --- AN L-STABLE METHOD,2 stages,order 2
1889!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1890
1891   double precision g
1892
1893    g = 1.0_dp + 1.0_dp/sqrt(2.0_dp)
1894    rosmethod = rs2
1895!~~~> name of the method
1896    ros_Name = 'ROS-2'
1897!~~~> number of stages
1898    ros_s = 2
1899
1900!~~~> the coefficient matrices a and c are strictly lower triangular.
1901!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1902!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
1903!   The general mapping formula is:
1904!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
1905!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
1906
1907    ros_a(1) = (1.0_dp) /g
1908    ros_c(1) = (- 2.0_dp) /g
1909!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
1910!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
1911    ros_newf(1) = .TRUE.
1912    ros_newf(2) = .TRUE.
1913!~~~> m_i = coefficients for new step solution
1914    ros_m(1) = (3.0_dp) /(2.0_dp* g)
1915    ros_m(2) = (1.0_dp) /(2.0_dp* g)
1916! E_i = Coefficients for error estimator
1917    ros_e(1) = 1.0_dp/(2.0_dp* g)
1918    ros_e(2) = 1.0_dp/(2.0_dp* g)
1919!~~~> ros_elo = estimator of local order - the minimum between the
1920!    main and the embedded scheme orders plus one
1921    ros_elo = 2.0_dp
1922!~~~> y_stage_i ~ y( t + h* alpha_i)
1923    ros_alpha(1) = 0.0_dp
1924    ros_alpha(2) = 1.0_dp
1925!~~~> gamma_i = \sum_j  gamma_{i, j}
1926    ros_gamma(1) = g
1927    ros_gamma(2) = -g
1928
1929 END SUBROUTINE ros2
1930
1931
1932!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1933  SUBROUTINE ros3
1934!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1935! --- AN L-STABLE METHOD,3 stages,order 3,2 function evaluations
1936!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1937
1938   rosmethod = rs3
1939!~~~> name of the method
1940   ros_Name = 'ROS-3'
1941!~~~> number of stages
1942   ros_s = 3
1943
1944!~~~> the coefficient matrices a and c are strictly lower triangular.
1945!   The lower triangular (subdiagonal) elements are stored in row-wise order:
1946!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
1947!   The general mapping formula is:
1948!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
1949!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
1950
1951   ros_a(1) = 1.0_dp
1952   ros_a(2) = 1.0_dp
1953   ros_a(3) = 0.0_dp
1954
1955   ros_c(1) = - 0.10156171083877702091975600115545e+01_dp
1956   ros_c(2) =  0.40759956452537699824805835358067e+01_dp
1957   ros_c(3) =  0.92076794298330791242156818474003e+01_dp
1958!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
1959!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
1960   ros_newf(1) = .TRUE.
1961   ros_newf(2) = .TRUE.
1962   ros_newf(3) = .FALSE.
1963!~~~> m_i = coefficients for new step solution
1964   ros_m(1) =  0.1e+01_dp
1965   ros_m(2) =  0.61697947043828245592553615689730e+01_dp
1966   ros_m(3) = - 0.42772256543218573326238373806514_dp
1967! E_i = Coefficients for error estimator
1968   ros_e(1) =  0.5_dp
1969   ros_e(2) = - 0.29079558716805469821718236208017e+01_dp
1970   ros_e(3) =  0.22354069897811569627360909276199_dp
1971!~~~> ros_elo = estimator of local order - the minimum between the
1972!    main and the embedded scheme orders plus 1
1973   ros_elo = 3.0_dp
1974!~~~> y_stage_i ~ y( t + h* alpha_i)
1975   ros_alpha(1) = 0.0_dp
1976   ros_alpha(2) = 0.43586652150845899941601945119356_dp
1977   ros_alpha(3) = 0.43586652150845899941601945119356_dp
1978!~~~> gamma_i = \sum_j  gamma_{i, j}
1979   ros_gamma(1) = 0.43586652150845899941601945119356_dp
1980   ros_gamma(2) = 0.24291996454816804366592249683314_dp
1981   ros_gamma(3) = 0.21851380027664058511513169485832e+01_dp
1982
1983  END SUBROUTINE ros3
1984
1985!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1986
1987
1988!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1989  SUBROUTINE ros4
1990!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1991!     L-STABLE ROSENBROCK METHOD OF ORDER 4,WITH 4 STAGES
1992!     L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
1993!
1994!      E. HAIRER AND G. WANNER,SOLVING ORDINARY DIFFERENTIAL
1995!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
1996!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
1997!      SPRINGER-VERLAG (1990)
1998!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1999
2000
2001   rosmethod = rs4
2002!~~~> name of the method
2003   ros_Name = 'ROS-4'
2004!~~~> number of stages
2005   ros_s = 4
2006
2007!~~~> the coefficient matrices a and c are strictly lower triangular.
2008!   The lower triangular (subdiagonal) elements are stored in row-wise order:
2009!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
2010!   The general mapping formula is:
2011!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
2012!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
2013
2014   ros_a(1) = 0.2000000000000000e+01_dp
2015   ros_a(2) = 0.1867943637803922e+01_dp
2016   ros_a(3) = 0.2344449711399156_dp
2017   ros_a(4) = ros_a(2)
2018   ros_a(5) = ros_a(3)
2019   ros_a(6) = 0.0_dp
2020
2021   ros_c(1) = -0.7137615036412310e+01_dp
2022   ros_c(2) = 0.2580708087951457e+01_dp
2023   ros_c(3) = 0.6515950076447975_dp
2024   ros_c(4) = -0.2137148994382534e+01_dp
2025   ros_c(5) = -0.3214669691237626_dp
2026   ros_c(6) = -0.6949742501781779_dp
2027!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
2028!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
2029   ros_newf(1) = .TRUE.
2030   ros_newf(2) = .TRUE.
2031   ros_newf(3) = .TRUE.
2032   ros_newf(4) = .FALSE.
2033!~~~> m_i = coefficients for new step solution
2034   ros_m(1) = 0.2255570073418735e+01_dp
2035   ros_m(2) = 0.2870493262186792_dp
2036   ros_m(3) = 0.4353179431840180_dp
2037   ros_m(4) = 0.1093502252409163e+01_dp
2038!~~~> e_i  = coefficients for error estimator
2039   ros_e(1) = -0.2815431932141155_dp
2040   ros_e(2) = -0.7276199124938920e-01_dp
2041   ros_e(3) = -0.1082196201495311_dp
2042   ros_e(4) = -0.1093502252409163e+01_dp
2043!~~~> ros_elo  = estimator of local order - the minimum between the
2044!    main and the embedded scheme orders plus 1
2045   ros_elo  = 4.0_dp
2046!~~~> y_stage_i ~ y( t + h* alpha_i)
2047   ros_alpha(1) = 0.0_dp
2048   ros_alpha(2) = 0.1145640000000000e+01_dp
2049   ros_alpha(3) = 0.6552168638155900_dp
2050   ros_alpha(4) = ros_alpha(3)
2051!~~~> gamma_i = \sum_j  gamma_{i, j}
2052   ros_gamma(1) = 0.5728200000000000_dp
2053   ros_gamma(2) = -0.1769193891319233e+01_dp
2054   ros_gamma(3) = 0.7592633437920482_dp
2055   ros_gamma(4) = -0.1049021087100450_dp
2056
2057  END SUBROUTINE ros4
2058
2059!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2060  SUBROUTINE rodas3
2061!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2062! --- A STIFFLY-STABLE METHOD,4 stages,order 3
2063!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2064
2065
2066   rosmethod = rd3
2067!~~~> name of the method
2068   ros_Name = 'RODAS-3'
2069!~~~> number of stages
2070   ros_s = 4
2071
2072!~~~> the coefficient matrices a and c are strictly lower triangular.
2073!   The lower triangular (subdiagonal) elements are stored in row-wise order:
2074!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
2075!   The general mapping formula is:
2076!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
2077!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
2078
2079   ros_a(1) = 0.0_dp
2080   ros_a(2) = 2.0_dp
2081   ros_a(3) = 0.0_dp
2082   ros_a(4) = 2.0_dp
2083   ros_a(5) = 0.0_dp
2084   ros_a(6) = 1.0_dp
2085
2086   ros_c(1) = 4.0_dp
2087   ros_c(2) = 1.0_dp
2088   ros_c(3) = -1.0_dp
2089   ros_c(4) = 1.0_dp
2090   ros_c(5) = -1.0_dp
2091   ros_c(6) = -(8.0_dp/3.0_dp)
2092
2093!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
2094!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
2095   ros_newf(1) = .TRUE.
2096   ros_newf(2) = .FALSE.
2097   ros_newf(3) = .TRUE.
2098   ros_newf(4) = .TRUE.
2099!~~~> m_i = coefficients for new step solution
2100   ros_m(1) = 2.0_dp
2101   ros_m(2) = 0.0_dp
2102   ros_m(3) = 1.0_dp
2103   ros_m(4) = 1.0_dp
2104!~~~> e_i  = coefficients for error estimator
2105   ros_e(1) = 0.0_dp
2106   ros_e(2) = 0.0_dp
2107   ros_e(3) = 0.0_dp
2108   ros_e(4) = 1.0_dp
2109!~~~> ros_elo  = estimator of local order - the minimum between the
2110!    main and the embedded scheme orders plus 1
2111   ros_elo  = 3.0_dp
2112!~~~> y_stage_i ~ y( t + h* alpha_i)
2113   ros_alpha(1) = 0.0_dp
2114   ros_alpha(2) = 0.0_dp
2115   ros_alpha(3) = 1.0_dp
2116   ros_alpha(4) = 1.0_dp
2117!~~~> gamma_i = \sum_j  gamma_{i, j}
2118   ros_gamma(1) = 0.5_dp
2119   ros_gamma(2) = 1.5_dp
2120   ros_gamma(3) = 0.0_dp
2121   ros_gamma(4) = 0.0_dp
2122
2123  END SUBROUTINE rodas3
2124
2125!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2126  SUBROUTINE rodas4
2127!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2128!     STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4,WITH 6 STAGES
2129!
2130!      E. HAIRER AND G. WANNER,SOLVING ORDINARY DIFFERENTIAL
2131!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
2132!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
2133!      SPRINGER-VERLAG (1996)
2134!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2135
2136
2137    rosmethod = rd4
2138!~~~> name of the method
2139    ros_Name = 'RODAS-4'
2140!~~~> number of stages
2141    ros_s = 6
2142
2143!~~~> y_stage_i ~ y( t + h* alpha_i)
2144    ros_alpha(1) = 0.000_dp
2145    ros_alpha(2) = 0.386_dp
2146    ros_alpha(3) = 0.210_dp
2147    ros_alpha(4) = 0.630_dp
2148    ros_alpha(5) = 1.000_dp
2149    ros_alpha(6) = 1.000_dp
2150
2151!~~~> gamma_i = \sum_j  gamma_{i, j}
2152    ros_gamma(1) = 0.2500000000000000_dp
2153    ros_gamma(2) = -0.1043000000000000_dp
2154    ros_gamma(3) = 0.1035000000000000_dp
2155    ros_gamma(4) = -0.3620000000000023e-01_dp
2156    ros_gamma(5) = 0.0_dp
2157    ros_gamma(6) = 0.0_dp
2158
2159!~~~> the coefficient matrices a and c are strictly lower triangular.
2160!   The lower triangular (subdiagonal) elements are stored in row-wise order:
2161!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
2162!   The general mapping formula is:  A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
2163!                  C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
2164
2165    ros_a(1) = 0.1544000000000000e+01_dp
2166    ros_a(2) = 0.9466785280815826_dp
2167    ros_a(3) = 0.2557011698983284_dp
2168    ros_a(4) = 0.3314825187068521e+01_dp
2169    ros_a(5) = 0.2896124015972201e+01_dp
2170    ros_a(6) = 0.9986419139977817_dp
2171    ros_a(7) = 0.1221224509226641e+01_dp
2172    ros_a(8) = 0.6019134481288629e+01_dp
2173    ros_a(9) = 0.1253708332932087e+02_dp
2174    ros_a(10) = -0.6878860361058950_dp
2175    ros_a(11) = ros_a(7)
2176    ros_a(12) = ros_a(8)
2177    ros_a(13) = ros_a(9)
2178    ros_a(14) = ros_a(10)
2179    ros_a(15) = 1.0_dp
2180
2181    ros_c(1) = -0.5668800000000000e+01_dp
2182    ros_c(2) = -0.2430093356833875e+01_dp
2183    ros_c(3) = -0.2063599157091915_dp
2184    ros_c(4) = -0.1073529058151375_dp
2185    ros_c(5) = -0.9594562251023355e+01_dp
2186    ros_c(6) = -0.2047028614809616e+02_dp
2187    ros_c(7) = 0.7496443313967647e+01_dp
2188    ros_c(8) = -0.1024680431464352e+02_dp
2189    ros_c(9) = -0.3399990352819905e+02_dp
2190    ros_c(10) = 0.1170890893206160e+02_dp
2191    ros_c(11) = 0.8083246795921522e+01_dp
2192    ros_c(12) = -0.7981132988064893e+01_dp
2193    ros_c(13) = -0.3152159432874371e+02_dp
2194    ros_c(14) = 0.1631930543123136e+02_dp
2195    ros_c(15) = -0.6058818238834054e+01_dp
2196
2197!~~~> m_i = coefficients for new step solution
2198    ros_m(1) = ros_a(7)
2199    ros_m(2) = ros_a(8)
2200    ros_m(3) = ros_a(9)
2201    ros_m(4) = ros_a(10)
2202    ros_m(5) = 1.0_dp
2203    ros_m(6) = 1.0_dp
2204
2205!~~~> e_i  = coefficients for error estimator
2206    ros_e(1) = 0.0_dp
2207    ros_e(2) = 0.0_dp
2208    ros_e(3) = 0.0_dp
2209    ros_e(4) = 0.0_dp
2210    ros_e(5) = 0.0_dp
2211    ros_e(6) = 1.0_dp
2212
2213!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
2214!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
2215    ros_newf(1) = .TRUE.
2216    ros_newf(2) = .TRUE.
2217    ros_newf(3) = .TRUE.
2218    ros_newf(4) = .TRUE.
2219    ros_newf(5) = .TRUE.
2220    ros_newf(6) = .TRUE.
2221
2222!~~~> ros_elo  = estimator of local order - the minimum between the
2223!        main and the embedded scheme orders plus 1
2224    ros_elo = 4.0_dp
2225
2226  END SUBROUTINE rodas4
2227!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2228  SUBROUTINE rang3
2229!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2230! STIFFLY-STABLE W METHOD OF ORDER 3,WITH 4 STAGES
2231!
2232! J. RANG and L. ANGERMANN
2233! NEW ROSENBROCK W-METHODS OF ORDER 3
2234! FOR PARTIAL DIFFERENTIAL ALGEBRAIC
2235!        EQUATIONS OF INDEX 1                                             
2236! BIT Numerical Mathematics (2005) 45: 761-787
2237!  DOI: 10.1007/s10543-005-0035-y
2238! Table 4.1-4.2
2239!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2240
2241
2242    rosmethod = rg3
2243!~~~> name of the method
2244    ros_Name = 'RANG-3'
2245!~~~> number of stages
2246    ros_s = 4
2247
2248    ros_a(1) = 5.09052051067020d+00;
2249    ros_a(2) = 5.09052051067020d+00;
2250    ros_a(3) = 0.0d0;
2251    ros_a(4) = 4.97628111010787d+00;
2252    ros_a(5) = 2.77268164715849d-02;
2253    ros_a(6) = 2.29428036027904d-01;
2254
2255    ros_c(1) = - 1.16790812312283d+01;
2256    ros_c(2) = - 1.64057326467367d+01;
2257    ros_c(3) = - 2.77268164715850d-01;
2258    ros_c(4) = - 8.38103960500476d+00;
2259    ros_c(5) = - 8.48328409199343d-01;
2260    ros_c(6) =  2.87009860433106d-01;
2261
2262    ros_m(1) =  5.22582761233094d+00;
2263    ros_m(2) = - 5.56971148154165d-01;
2264    ros_m(3) =  3.57979469353645d-01;
2265    ros_m(4) =  1.72337398521064d+00;
2266
2267    ros_e(1) = - 5.16845212784040d+00;
2268    ros_e(2) = - 1.26351942603842d+00;
2269    ros_e(3) = - 1.11022302462516d-16;
2270    ros_e(4) =  2.22044604925031d-16;
2271
2272    ros_alpha(1) = 0.0d00;
2273    ros_alpha(2) = 2.21878746765329d+00;
2274    ros_alpha(3) = 2.21878746765329d+00;
2275    ros_alpha(4) = 1.55392337535788d+00;
2276
2277    ros_gamma(1) =  4.35866521508459d-01;
2278    ros_gamma(2) = - 1.78292094614483d+00;
2279    ros_gamma(3) = - 2.46541900496934d+00;
2280    ros_gamma(4) = - 8.05529997906370d-01;
2281
2282
2283!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
2284!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
2285    ros_newf(1) = .TRUE.
2286    ros_newf(2) = .TRUE.
2287    ros_newf(3) = .TRUE.
2288    ros_newf(4) = .TRUE.
2289
2290!~~~> ros_elo  = estimator of local order - the minimum between the
2291!        main and the embedded scheme orders plus 1
2292    ros_elo = 3.0_dp
2293
2294  END SUBROUTINE rang3
2295!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2296
2297!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2298!   End of the set of internal Rosenbrock subroutines
2299!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2300END SUBROUTINE rosenbrock
2301 
2302SUBROUTINE funtemplate( t, y, ydot)
2303!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2304!  Template for the ODE function call.
2305!  Updates the rate coefficients (and possibly the fixed species) at each call
2306!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2307!~~~> input variables
2308   REAL(kind=dp):: t, y(nvar)
2309!~~~> output variables
2310   REAL(kind=dp):: ydot(nvar)
2311!~~~> local variables
2312   REAL(kind=dp):: told
2313
2314   told = time
2315   time = t
2316   CALL fun( y, fix, rconst, ydot)
2317   time = told
2318
2319END SUBROUTINE funtemplate
2320 
2321SUBROUTINE jactemplate( t, y, jcb)
2322!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2323!  Template for the ODE Jacobian call.
2324!  Updates the rate coefficients (and possibly the fixed species) at each call
2325!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2326!~~~> input variables
2327    REAL(kind=dp):: t, y(nvar)
2328!~~~> output variables
2329#ifdef full_algebra   
2330    REAL(kind=dp):: jv(lu_nonzero), jcb(nvar, nvar)
2331#else
2332    REAL(kind=dp):: jcb(lu_nonzero)
2333#endif   
2334!~~~> local variables
2335    REAL(kind=dp):: told
2336#ifdef full_algebra   
2337    INTEGER :: i, j
2338#endif   
2339
2340    told = time
2341    time = t
2342#ifdef full_algebra   
2343    CALL jac_sp(y, fix, rconst, jv)
2344    DO j=1, nvar
2345      DO i=1, nvar
2346         jcb(i, j) = 0.0_dp
2347      ENDDO
2348    ENDDO
2349    DO i=1, lu_nonzero
2350       jcb(lu_irow(i), lu_icol(i)) = jv(i)
2351    ENDDO
2352#else
2353    CALL jac_sp( y, fix, rconst, jcb)
2354#endif   
2355    time = told
2356
2357END SUBROUTINE jactemplate
2358 
2359  SUBROUTINE kppdecomp( jvs, ier)                                   
2360  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2361  !        sparse lu factorization                                   
2362  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2363  ! loop expansion generated by kp4                                   
2364                                                                     
2365    INTEGER  :: ier                                                   
2366    REAL(kind=dp):: jvs(lu_nonzero), w(nvar), a             
2367    INTEGER  :: k, kk, j, jj                                         
2368                                                                     
2369    a = 0.                                                           
2370    ier = 0                                                           
2371                                                                     
2372!   i = 1
2373!   i = 2
2374!   i = 3
2375!   i = 4
2376!   i = 5
2377!   i = 6
2378!   i = 7
2379    jvs(16) = (jvs(16)) / jvs(8) 
2380    jvs(17) = (jvs(17)) / jvs(10) 
2381    jvs(18) = (jvs(18)) / jvs(13) 
2382    jvs(19) = jvs(19) - jvs(9) * jvs(16)
2383    jvs(20) = jvs(20) - jvs(11) * jvs(17)
2384    jvs(21) = jvs(21) - jvs(12) * jvs(17) - jvs(14) * jvs(18)
2385    jvs(22) = jvs(22) - jvs(15) * jvs(18)
2386!   i = 8
2387    jvs(23) = (jvs(23)) / jvs(8) 
2388    a = 0.0; a = a  - jvs(9) * jvs(23)
2389    jvs(24) = (jvs(24) + a) / jvs(19) 
2390    jvs(25) = jvs(25) - jvs(20) * jvs(24)
2391    jvs(26) = jvs(26) - jvs(21) * jvs(24)
2392    jvs(27) = jvs(27) - jvs(22) * jvs(24)
2393!   i = 9
2394    jvs(28) = (jvs(28)) / jvs(10) 
2395    jvs(29) = (jvs(29)) / jvs(13) 
2396    a = 0.0; a = a  - jvs(11) * jvs(28)
2397    jvs(30) = (jvs(30) + a) / jvs(25) 
2398    jvs(31) = jvs(31) - jvs(12) * jvs(28) - jvs(14) * jvs(29) - jvs(26) * jvs(30)
2399    jvs(32) = jvs(32) - jvs(15) * jvs(29) - jvs(27) * jvs(30)
2400!   i = 10
2401    jvs(33) = (jvs(33)) / jvs(10) 
2402    jvs(34) = (jvs(34)) / jvs(13) 
2403    jvs(35) = (jvs(35)) / jvs(19) 
2404    a = 0.0; a = a  - jvs(11) * jvs(33) - jvs(20) * jvs(35)
2405    jvs(36) = (jvs(36) + a) / jvs(25) 
2406    a = 0.0; a = a  - jvs(12) * jvs(33) - jvs(14) * jvs(34) - jvs(21) * jvs(35) - jvs(26) * jvs(36)
2407    jvs(37) = (jvs(37) + a) / jvs(31) 
2408    jvs(38) = jvs(38) - jvs(15) * jvs(34) - jvs(22) * jvs(35) - jvs(27) * jvs(36) - jvs(32) * jvs(37)
2409    RETURN                                                           
2410                                                                     
2411  END SUBROUTINE kppdecomp                                           
2412 
2413SUBROUTINE chem_gasphase_integrate (time_step_len, conc, tempi, qvapi, fakti, photo, ierrf, xnacc, xnrej, istatus, l_debug, pe, &
2414                     icntrl_i, rcntrl_i) 
2415                                                                   
2416  IMPLICIT NONE                                                     
2417                                                                   
2418  REAL(dp), INTENT(IN)                  :: time_step_len           
2419  REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: conc                   
2420  REAL(dp), DIMENSION(:, :), INTENT(IN)   :: photo                   
2421  REAL(dp), DIMENSION(:), INTENT(IN)     :: tempi                   
2422  REAL(dp), DIMENSION(:), INTENT(IN)     :: qvapi                   
2423  REAL(dp), DIMENSION(:), INTENT(IN)     :: fakti                   
2424  INTEGER, INTENT(OUT), OPTIONAL        :: ierrf(:)               
2425  INTEGER, INTENT(OUT), OPTIONAL        :: xnacc(:)               
2426  INTEGER, INTENT(OUT), OPTIONAL        :: xnrej(:)               
2427  INTEGER, INTENT(INOUT), OPTIONAL      :: istatus(:)             
2428  INTEGER, INTENT(IN), OPTIONAL         :: pe                     
2429  LOGICAL, INTENT(IN), OPTIONAL         :: l_debug                 
2430  INTEGER, DIMENSION(nkppctrl), INTENT(IN), OPTIONAL  :: icntrl_i         
2431  REAL(dp), DIMENSION(nkppctrl), INTENT(IN), OPTIONAL  :: rcntrl_i         
2432                                                                   
2433  INTEGER                                 :: k   ! loop variable     
2434  REAL(dp)                               :: dt                     
2435  INTEGER, DIMENSION(20)                :: istatus_u               
2436  INTEGER                                :: ierr_u                 
2437  INTEGER                                :: istatf                 
2438  INTEGER                                :: vl_dim_lo               
2439                                                                   
2440                                                                   
2441  IF (PRESENT (istatus)) istatus = 0                             
2442  IF (PRESENT (icntrl_i)) icntrl  = icntrl_i                     
2443  IF (PRESENT (rcntrl_i)) rcntrl  = rcntrl_i                     
2444                                                                   
2445  vl_glo = size(tempi, 1)                                           
2446                                                                   
2447  vl_dim_lo = vl_dim                                               
2448  DO k=1, vl_glo, vl_dim_lo                                           
2449    is = k                                                         
2450    ie = min(k+ vl_dim_lo-1, vl_glo)                                 
2451    vl = ie-is+ 1                                                   
2452                                                                   
2453    c(:) = conc(is, :)                                           
2454                                                                   
2455    temp = tempi(is)                                             
2456                                                                   
2457    qvap = qvapi(is)                                             
2458                                                                   
2459    fakt = fakti(is)                                             
2460                                                                   
2461    CALL initialize                                                 
2462                                                                   
2463    phot(:) = photo(is, :)                                           
2464                                                                   
2465    CALL update_rconst                                             
2466                                                                   
2467    dt = time_step_len                                             
2468                                                                   
2469    ! integrate from t=0 to t=dt                                   
2470    CALL integrate(0._dp, dt, icntrl, rcntrl, istatus_u = istatus_u, ierr_u=ierr_u)
2471                                                                   
2472                                                                   
2473   IF (PRESENT(l_debug) .AND. PRESENT(pe)) THEN                       
2474      IF (l_debug) CALL error_output(conc(is, :), ierr_u, pe)         
2475   ENDIF                                                             
2476                                                                     
2477    conc(is, :) = c(:)                                               
2478                                                                   
2479    ! RETURN diagnostic information                                 
2480                                                                   
2481    IF (PRESENT(ierrf))   ierrf(is) = ierr_u                     
2482    IF (PRESENT(xnacc))   xnacc(is) = istatus_u(4)               
2483    IF (PRESENT(xnrej))   xnrej(is) = istatus_u(5)               
2484                                                                   
2485    IF (PRESENT (istatus)) THEN                                   
2486      istatus(1:8) = istatus(1:8) + istatus_u(1:8)               
2487    ENDIF                                                         
2488                                                                   
2489  END DO                                                           
2490 
2491                                                                   
2492! Deallocate input arrays                                           
2493                                                                   
2494                                                                   
2495  data_loaded = .FALSE.                                             
2496                                                                   
2497  RETURN                                                           
2498END SUBROUTINE chem_gasphase_integrate                                       
2499
2500END MODULE chem_gasphase_mod
2501 
Note: See TracBrowser for help on using the repository browser.