source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_smog/chem_gasphase_mod.f90 @ 3697

Last change on this file since 3697 was 3681, checked in by hellstea, 6 years ago

Major update of pmc_interface_mod

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