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

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

Comments in all chem_gasphase_mod.kpp, add def_salsagas/*.f90, removed executables

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