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

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

renaming of get_mechanismname, do_emiss and do_depo, sorting in chem_modules

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