source: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_salsa+simple/chem_gasphase_mod.f90 @ 3780

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

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

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