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

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

Removed unused variables from chem_gasphase_mod.f90

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