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

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

removed read from unit 10 in chemistry_model_mod.f90, added get_mechanismname

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