source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_simple/chem_gasphase_mod.f90 @ 3639

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

Bug fix of rate of eq.2 in simple.eqn and derived mechanisms, update of corresponding chem_gasphase_mod.f90

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