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

Last change on this file since 3458 was 3458, checked in by kanani, 6 years ago

Reintegrated fixes/changes from branch chemistry

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