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

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

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

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