Index: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_simple/chem_gasphase_mod.f90
===================================================================
--- palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_simple/chem_gasphase_mod.f90	(revision 3298)
+++ palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_simple/chem_gasphase_mod.f90	(revision 3298)
@@ -0,0 +1,2513 @@
+MODULE chem_gasphase_mod
+ 
+!   Mechanism: simple
+!
+!------------------------------------------------------------------------------!
+!
+! ******Module chem_gasphase_mod is automatically generated by kpp4palm ******
+!
+!   *********Please do NOT change this Code,it will be ovewritten *********
+!
+!------------------------------------------------------------------------------!
+! This file was created by KPP (http://people.cs.vt.edu/asandu/Software/Kpp/)
+! and kpp4palm (created by Klaus Ketelsen). kpp4palm is an adapted version 
+! of KP4 (JÃ¶ckel,P.,Kerkweg,A.,Pozzer,A.,Sander,R.,Tost,H.,Riede,
+! H.,Baumgaertner,A.,Gromov,S.,and Kern,B.,2010: Development cycle 2 of
+! the Modular Earth Submodel System (MESSy2),Geosci. Model Dev.,3,717-752,
+! https://doi.org/10.5194/gmd-3-717-2010). KP4 is part of the Modular Earth 
+! Submodel System (MESSy),which is is available under the  GNU General Public 
+! License (GPL).
+!
+! KPP is free software; you can redistribute it and/or modify it under the terms 
+! of the General Public Licence as published by the Free Software Foundation; 
+! either version 2 of the License,or (at your option) any later version. 
+! KPP is distributed in the hope that it will be useful,but WITHOUT ANY WARRANTY; 
+! without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
+! PURPOSE. See the GNU General Public Licence for more details. 
+!
+!------------------------------------------------------------------------------!
+! This file is part of the PALM model system.
+!
+! PALM is free software: you can redistribute it and/or modify it under the
+! terms of the GNU General Public License as published by the Free Software
+! Foundation,either version 3 of the License,or (at your option) any later
+! version.
+!
+! PALM is distributed in the hope that it will be useful,but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with
+! PALM. If not,see <http://www.gnu.org/licenses/>.
+!
+! Copyright 1997-2018 Leibniz Universitaet Hannover
+!--------------------------------------------------------------------------------!
+!
+! Current revisions:
+! ------------------
+! 
+! 
+! Former revisions:
+! -----------------
+! $Id: module_header 2460 2017-09-13 14:47:48Z forkel $
+! forkel June 2018: qvap,fakt added
+! forkel June 2018: reset case in  Initialize,Integrate,Update_rconst
+! 
+! 
+! 2460 2017-09-13 14:47:48Z forkel
+!
+! forkel Sept. 2017: Variables for photolyis added
+!
+!
+! Nov. 2016: Intial version (Klaus Ketelsen)
+!
+!------------------------------------------------------------------------------!
+!
+
+
+! Set kpp Double Precision to PALM Default Precision
+
+  USE kinds,           ONLY: dp=>wp
+
+  USE pegrid,          ONLY: myid, threads_per_task
+
+  IMPLICIT        NONE
+  PRIVATE
+  !SAVE  ! note: occurs again in automatically generated code ...
+
+!  PUBLIC :: IERR_NAMES
+ 
+! PUBLIC :: SPC_NAMES,EQN_NAMES,EQN_TAGS,REQ_HET,REQ_AEROSOL,REQ_PHOTRAT &
+!         ,REQ_MCFCT,IP_MAX,jname
+
+  PUBLIC :: eqn_names, phot_names, spc_names
+  PUBLIC :: nmaxfixsteps
+  PUBLIC :: atol, rtol
+  PUBLIC :: nspec, nreact
+  PUBLIC :: temp
+  PUBLIC :: qvap
+  PUBLIC :: fakt
+  PUBLIC :: phot 
+  PUBLIC :: rconst
+  PUBLIC :: nvar
+  PUBLIC :: nphot
+  PUBLIC :: vl_dim                     ! PUBLIC to ebable other MODULEs to distiguish between scalar and vec
+ 
+  PUBLIC :: initialize, integrate, update_rconst
+  PUBLIC :: chem_gasphase_integrate
+  PUBLIC :: initialize_kpp_ctrl
+
+! END OF MODULE HEADER TEMPLATE
+                                                                 
+! Variables used for vector mode                                 
+                                                                 
+  LOGICAL, PARAMETER          :: l_vector = .FALSE.            
+  INTEGER, PARAMETER          :: i_lu_di = 2
+  INTEGER, PARAMETER          :: vl_dim = 1
+  INTEGER                     :: vl                              
+                                                                 
+  INTEGER                     :: vl_glo                          
+  INTEGER                     :: is, ie                           
+                                                                 
+                                                                 
+  INTEGER, DIMENSION(vl_dim)  :: kacc, krej                       
+  INTEGER, DIMENSION(vl_dim)  :: ierrv                           
+  LOGICAL                     :: data_loaded = .FALSE.             
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! 
+! Parameter Module File
+! 
+! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
+!       (http://www.cs.vt.edu/~asandu/Software/KPP)
+! KPP is distributed under GPL,the general public licence
+!       (http://www.gnu.org/copyleft/gpl.html)
+! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
+! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
+!     With important contributions from:
+!        M. Damian,Villanova University,USA
+!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
+! 
+! File                 : chem_gasphase_mod_Parameters.f90
+! Time                 : Tue Sep 25 18:35:44 2018
+! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
+! Equation file        : chem_gasphase_mod.kpp
+! Output root filename : chem_gasphase_mod
+! 
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+
+
+
+
+! NSPEC - Number of chemical species
+  INTEGER, PARAMETER :: nspec = 11 
+! NVAR - Number of Variable species
+  INTEGER, PARAMETER :: nvar = 9 
+! NVARACT - Number of Active species
+  INTEGER, PARAMETER :: nvaract = 7 
+! NFIX - Number of Fixed species
+  INTEGER, PARAMETER :: nfix = 2 
+! NREACT - Number of reactions
+  INTEGER, PARAMETER :: nreact = 7 
+! NVARST - Starting of variables in conc. vect.
+  INTEGER, PARAMETER :: nvarst = 1 
+! NFIXST - Starting of fixed in conc. vect.
+  INTEGER, PARAMETER :: nfixst = 10 
+! NONZERO - Number of nonzero entries in Jacobian
+  INTEGER, PARAMETER :: nonzero = 35 
+! LU_NONZERO - Number of nonzero entries in LU factoriz. of Jacobian
+  INTEGER, PARAMETER :: lu_nonzero = 37 
+! CNVAR - (NVAR+1) Number of elements in compressed row format
+  INTEGER, PARAMETER :: cnvar = 10 
+! CNEQN - (NREACT+1) Number stoicm elements in compressed col format
+  INTEGER, PARAMETER :: cneqn = 8 
+! NHESS - Length of Sparse Hessian
+  INTEGER, PARAMETER :: nhess = 18 
+! NMASS - Number of atoms to check mass balance
+  INTEGER, PARAMETER :: nmass = 1 
+
+! Index declaration for variable species in C and VAR
+!   VAR(ind_spc) = C(ind_spc)
+
+  INTEGER, PARAMETER, PUBLIC :: ind_hno3 = 1 
+  INTEGER, PARAMETER, PUBLIC :: ind_rcho = 2 
+  INTEGER, PARAMETER, PUBLIC :: ind_rh = 3 
+  INTEGER, PARAMETER, PUBLIC :: ind_ho2 = 4 
+  INTEGER, PARAMETER, PUBLIC :: ind_ro2 = 5 
+  INTEGER, PARAMETER, PUBLIC :: ind_oh = 6 
+  INTEGER, PARAMETER, PUBLIC :: ind_no2 = 7 
+  INTEGER, PARAMETER, PUBLIC :: ind_o3 = 8 
+  INTEGER, PARAMETER, PUBLIC :: ind_no = 9 
+
+! Index declaration for fixed species in C
+!   C(ind_spc)
+
+  INTEGER, PARAMETER, PUBLIC :: ind_h2o = 10 
+  INTEGER, PARAMETER, PUBLIC :: ind_o2 = 11 
+
+! Index declaration for fixed species in FIX
+!    FIX(indf_spc) = C(ind_spc) = C(NVAR+indf_spc)
+
+  INTEGER, PARAMETER :: indf_h2o = 1 
+  INTEGER, PARAMETER :: indf_o2 = 2 
+
+! NJVRP - Length of sparse Jacobian JVRP
+  INTEGER, PARAMETER :: njvrp = 12 
+
+! NSTOICM - Length of Sparse Stoichiometric Matrix
+  INTEGER, PARAMETER :: nstoicm = 23 
+
+
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! 
+! Global Data Module File
+! 
+! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
+!       (http://www.cs.vt.edu/~asandu/Software/KPP)
+! KPP is distributed under GPL,the general public licence
+!       (http://www.gnu.org/copyleft/gpl.html)
+! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
+! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
+!     With important contributions from:
+!        M. Damian,Villanova University,USA
+!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
+! 
+! File                 : chem_gasphase_mod_Global.f90
+! Time                 : Tue Sep 25 18:35:44 2018
+! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
+! Equation file        : chem_gasphase_mod.kpp
+! Output root filename : chem_gasphase_mod
+! 
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+
+
+
+
+! Declaration of global variables
+
+! C - Concentration of all species
+  REAL(kind=dp):: c(nspec)
+! VAR - Concentrations of variable species (global)
+  REAL(kind=dp):: var(nvar)
+! FIX - Concentrations of fixed species (global)
+  REAL(kind=dp):: fix(nfix)
+! VAR,FIX are chunks of array C
+      EQUIVALENCE( c(1), var(1))
+      EQUIVALENCE( c(10), fix(1))
+! RCONST - Rate constants (global)
+  REAL(kind=dp):: rconst(nreact)
+! TIME - Current integration time
+  REAL(kind=dp):: time
+! TEMP - Temperature
+  REAL(kind=dp):: temp
+! TSTART - Integration start time
+  REAL(kind=dp):: tstart
+! ATOL - Absolute tolerance
+  REAL(kind=dp):: atol(nvar)
+! RTOL - Relative tolerance
+  REAL(kind=dp):: rtol(nvar)
+! STEPMIN - Lower bound for integration step
+  REAL(kind=dp):: stepmin
+! CFACTOR - Conversion factor for concentration units
+  REAL(kind=dp):: cfactor
+
+! INLINED global variable declarations
+
+! QVAP - Water vapor
+  REAL(kind=dp):: qvap
+! FAKT - Conversion factor
+  REAL(kind=dp):: fakt
+
+
+! INLINED global variable declarations
+
+
+
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! 
+! Sparse Jacobian Data Structures File
+! 
+! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
+!       (http://www.cs.vt.edu/~asandu/Software/KPP)
+! KPP is distributed under GPL,the general public licence
+!       (http://www.gnu.org/copyleft/gpl.html)
+! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
+! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
+!     With important contributions from:
+!        M. Damian,Villanova University,USA
+!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
+! 
+! File                 : chem_gasphase_mod_JacobianSP.f90
+! Time                 : Tue Sep 25 18:35:44 2018
+! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
+! Equation file        : chem_gasphase_mod.kpp
+! Output root filename : chem_gasphase_mod
+! 
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+
+
+
+
+! Sparse Jacobian Data
+
+
+  INTEGER, PARAMETER, DIMENSION(37):: lu_irow =  (/ &
+       1, 1, 1, 2, 2, 2, 3, 3, 4, 4, 4, 5, &
+       5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 7, 7, &
+       7, 7, 7, 7, 8, 8, 8, 9, 9, 9, 9, 9, &
+       9 /) 
+
+  INTEGER, PARAMETER, DIMENSION(37):: lu_icol =  (/ &
+       1, 6, 7, 2, 5, 9, 3, 6, 4, 5, 9, 3, &
+       5, 6, 9, 3, 4, 5, 6, 7, 8, 9, 4, 5, &
+       6, 7, 8, 9, 7, 8, 9, 4, 5, 6, 7, 8, &
+       9 /) 
+
+  INTEGER, PARAMETER, DIMENSION(10):: lu_crow =  (/ &
+       1, 4, 7, 9, 12, 16, 23, 29, 32, 38 /) 
+
+  INTEGER, PARAMETER, DIMENSION(10):: lu_diag =  (/ &
+       1, 4, 7, 9, 13, 19, 26, 30, 37, 38 /) 
+
+
+
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! 
+! Utility Data Module File
+! 
+! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
+!       (http://www.cs.vt.edu/~asandu/Software/KPP)
+! KPP is distributed under GPL,the general public licence
+!       (http://www.gnu.org/copyleft/gpl.html)
+! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
+! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
+!     With important contributions from:
+!        M. Damian,Villanova University,USA
+!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
+! 
+! File                 : chem_gasphase_mod_Monitor.f90
+! Time                 : Tue Sep 25 18:35:44 2018
+! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
+! Equation file        : chem_gasphase_mod.kpp
+! Output root filename : chem_gasphase_mod
+! 
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+
+
+
+  CHARACTER(len=15), PARAMETER, DIMENSION(11):: spc_names =  (/ &
+     'HNO3           ','RCHO           ','RH             ',&
+     'HO2            ','RO2            ','OH             ',&
+     'NO2            ','O3             ','NO             ',&
+     'H2O            ','O2             ' /)
+
+  CHARACTER(len=100), PARAMETER, DIMENSION(7):: eqn_names =  (/ &
+     '     NO2 --> O3 + NO                                                                                ',&
+     '      O3 --> 2 OH + O2                                                                              ',&
+     ' O3 + NO --> NO2                                                                                    ',&
+     ' RH + OH --> RO2 + H2O                                                                              ',&
+     'RO2 + NO --> RCHO + HO2 + NO2                                                                       ',&
+     'HO2 + NO --> OH + NO2                                                                               ',&
+     'OH + NO2 --> HNO3                                                                                   ' /)
+
+! INLINED global variables
+
+  !   inline f90_data: declaration of global variables for photolysis
+  !   REAL(kind=dp):: phot(nphot)must eventually be moved to global later for 
+  INTEGER, PARAMETER :: nphot = 2
+  !   phot photolysis frequencies 
+  REAL(kind=dp):: phot(nphot)
+
+  INTEGER, PARAMETER, PUBLIC :: j_no2 = 1
+  INTEGER, PARAMETER, PUBLIC :: j_o31d = 2
+
+  CHARACTER(len=15), PARAMETER, DIMENSION(nphot):: phot_names =   (/ &
+     'J_NO2          ','J_O31D         '/)
+
+! End INLINED global variables
+
+
+! Automatic generated PUBLIC Statements for ip_ and ihs_ variables 
+ 
+! Automatic generated PUBLIC Statements for ip_ and ihs_ variables 
+ 
+! Automatic generated PUBLIC Statements for ip_ and ihs_ variables 
+ 
+! Automatic generated PUBLIC Statements for ip_ and ihs_ variables 
+ 
+ 
+!  variable definations from  individual module headers 
+ 
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! 
+! Initialization File
+! 
+! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
+!       (http://www.cs.vt.edu/~asandu/Software/KPP)
+! KPP is distributed under GPL,the general public licence
+!       (http://www.gnu.org/copyleft/gpl.html)
+! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
+! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
+!     With important contributions from:
+!        M. Damian,Villanova University,USA
+!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
+! 
+! File                 : chem_gasphase_mod_Initialize.f90
+! Time                 : Tue Sep 25 18:35:44 2018
+! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
+! Equation file        : chem_gasphase_mod.kpp
+! Output root filename : chem_gasphase_mod
+! 
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+
+
+
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! 
+! Numerical Integrator (Time-Stepping) File
+! 
+! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
+!       (http://www.cs.vt.edu/~asandu/Software/KPP)
+! KPP is distributed under GPL,the general public licence
+!       (http://www.gnu.org/copyleft/gpl.html)
+! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
+! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
+!     With important contributions from:
+!        M. Damian,Villanova University,USA
+!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
+! 
+! File                 : chem_gasphase_mod_Integrator.f90
+! Time                 : Tue Sep 25 18:35:44 2018
+! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
+! Equation file        : chem_gasphase_mod.kpp
+! Output root filename : chem_gasphase_mod
+! 
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! 
+! INTEGRATE - Integrator routine
+!   Arguments :
+!      TIN       - Start Time for Integration
+!      TOUT      - End Time for Integration
+! 
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
+!  Rosenbrock - Implementation of several Rosenbrock methods:             !
+!               *Ros2                                                    !
+!               *Ros3                                                    !
+!               *Ros4                                                    !
+!               *Rodas3                                                  !
+!               *Rodas4                                                  !
+!  By default the code employs the KPP sparse linear algebra routines     !
+!  Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK)        !
+!                                                                         !
+!    (C)  Adrian Sandu,August 2004                                       !
+!    Virginia Polytechnic Institute and State University                  !
+!    Contact: sandu@cs.vt.edu                                             !
+!    Revised by Philipp Miehe and Adrian Sandu,May 2006                  !                               !
+!    This implementation is part of KPP - the Kinetic PreProcessor        !
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!
+
+
+  SAVE
+  
+!~~~>  statistics on the work performed by the rosenbrock method
+  INTEGER, PARAMETER :: nfun=1, njac=2, nstp=3, nacc=4, &
+                        nrej=5, ndec=6, nsol=7, nsng=8, &
+                        ntexit=1, nhexit=2, nhnew = 3
+
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! 
+! Linear Algebra Data and Routines File
+! 
+! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
+!       (http://www.cs.vt.edu/~asandu/Software/KPP)
+! KPP is distributed under GPL,the general public licence
+!       (http://www.gnu.org/copyleft/gpl.html)
+! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
+! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
+!     With important contributions from:
+!        M. Damian,Villanova University,USA
+!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
+! 
+! File                 : chem_gasphase_mod_LinearAlgebra.f90
+! Time                 : Tue Sep 25 18:35:44 2018
+! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
+! Equation file        : chem_gasphase_mod.kpp
+! Output root filename : chem_gasphase_mod
+! 
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+
+
+
+
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! 
+! The ODE Jacobian of Chemical Model File
+! 
+! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
+!       (http://www.cs.vt.edu/~asandu/Software/KPP)
+! KPP is distributed under GPL,the general public licence
+!       (http://www.gnu.org/copyleft/gpl.html)
+! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
+! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
+!     With important contributions from:
+!        M. Damian,Villanova University,USA
+!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
+! 
+! File                 : chem_gasphase_mod_Jacobian.f90
+! Time                 : Tue Sep 25 18:35:44 2018
+! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
+! Equation file        : chem_gasphase_mod.kpp
+! Output root filename : chem_gasphase_mod
+! 
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+
+
+
+
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! 
+! The ODE Function of Chemical Model File
+! 
+! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
+!       (http://www.cs.vt.edu/~asandu/Software/KPP)
+! KPP is distributed under GPL,the general public licence
+!       (http://www.gnu.org/copyleft/gpl.html)
+! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
+! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
+!     With important contributions from:
+!        M. Damian,Villanova University,USA
+!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
+! 
+! File                 : chem_gasphase_mod_Function.f90
+! Time                 : Tue Sep 25 18:35:44 2018
+! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
+! Equation file        : chem_gasphase_mod.kpp
+! Output root filename : chem_gasphase_mod
+! 
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+
+
+
+! A - Rate for each equation
+  REAL(kind=dp):: a(nreact)
+
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! 
+! The Reaction Rates File
+! 
+! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
+!       (http://www.cs.vt.edu/~asandu/Software/KPP)
+! KPP is distributed under GPL,the general public licence
+!       (http://www.gnu.org/copyleft/gpl.html)
+! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
+! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
+!     With important contributions from:
+!        M. Damian,Villanova University,USA
+!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
+! 
+! File                 : chem_gasphase_mod_Rates.f90
+! Time                 : Tue Sep 25 18:35:44 2018
+! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
+! Equation file        : chem_gasphase_mod.kpp
+! Output root filename : chem_gasphase_mod
+! 
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+
+
+
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! 
+! Auxiliary Routines File
+! 
+! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
+!       (http://www.cs.vt.edu/~asandu/Software/KPP)
+! KPP is distributed under GPL,the general public licence
+!       (http://www.gnu.org/copyleft/gpl.html)
+! (C) 1995-1997,V. Damian & A. Sandu,CGRER,Univ. Iowa
+! (C) 1997-2005,A. Sandu,Michigan Tech,Virginia Tech
+!     With important contributions from:
+!        M. Damian,Villanova University,USA
+!        R. Sander,Max-Planck Institute for Chemistry,Mainz,Germany
+! 
+! File                 : chem_gasphase_mod_Util.f90
+! Time                 : Tue Sep 25 18:35:44 2018
+! Working directory    : /home/forkel-r/palmstuff/work/chemistry20180925/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm
+! Equation file        : chem_gasphase_mod.kpp
+! Output root filename : chem_gasphase_mod
+! 
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+
+
+
+
+  ! header MODULE initialize_kpp_ctrl_template
+
+  ! notes:
+  ! - l_vector is automatically defined by kp4
+  ! - vl_dim is automatically defined by kp4
+  ! - i_lu_di is automatically defined by kp4
+  ! - wanted is automatically defined by xmecca
+  ! - icntrl rcntrl are automatically defined by kpp
+  ! - "USE messy_main_tools" is in MODULE_header of messy_mecca_kpp.f90
+  ! - SAVE will be automatically added by kp4
+
+  !SAVE
+
+  ! for fixed time step control
+  ! ... max. number of fixed time steps (sum must be 1)
+  INTEGER, PARAMETER         :: nmaxfixsteps = 50
+  ! ... switch for fixed time stepping
+  LOGICAL, PUBLIC            :: l_fixed_step = .FALSE.
+  INTEGER, PUBLIC            :: nfsteps = 1
+  ! ... number of kpp control PARAMETERs
+  INTEGER, PARAMETER, PUBLIC :: nkppctrl = 20
+  !
+  INTEGER, DIMENSION(nkppctrl), PUBLIC     :: icntrl = 0
+  REAL(dp), DIMENSION(nkppctrl), PUBLIC     :: rcntrl = 0.0_dp
+  REAL(dp), DIMENSION(nmaxfixsteps), PUBLIC :: t_steps = 0.0_dp
+
+  ! END header MODULE initialize_kpp_ctrl_template
+
+ 
+! Interface Block 
+ 
+  INTERFACE            initialize
+    MODULE PROCEDURE   initialize
+  END INTERFACE        initialize
+ 
+  INTERFACE            integrate
+    MODULE PROCEDURE   integrate
+  END INTERFACE        integrate
+ 
+  INTERFACE            fun
+    MODULE PROCEDURE   fun
+  END INTERFACE        fun
+ 
+  INTERFACE            kppsolve
+    MODULE PROCEDURE   kppsolve
+  END INTERFACE        kppsolve
+ 
+  INTERFACE            jac_sp
+    MODULE PROCEDURE   jac_sp
+  END INTERFACE        jac_sp
+ 
+  INTERFACE            k_arr
+    MODULE PROCEDURE   k_arr
+  END INTERFACE        k_arr
+ 
+  INTERFACE            update_rconst
+    MODULE PROCEDURE   update_rconst
+  END INTERFACE        update_rconst
+ 
+  INTERFACE            arr2
+    MODULE PROCEDURE   arr2
+  END INTERFACE        arr2
+ 
+  INTERFACE            initialize_kpp_ctrl
+    MODULE PROCEDURE   initialize_kpp_ctrl
+  END INTERFACE        initialize_kpp_ctrl
+ 
+  INTERFACE            error_output
+    MODULE PROCEDURE   error_output
+  END INTERFACE        error_output
+ 
+  INTERFACE            wscal
+    MODULE PROCEDURE   wscal
+  END INTERFACE        wscal
+ 
+!INTERFACE not working  INTERFACE            waxpy
+!INTERFACE not working    MODULE PROCEDURE   waxpy
+!INTERFACE not working  END INTERFACE        waxpy
+ 
+  INTERFACE            rosenbrock
+    MODULE PROCEDURE   rosenbrock
+  END INTERFACE        rosenbrock
+ 
+  INTERFACE            funtemplate
+    MODULE PROCEDURE   funtemplate
+  END INTERFACE        funtemplate
+ 
+  INTERFACE            jactemplate
+    MODULE PROCEDURE   jactemplate
+  END INTERFACE        jactemplate
+ 
+  INTERFACE            kppdecomp
+    MODULE PROCEDURE   kppdecomp
+  END INTERFACE        kppdecomp
+ 
+  INTERFACE            chem_gasphase_integrate
+    MODULE PROCEDURE   chem_gasphase_integrate
+  END INTERFACE        chem_gasphase_integrate
+ 
+ 
+ CONTAINS
+ 
+SUBROUTINE initialize()
+
+
+  INTEGER         :: j, k
+
+  INTEGER :: i
+  REAL(kind=dp):: x
+  k = is
+  cfactor = 1.000000e+00_dp
+
+  x = (0.) * cfactor
+  DO i = 1 , nvar
+  ENDDO
+
+  x = (0.) * cfactor
+  DO i = 1 , nfix
+    fix(i) = x
+  ENDDO
+
+! constant rate coefficients
+! END constant rate coefficients
+
+! INLINED initializations
+
+! End INLINED initializations
+
+      
+END SUBROUTINE initialize
+ 
+SUBROUTINE integrate( tin, tout, &
+  icntrl_u, rcntrl_u, istatus_u, rstatus_u, ierr_u)
+
+
+   REAL(kind=dp), INTENT(IN):: tin  ! start time
+   REAL(kind=dp), INTENT(IN):: tout ! END time
+   ! OPTIONAL input PARAMETERs and statistics
+   INTEGER,      INTENT(IN), OPTIONAL :: icntrl_u(20)
+   REAL(kind=dp), INTENT(IN), OPTIONAL :: rcntrl_u(20)
+   INTEGER,      INTENT(OUT), OPTIONAL :: istatus_u(20)
+   REAL(kind=dp), INTENT(OUT), OPTIONAL :: rstatus_u(20)
+   INTEGER,      INTENT(OUT), OPTIONAL :: ierr_u
+
+   REAL(kind=dp):: rcntrl(20), rstatus(20)
+   INTEGER       :: icntrl(20), istatus(20), ierr
+
+   INTEGER, SAVE :: ntotal = 0
+
+   icntrl(:) = 0
+   rcntrl(:) = 0.0_dp
+   istatus(:) = 0
+   rstatus(:) = 0.0_dp
+
+    !~~~> fine-tune the integrator:
+   icntrl(1) = 0      ! 0 - non- autonomous, 1 - autonomous
+   icntrl(2) = 0      ! 0 - vector tolerances, 1 - scalars
+
+   ! IF OPTIONAL PARAMETERs are given, and IF they are >0, 
+   ! THEN they overwrite default settings. 
+   IF (PRESENT(icntrl_u))THEN
+     WHERE(icntrl_u(:)> 0)icntrl(:) = icntrl_u(:)
+   ENDIF
+   IF (PRESENT(rcntrl_u))THEN
+     WHERE(rcntrl_u(:)> 0)rcntrl(:) = rcntrl_u(:)
+   ENDIF
+
+
+   CALL rosenbrock(nvar, var, tin, tout,  &
+         atol, rtol,               &
+         rcntrl, icntrl, rstatus, istatus, ierr)
+
+   !~~~> debug option: show no of steps
+   ! ntotal = ntotal + istatus(nstp)
+   ! PRINT*,'NSTEPS=',ISTATUS(Nstp),' (',Ntotal,')','  O3=',VAR(ind_O3)
+
+   stepmin = rstatus(nhexit)
+   ! IF OPTIONAL PARAMETERs are given for output they 
+   ! are updated with the RETURN information
+   IF (PRESENT(istatus_u))istatus_u(:) = istatus(:)
+   IF (PRESENT(rstatus_u))rstatus_u(:) = rstatus(:)
+   IF (PRESENT(ierr_u))  ierr_u       = ierr
+
+END SUBROUTINE integrate
+ 
+SUBROUTINE fun(v, f, rct, vdot)
+
+! V - Concentrations of variable species (local)
+  REAL(kind=dp):: v(nvar)
+! F - Concentrations of fixed species (local)
+  REAL(kind=dp):: f(nfix)
+! RCT - Rate constants (local)
+  REAL(kind=dp):: rct(nreact)
+! Vdot - Time derivative of variable species concentrations
+  REAL(kind=dp):: vdot(nvar)
+
+
+! Computation of equation rates
+  a(1) = rct(1) * v(7)
+  a(2) = rct(2) * v(8)
+  a(3) = rct(3) * v(8) * v(9)
+  a(4) = rct(4) * v(3) * v(6)
+  a(5) = rct(5) * v(5) * v(9)
+  a(6) = rct(6) * v(4) * v(9)
+  a(7) = rct(7) * v(6) * v(7)
+
+! Aggregate function
+  vdot(1) = a(7)
+  vdot(2) = a(5)
+  vdot(3) = - a(4)
+  vdot(4) = a(5) - a(6)
+  vdot(5) = a(4) - a(5)
+  vdot(6) = 2* a(2) - a(4) + a(6) - a(7)
+  vdot(7) = - a(1) + a(3) + a(5) + a(6) - a(7)
+  vdot(8) = a(1) - a(2) - a(3)
+  vdot(9) = a(1) - a(3) - a(5) - a(6)
+      
+END SUBROUTINE fun
+ 
+SUBROUTINE kppsolve(jvs, x)
+
+! JVS - sparse Jacobian of variables
+  REAL(kind=dp):: jvs(lu_nonzero)
+! X - Vector for variables
+  REAL(kind=dp):: x(nvar)
+
+  x(5) = x(5) - jvs(12) * x(3)
+  x(6) = x(6) - jvs(16) * x(3) - jvs(17) * x(4) - jvs(18) * x(5)
+  x(7) = x(7) - jvs(23) * x(4) - jvs(24) * x(5) - jvs(25) * x(6)
+  x(8) = x(8) - jvs(29) * x(7)
+  x(9) = x(9) - jvs(32) * x(4) - jvs(33) * x(5) - jvs(34) * x(6) - jvs(35) * x(7) - jvs(36) * x(8)
+  x(9) = x(9) / jvs(37)
+  x(8) = (x(8) - jvs(31) * x(9)) /(jvs(30))
+  x(7) = (x(7) - jvs(27) * x(8) - jvs(28) * x(9)) /(jvs(26))
+  x(6) = (x(6) - jvs(20) * x(7) - jvs(21) * x(8) - jvs(22) * x(9)) /(jvs(19))
+  x(5) = (x(5) - jvs(14) * x(6) - jvs(15) * x(9)) /(jvs(13))
+  x(4) = (x(4) - jvs(10) * x(5) - jvs(11) * x(9)) /(jvs(9))
+  x(3) = (x(3) - jvs(8) * x(6)) /(jvs(7))
+  x(2) = (x(2) - jvs(5) * x(5) - jvs(6) * x(9)) /(jvs(4))
+  x(1) = (x(1) - jvs(2) * x(6) - jvs(3) * x(7)) /(jvs(1))
+      
+END SUBROUTINE kppsolve
+ 
+SUBROUTINE jac_sp(v, f, rct, jvs)
+
+! V - Concentrations of variable species (local)
+  REAL(kind=dp):: v(nvar)
+! F - Concentrations of fixed species (local)
+  REAL(kind=dp):: f(nfix)
+! RCT - Rate constants (local)
+  REAL(kind=dp):: rct(nreact)
+! JVS - sparse Jacobian of variables
+  REAL(kind=dp):: jvs(lu_nonzero)
+
+
+! Local variables
+! B - Temporary array
+  REAL(kind=dp):: b(12)
+
+! B(1) = dA(1)/dV(7)
+  b(1) = rct(1)
+! B(2) = dA(2)/dV(8)
+  b(2) = rct(2)
+! B(3) = dA(3)/dV(8)
+  b(3) = rct(3) * v(9)
+! B(4) = dA(3)/dV(9)
+  b(4) = rct(3) * v(8)
+! B(5) = dA(4)/dV(3)
+  b(5) = rct(4) * v(6)
+! B(6) = dA(4)/dV(6)
+  b(6) = rct(4) * v(3)
+! B(7) = dA(5)/dV(5)
+  b(7) = rct(5) * v(9)
+! B(8) = dA(5)/dV(9)
+  b(8) = rct(5) * v(5)
+! B(9) = dA(6)/dV(4)
+  b(9) = rct(6) * v(9)
+! B(10) = dA(6)/dV(9)
+  b(10) = rct(6) * v(4)
+! B(11) = dA(7)/dV(6)
+  b(11) = rct(7) * v(7)
+! B(12) = dA(7)/dV(7)
+  b(12) = rct(7) * v(6)
+
+! Construct the Jacobian terms from B's
+! JVS(1) = Jac_FULL(1,1)
+  jvs(1) = 0
+! JVS(2) = Jac_FULL(1,6)
+  jvs(2) = b(11)
+! JVS(3) = Jac_FULL(1,7)
+  jvs(3) = b(12)
+! JVS(4) = Jac_FULL(2,2)
+  jvs(4) = 0
+! JVS(5) = Jac_FULL(2,5)
+  jvs(5) = b(7)
+! JVS(6) = Jac_FULL(2,9)
+  jvs(6) = b(8)
+! JVS(7) = Jac_FULL(3,3)
+  jvs(7) = - b(5)
+! JVS(8) = Jac_FULL(3,6)
+  jvs(8) = - b(6)
+! JVS(9) = Jac_FULL(4,4)
+  jvs(9) = - b(9)
+! JVS(10) = Jac_FULL(4,5)
+  jvs(10) = b(7)
+! JVS(11) = Jac_FULL(4,9)
+  jvs(11) = b(8) - b(10)
+! JVS(12) = Jac_FULL(5,3)
+  jvs(12) = b(5)
+! JVS(13) = Jac_FULL(5,5)
+  jvs(13) = - b(7)
+! JVS(14) = Jac_FULL(5,6)
+  jvs(14) = b(6)
+! JVS(15) = Jac_FULL(5,9)
+  jvs(15) = - b(8)
+! JVS(16) = Jac_FULL(6,3)
+  jvs(16) = - b(5)
+! JVS(17) = Jac_FULL(6,4)
+  jvs(17) = b(9)
+! JVS(18) = Jac_FULL(6,5)
+  jvs(18) = 0
+! JVS(19) = Jac_FULL(6,6)
+  jvs(19) = - b(6) - b(11)
+! JVS(20) = Jac_FULL(6,7)
+  jvs(20) = - b(12)
+! JVS(21) = Jac_FULL(6,8)
+  jvs(21) = 2* b(2)
+! JVS(22) = Jac_FULL(6,9)
+  jvs(22) = b(10)
+! JVS(23) = Jac_FULL(7,4)
+  jvs(23) = b(9)
+! JVS(24) = Jac_FULL(7,5)
+  jvs(24) = b(7)
+! JVS(25) = Jac_FULL(7,6)
+  jvs(25) = - b(11)
+! JVS(26) = Jac_FULL(7,7)
+  jvs(26) = - b(1) - b(12)
+! JVS(27) = Jac_FULL(7,8)
+  jvs(27) = b(3)
+! JVS(28) = Jac_FULL(7,9)
+  jvs(28) = b(4) + b(8) + b(10)
+! JVS(29) = Jac_FULL(8,7)
+  jvs(29) = b(1)
+! JVS(30) = Jac_FULL(8,8)
+  jvs(30) = - b(2) - b(3)
+! JVS(31) = Jac_FULL(8,9)
+  jvs(31) = - b(4)
+! JVS(32) = Jac_FULL(9,4)
+  jvs(32) = - b(9)
+! JVS(33) = Jac_FULL(9,5)
+  jvs(33) = - b(7)
+! JVS(34) = Jac_FULL(9,6)
+  jvs(34) = 0
+! JVS(35) = Jac_FULL(9,7)
+  jvs(35) = b(1)
+! JVS(36) = Jac_FULL(9,8)
+  jvs(36) = - b(3)
+! JVS(37) = Jac_FULL(9,9)
+  jvs(37) = - b(4) - b(8) - b(10)
+      
+END SUBROUTINE jac_sp
+ 
+  elemental REAL(kind=dp)FUNCTION k_arr (k_298, tdep, temp)
+    ! arrhenius FUNCTION
+
+    REAL,    INTENT(IN):: k_298 ! k at t = 298.15k
+    REAL,    INTENT(IN):: tdep  ! temperature dependence
+    REAL(kind=dp), INTENT(IN):: temp  ! temperature
+
+    intrinsic exp
+
+    k_arr = k_298 * exp(tdep* (1._dp/temp- 3.3540e-3_dp))! 1/298.15=3.3540e-3
+
+  END FUNCTION k_arr
+ 
+SUBROUTINE update_rconst()
+ INTEGER         :: k
+
+  k = is
+
+! Begin INLINED RCONST
+
+
+! End INLINED RCONST
+
+  rconst(1) = (phot(j_no2))
+  rconst(2) = (phot(j_o31d))
+  rconst(3) = (arr2(1.8e-12_dp , 1370.0_dp , temp))
+  rconst(4) = (arr2(2.e-11_dp , 500.0_dp , temp))
+  rconst(5) = (arr2(4.2e-12_dp , -180.0_dp , temp))
+  rconst(6) = (arr2(3.7e-12_dp , -240.0_dp , temp))
+  rconst(7) = (arr2(1.15e-11_dp , 0.0_dp , temp))
+      
+END SUBROUTINE update_rconst
+ 
+!  END FUNCTION ARR2
+REAL(kind=dp)FUNCTION arr2( a0, b0, temp)
+    REAL(kind=dp):: temp
+    REAL(kind=dp):: a0, b0
+    arr2 = a0 * exp( - b0 / temp)
+END FUNCTION arr2
+ 
+SUBROUTINE initialize_kpp_ctrl(status)
+
+
+  ! i/o
+  INTEGER,         INTENT(OUT):: status
+
+  ! local
+  REAL(dp):: tsum
+  INTEGER  :: i
+
+  ! check fixed time steps
+  tsum = 0.0_dp
+  DO i=1, nmaxfixsteps
+     IF (t_steps(i)< tiny(0.0_dp))exit
+     tsum = tsum + t_steps(i)
+  ENDDO
+
+  nfsteps = i- 1
+
+  l_fixed_step = (nfsteps > 0).and.((tsum - 1.0)< tiny(0.0_dp))
+
+  IF (l_vector)THEN
+     WRITE(*,*) ' MODE             : VECTOR (LENGTH=',VL_DIM,')'
+  ELSE
+     WRITE(*,*) ' MODE             : SCALAR'
+  ENDIF
+  !
+  WRITE(*,*) ' DE-INDEXING MODE :',I_LU_DI
+  !
+  WRITE(*,*) ' ICNTRL           : ',icntrl
+  WRITE(*,*) ' RCNTRL           : ',rcntrl
+  !
+  ! note: this is ONLY meaningful for vectorized (kp4)rosenbrock- methods
+  IF (l_vector)THEN
+     IF (l_fixed_step)THEN
+        WRITE(*,*) ' TIME STEPS       : FIXED (',t_steps(1:nfsteps),')'
+     ELSE
+        WRITE(*,*) ' TIME STEPS       : AUTOMATIC'
+     ENDIF
+  ELSE
+     WRITE(*,*) ' TIME STEPS       : AUTOMATIC '//&
+          &'(t_steps (CTRL_KPP) ignored in SCALAR MODE)'
+  ENDIF
+  ! mz_pj_20070531- 
+
+  status = 0
+
+
+END SUBROUTINE initialize_kpp_ctrl
+ 
+SUBROUTINE error_output(c, ierr, pe)
+
+
+  INTEGER, INTENT(IN):: ierr
+  INTEGER, INTENT(IN):: pe
+  REAL(dp), DIMENSION(:), INTENT(IN):: c
+
+  write(6,*) 'ERROR in chem_gasphase_mod ',ierr,C(1)
+
+
+END SUBROUTINE error_output
+ 
+      SUBROUTINE wscal(n, alpha, x, incx)
+!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
+!     constant times a vector: x(1:N) <- Alpha*x(1:N) 
+!     only for incX=incY=1
+!     after BLAS
+!     replace this by the function from the optimized BLAS implementation:
+!         CALL SSCAL(N,Alpha,X,1) or  CALL DSCAL(N,Alpha,X,1)
+!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
+
+      INTEGER  :: i, incx, m, mp1, n
+      REAL(kind=dp) :: x(n), alpha
+      REAL(kind=dp), PARAMETER  :: zero=0.0_dp, one=1.0_dp
+
+      IF (alpha .eq. one)RETURN
+      IF (n .le. 0)RETURN
+
+      m = mod(n, 5)
+      IF ( m .ne. 0)THEN
+        IF (alpha .eq. (- one))THEN
+          DO i = 1, m
+            x(i) = - x(i)
+          ENDDO
+        ELSEIF (alpha .eq. zero)THEN
+          DO i = 1, m
+            x(i) = zero
+          ENDDO
+        ELSE
+          DO i = 1, m
+            x(i) = alpha* x(i)
+          ENDDO
+        ENDIF
+        IF ( n .lt. 5)RETURN
+      ENDIF
+      mp1 = m + 1
+      IF (alpha .eq. (- one))THEN
+        DO i = mp1, n, 5
+          x(i)   = - x(i)
+          x(i + 1) = - x(i + 1)
+          x(i + 2) = - x(i + 2)
+          x(i + 3) = - x(i + 3)
+          x(i + 4) = - x(i + 4)
+        ENDDO
+      ELSEIF (alpha .eq. zero)THEN
+        DO i = mp1, n, 5
+          x(i)   = zero
+          x(i + 1) = zero
+          x(i + 2) = zero
+          x(i + 3) = zero
+          x(i + 4) = zero
+        ENDDO
+      ELSE
+        DO i = mp1, n, 5
+          x(i)   = alpha* x(i)
+          x(i + 1) = alpha* x(i + 1)
+          x(i + 2) = alpha* x(i + 2)
+          x(i + 3) = alpha* x(i + 3)
+          x(i + 4) = alpha* x(i + 4)
+        ENDDO
+      ENDIF
+
+      END SUBROUTINE wscal
+ 
+      SUBROUTINE waxpy(n, alpha, x, incx, y, incy)
+!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
+!     constant times a vector plus a vector: y <- y + Alpha*x
+!     only for incX=incY=1
+!     after BLAS
+!     replace this by the function from the optimized BLAS implementation:
+!         CALL SAXPY(N,Alpha,X,1,Y,1) or  CALL DAXPY(N,Alpha,X,1,Y,1)
+!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
+
+      INTEGER  :: i, incx, incy, m, mp1, n
+      REAL(kind=dp):: x(n), y(n), alpha
+      REAL(kind=dp), PARAMETER :: zero = 0.0_dp
+
+      IF (alpha .eq. zero)RETURN
+      IF (n .le. 0)RETURN
+
+      m = mod(n, 4)
+      IF ( m .ne. 0)THEN
+        DO i = 1, m
+          y(i) = y(i) + alpha* x(i)
+        ENDDO
+        IF ( n .lt. 4)RETURN
+      ENDIF
+      mp1 = m + 1
+      DO i = mp1, n, 4
+        y(i) = y(i) + alpha* x(i)
+        y(i + 1) = y(i + 1) + alpha* x(i + 1)
+        y(i + 2) = y(i + 2) + alpha* x(i + 2)
+        y(i + 3) = y(i + 3) + alpha* x(i + 3)
+      ENDDO
+      
+      END SUBROUTINE waxpy
+ 
+SUBROUTINE rosenbrock(n, y, tstart, tend, &
+           abstol, reltol,             &
+           rcntrl, icntrl, rstatus, istatus, ierr)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!
+!    Solves the system y'=F(t,y) using a Rosenbrock method defined by:
+!
+!     G = 1/(H*gamma(1)) - Jac(t0,Y0)
+!     T_i = t0 + Alpha(i)*H
+!     Y_i = Y0 + \sum_{j=1}^{i-1} A(i,j)*K_j
+!     G *K_i = Fun( T_i,Y_i)+ \sum_{j=1}^S C(i,j)/H *K_j +
+!         gamma(i)*dF/dT(t0,Y0)
+!     Y1 = Y0 + \sum_{j=1}^S M(j)*K_j
+!
+!    For details on Rosenbrock methods and their implementation consult:
+!      E. Hairer and G. Wanner
+!      "Solving ODEs II. Stiff and differential-algebraic problems".
+!      Springer series in computational mathematics,Springer-Verlag,1996.
+!    The codes contained in the book inspired this implementation.
+!
+!    (C)  Adrian Sandu,August 2004
+!    Virginia Polytechnic Institute and State University
+!    Contact: sandu@cs.vt.edu
+!    Revised by Philipp Miehe and Adrian Sandu,May 2006                  
+!    This implementation is part of KPP - the Kinetic PreProcessor
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!
+!~~~>   input arguments:
+!
+!-     y(n)  = vector of initial conditions (at t=tstart)
+!-    [tstart, tend]  = time range of integration
+!     (if Tstart>Tend the integration is performed backwards in time)
+!-    reltol, abstol = user precribed accuracy
+!- SUBROUTINE  fun( t, y, ydot) = ode FUNCTION, 
+!                       returns Ydot = Y' = F(T,Y)
+!- SUBROUTINE  jac( t, y, jcb) = jacobian of the ode FUNCTION, 
+!                       returns Jcb = dFun/dY
+!-    icntrl(1:20)  = INTEGER inputs PARAMETERs
+!-    rcntrl(1:20)  = REAL inputs PARAMETERs
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!
+!~~~>     output arguments:
+!
+!-    y(n)  - > vector of final states (at t- >tend)
+!-    istatus(1:20) - > INTEGER output PARAMETERs
+!-    rstatus(1:20) - > REAL output PARAMETERs
+!-    ierr            - > job status upon RETURN
+!                        success (positive value) or
+!                        failure (negative value)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!
+!~~~>     input PARAMETERs:
+!
+!    Note: For input parameters equal to zero the default values of the
+!       corresponding variables are used.
+!
+!    ICNTRL(1) = 1: F = F(y)   Independent of T (AUTONOMOUS)
+!              = 0: F = F(t,y) Depends on T (NON-AUTONOMOUS)
+!
+!    ICNTRL(2) = 0: AbsTol,RelTol are N-dimensional vectors
+!              = 1: AbsTol,RelTol are scalars
+!
+!    ICNTRL(3)  -> selection of a particular Rosenbrock method
+!        = 0 :    Rodas3 (default)
+!        = 1 :    Ros2
+!        = 2 :    Ros3
+!        = 3 :    Ros4
+!        = 4 :    Rodas3
+!        = 5 :    Rodas4
+!
+!    ICNTRL(4)  -> maximum number of integration steps
+!        For ICNTRL(4) =0) the default value of 100000 is used
+!
+!    RCNTRL(1)  -> Hmin,lower bound for the integration step size
+!          It is strongly recommended to keep Hmin = ZERO
+!    RCNTRL(2)  -> Hmax,upper bound for the integration step size
+!    RCNTRL(3)  -> Hstart,starting value for the integration step size
+!
+!    RCNTRL(4)  -> FacMin,lower bound on step decrease factor (default=0.2)
+!    RCNTRL(5)  -> FacMax,upper bound on step increase factor (default=6)
+!    RCNTRL(6)  -> FacRej,step decrease factor after multiple rejections
+!                          (default=0.1)
+!    RCNTRL(7)  -> FacSafe,by which the new step is slightly smaller
+!         than the predicted value  (default=0.9)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!
+!
+!    OUTPUT ARGUMENTS:
+!    -----------------
+!
+!    T           -> T value for which the solution has been computed
+!                     (after successful return T=Tend).
+!
+!    Y(N)        -> Numerical solution at T
+!
+!    IDID        -> Reports on successfulness upon return:
+!                    = 1 for success
+!                    < 0 for error (value equals error code)
+!
+!    ISTATUS(1)  -> No. of function calls
+!    ISTATUS(2)  -> No. of jacobian calls
+!    ISTATUS(3)  -> No. of steps
+!    ISTATUS(4)  -> No. of accepted steps
+!    ISTATUS(5)  -> No. of rejected steps (except at very beginning)
+!    ISTATUS(6)  -> No. of LU decompositions
+!    ISTATUS(7)  -> No. of forward/backward substitutions
+!    ISTATUS(8)  -> No. of singular matrix decompositions
+!
+!    RSTATUS(1)  -> Texit,the time corresponding to the
+!                     computed Y upon return
+!    RSTATUS(2)  -> Hexit,last accepted step before exit
+!    RSTATUS(3)  -> Hnew,last predicted step (not yet taken)
+!                   For multiple restarts,use Hnew as Hstart 
+!                     in the subsequent run
+!
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+!~~~>  arguments
+   INTEGER,      INTENT(IN)  :: n
+   REAL(kind=dp), INTENT(INOUT):: y(n)
+   REAL(kind=dp), INTENT(IN)  :: tstart, tend
+   REAL(kind=dp), INTENT(IN)  :: abstol(n), reltol(n)
+   INTEGER,      INTENT(IN)  :: icntrl(20)
+   REAL(kind=dp), INTENT(IN)  :: rcntrl(20)
+   INTEGER,      INTENT(INOUT):: istatus(20)
+   REAL(kind=dp), INTENT(INOUT):: rstatus(20)
+   INTEGER, INTENT(OUT) :: ierr
+!~~~>  PARAMETERs of the rosenbrock method, up to 6 stages
+   INTEGER ::  ros_s, rosmethod
+   INTEGER, PARAMETER :: rs2=1, rs3=2, rs4=3, rd3=4, rd4=5, rg3=6
+   REAL(kind=dp):: ros_a(15), ros_c(15), ros_m(6), ros_e(6), &
+                    ros_alpha(6), ros_gamma(6), ros_elo
+   LOGICAL :: ros_newf(6)
+   CHARACTER(len=12):: ros_name
+!~~~>  local variables
+   REAL(kind=dp):: roundoff, facmin, facmax, facrej, facsafe
+   REAL(kind=dp):: hmin, hmax, hstart
+   REAL(kind=dp):: texit
+   INTEGER       :: i, uplimtol, max_no_steps
+   LOGICAL       :: autonomous, vectortol
+!~~~>   PARAMETERs
+   REAL(kind=dp), PARAMETER :: zero = 0.0_dp, one  = 1.0_dp
+   REAL(kind=dp), PARAMETER :: deltamin = 1.0e-5_dp
+
+!~~~>  initialize statistics
+   istatus(1:8) = 0
+   rstatus(1:3) = zero
+
+!~~~>  autonomous or time dependent ode. default is time dependent.
+   autonomous = .not.(icntrl(1) == 0)
+
+!~~~>  for scalar tolerances (icntrl(2).ne.0) the code uses abstol(1)and reltol(1)
+!   For Vector tolerances (ICNTRL(2) == 0) the code uses AbsTol(1:N) and RelTol(1:N)
+   IF (icntrl(2) == 0)THEN
+      vectortol = .TRUE.
+      uplimtol  = n
+   ELSE
+      vectortol = .FALSE.
+      uplimtol  = 1
+   ENDIF
+
+!~~~>   initialize the particular rosenbrock method selected
+   select CASE (icntrl(3))
+     CASE (1)
+       CALL ros2
+     CASE (2)
+       CALL ros3
+     CASE (3)
+       CALL ros4
+     CASE (0, 4)
+       CALL rodas3
+     CASE (5)
+       CALL rodas4
+     CASE (6)
+       CALL rang3
+     CASE default
+       PRINT *,'Unknown Rosenbrock method: ICNTRL(3) =',ICNTRL(3) 
+       CALL ros_errormsg(- 2, tstart, zero, ierr)
+       RETURN
+   END select
+
+!~~~>   the maximum number of steps admitted
+   IF (icntrl(4) == 0)THEN
+      max_no_steps = 200000
+   ELSEIF (icntrl(4)> 0)THEN
+      max_no_steps=icntrl(4)
+   ELSE
+      PRINT *,'User-selected max no. of steps: ICNTRL(4) =',ICNTRL(4)
+      CALL ros_errormsg(- 1, tstart, zero, ierr)
+      RETURN
+   ENDIF
+
+!~~~>  unit roundoff (1+ roundoff>1)
+   roundoff = epsilon(one)
+
+!~~~>  lower bound on the step size: (positive value)
+   IF (rcntrl(1) == zero)THEN
+      hmin = zero
+   ELSEIF (rcntrl(1)> zero)THEN
+      hmin = rcntrl(1)
+   ELSE
+      PRINT *,'User-selected Hmin: RCNTRL(1) =',RCNTRL(1)
+      CALL ros_errormsg(- 3, tstart, zero, ierr)
+      RETURN
+   ENDIF
+!~~~>  upper bound on the step size: (positive value)
+   IF (rcntrl(2) == zero)THEN
+      hmax = abs(tend-tstart)
+   ELSEIF (rcntrl(2)> zero)THEN
+      hmax = min(abs(rcntrl(2)), abs(tend-tstart))
+   ELSE
+      PRINT *,'User-selected Hmax: RCNTRL(2) =',RCNTRL(2)
+      CALL ros_errormsg(- 3, tstart, zero, ierr)
+      RETURN
+   ENDIF
+!~~~>  starting step size: (positive value)
+   IF (rcntrl(3) == zero)THEN
+      hstart = max(hmin, deltamin)
+   ELSEIF (rcntrl(3)> zero)THEN
+      hstart = min(abs(rcntrl(3)), abs(tend-tstart))
+   ELSE
+      PRINT *,'User-selected Hstart: RCNTRL(3) =',RCNTRL(3)
+      CALL ros_errormsg(- 3, tstart, zero, ierr)
+      RETURN
+   ENDIF
+!~~~>  step size can be changed s.t.  facmin < hnew/hold < facmax
+   IF (rcntrl(4) == zero)THEN
+      facmin = 0.2_dp
+   ELSEIF (rcntrl(4)> zero)THEN
+      facmin = rcntrl(4)
+   ELSE
+      PRINT *,'User-selected FacMin: RCNTRL(4) =',RCNTRL(4)
+      CALL ros_errormsg(- 4, tstart, zero, ierr)
+      RETURN
+   ENDIF
+   IF (rcntrl(5) == zero)THEN
+      facmax = 6.0_dp
+   ELSEIF (rcntrl(5)> zero)THEN
+      facmax = rcntrl(5)
+   ELSE
+      PRINT *,'User-selected FacMax: RCNTRL(5) =',RCNTRL(5)
+      CALL ros_errormsg(- 4, tstart, zero, ierr)
+      RETURN
+   ENDIF
+!~~~>   facrej: factor to decrease step after 2 succesive rejections
+   IF (rcntrl(6) == zero)THEN
+      facrej = 0.1_dp
+   ELSEIF (rcntrl(6)> zero)THEN
+      facrej = rcntrl(6)
+   ELSE
+      PRINT *,'User-selected FacRej: RCNTRL(6) =',RCNTRL(6)
+      CALL ros_errormsg(- 4, tstart, zero, ierr)
+      RETURN
+   ENDIF
+!~~~>   facsafe: safety factor in the computation of new step size
+   IF (rcntrl(7) == zero)THEN
+      facsafe = 0.9_dp
+   ELSEIF (rcntrl(7)> zero)THEN
+      facsafe = rcntrl(7)
+   ELSE
+      PRINT *,'User-selected FacSafe: RCNTRL(7) =',RCNTRL(7)
+      CALL ros_errormsg(- 4, tstart, zero, ierr)
+      RETURN
+   ENDIF
+!~~~>  check IF tolerances are reasonable
+    DO i=1, uplimtol
+      IF ((abstol(i)<= zero).or. (reltol(i)<= 10.0_dp* roundoff)&
+         .or. (reltol(i)>= 1.0_dp))THEN
+        PRINT *,' AbsTol(',i,') = ',AbsTol(i)
+        PRINT *,' RelTol(',i,') = ',RelTol(i)
+        CALL ros_errormsg(- 5, tstart, zero, ierr)
+        RETURN
+      ENDIF
+    ENDDO
+
+
+!~~~>  CALL rosenbrock method
+   CALL ros_integrator(y, tstart, tend, texit,  &
+        abstol, reltol,                         &
+!  Integration parameters
+        autonomous, vectortol, max_no_steps,    &
+        roundoff, hmin, hmax, hstart,           &
+        facmin, facmax, facrej, facsafe,        &
+!  Error indicator
+        ierr)
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+CONTAINS !  SUBROUTINEs internal to rosenbrock
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+   
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
+ SUBROUTINE ros_errormsg(code, t, h, ierr)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
+!    Handles all error messages
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~   
+  
+   REAL(kind=dp), INTENT(IN):: t, h
+   INTEGER, INTENT(IN) :: code
+   INTEGER, INTENT(OUT):: ierr
+   
+   ierr = code
+   print * , &
+     'Forced exit from Rosenbrock due to the following error:' 
+     
+   select CASE (code)
+    CASE (- 1)  
+      PRINT *,'--> Improper value for maximal no of steps'
+    CASE (- 2)  
+      PRINT *,'--> Selected Rosenbrock method not implemented'
+    CASE (- 3)  
+      PRINT *,'--> Hmin/Hmax/Hstart must be positive'
+    CASE (- 4)  
+      PRINT *,'--> FacMin/FacMax/FacRej must be positive'
+    CASE (- 5)
+      PRINT *,'--> Improper tolerance values'
+    CASE (- 6)
+      PRINT *,'--> No of steps exceeds maximum bound'
+    CASE (- 7)
+      PRINT *,'--> Step size too small: T + 10*H = T',&
+            ' or H < Roundoff'
+    CASE (- 8)  
+      PRINT *,'--> Matrix is repeatedly singular'
+    CASE default
+      PRINT *,'Unknown Error code: ',Code
+   END select
+   
+   print * , "t=", t, "and h=", h
+     
+ END SUBROUTINE ros_errormsg
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ SUBROUTINE ros_integrator (y, tstart, tend, t, &
+        abstol, reltol,                         &
+!~~~> integration PARAMETERs
+        autonomous, vectortol, max_no_steps,    &
+        roundoff, hmin, hmax, hstart,           &
+        facmin, facmax, facrej, facsafe,        &
+!~~~> error indicator
+        ierr)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!   Template for the implementation of a generic Rosenbrock method
+!      defined by ros_S (no of stages)
+!      and its coefficients ros_{A,C,M,E,Alpha,Gamma}
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+!~~~> input: the initial condition at tstart; output: the solution at t
+   REAL(kind=dp), INTENT(INOUT):: y(n)
+!~~~> input: integration interval
+   REAL(kind=dp), INTENT(IN):: tstart, tend
+!~~~> output: time at which the solution is RETURNed (t=tendIF success)
+   REAL(kind=dp), INTENT(OUT)::  t
+!~~~> input: tolerances
+   REAL(kind=dp), INTENT(IN)::  abstol(n), reltol(n)
+!~~~> input: integration PARAMETERs
+   LOGICAL, INTENT(IN):: autonomous, vectortol
+   REAL(kind=dp), INTENT(IN):: hstart, hmin, hmax
+   INTEGER, INTENT(IN):: max_no_steps
+   REAL(kind=dp), INTENT(IN):: roundoff, facmin, facmax, facrej, facsafe
+!~~~> output: error indicator
+   INTEGER, INTENT(OUT):: ierr
+! ~~~~ Local variables
+   REAL(kind=dp):: ynew(n), fcn0(n), fcn(n)
+   REAL(kind=dp):: k(n* ros_s), dfdt(n)
+#ifdef full_algebra    
+   REAL(kind=dp):: jac0(n, n), ghimj(n, n)
+#else
+   REAL(kind=dp):: jac0(lu_nonzero), ghimj(lu_nonzero)
+#endif
+   REAL(kind=dp):: h, hnew, hc, hg, fac, tau
+   REAL(kind=dp):: err, yerr(n)
+   INTEGER :: pivot(n), direction, ioffset, j, istage
+   LOGICAL :: rejectlasth, rejectmoreh, singular
+!~~~>  local PARAMETERs
+   REAL(kind=dp), PARAMETER :: zero = 0.0_dp, one  = 1.0_dp
+   REAL(kind=dp), PARAMETER :: deltamin = 1.0e-5_dp
+!~~~>  locally called FUNCTIONs
+!    REAL(kind=dp) WLAMCH
+!    EXTERNAL WLAMCH
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+!~~~>  initial preparations
+   t = tstart
+   rstatus(nhexit) = zero
+   h = min( max(abs(hmin), abs(hstart)), abs(hmax))
+   IF (abs(h)<= 10.0_dp* roundoff)h = deltamin
+
+   IF (tend  >=  tstart)THEN
+     direction = + 1
+   ELSE
+     direction = - 1
+   ENDIF
+   h = direction* h
+
+   rejectlasth=.FALSE.
+   rejectmoreh=.FALSE.
+
+!~~~> time loop begins below
+
+timeloop: DO WHILE((direction > 0).and.((t- tend) + roundoff <= zero)&
+       .or. (direction < 0).and.((tend-t) + roundoff <= zero))
+
+   IF (istatus(nstp)> max_no_steps)THEN  ! too many steps
+      CALL ros_errormsg(- 6, t, h, ierr)
+      RETURN
+   ENDIF
+   IF (((t+ 0.1_dp* h) == t).or.(h <= roundoff))THEN  ! step size too small
+      CALL ros_errormsg(- 7, t, h, ierr)
+      RETURN
+   ENDIF
+
+!~~~>  limit h IF necessary to avoid going beyond tend
+   h = min(h, abs(tend-t))
+
+!~~~>   compute the FUNCTION at current time
+   CALL funtemplate(t, y, fcn0)
+   istatus(nfun) = istatus(nfun) + 1
+
+!~~~>  compute the FUNCTION derivative with respect to t
+   IF (.not.autonomous)THEN
+      CALL ros_funtimederivative(t, roundoff, y, &
+                fcn0, dfdt)
+   ENDIF
+
+!~~~>   compute the jacobian at current time
+   CALL jactemplate(t, y, jac0)
+   istatus(njac) = istatus(njac) + 1
+
+!~~~>  repeat step calculation until current step accepted
+untilaccepted: do
+
+   CALL ros_preparematrix(h, direction, ros_gamma(1), &
+          jac0, ghimj, pivot, singular)
+   IF (singular)THEN ! more than 5 consecutive failed decompositions
+       CALL ros_errormsg(- 8, t, h, ierr)
+       RETURN
+   ENDIF
+
+!~~~>   compute the stages
+stage: DO istage = 1, ros_s
+
+      ! current istage offset. current istage vector is k(ioffset+ 1:ioffset+ n)
+       ioffset = n* (istage-1)
+
+      ! for the 1st istage the FUNCTION has been computed previously
+       IF (istage == 1)THEN
+         !slim: CALL wcopy(n, fcn0, 1, fcn, 1)
+       fcn(1:n) = fcn0(1:n)
+      ! istage>1 and a new FUNCTION evaluation is needed at the current istage
+       ELSEIF(ros_newf(istage))THEN
+         !slim: CALL wcopy(n, y, 1, ynew, 1)
+       ynew(1:n) = y(1:n)
+         DO j = 1, istage-1
+           CALL waxpy(n, ros_a((istage-1) * (istage-2) /2+ j), &
+            k(n* (j- 1) + 1), 1, ynew, 1)
+         ENDDO
+         tau = t + ros_alpha(istage) * direction* h
+         CALL funtemplate(tau, ynew, fcn)
+         istatus(nfun) = istatus(nfun) + 1
+       ENDIF ! IF istage == 1 ELSEIF ros_newf(istage)
+       !slim: CALL wcopy(n, fcn, 1, k(ioffset+ 1), 1)
+       k(ioffset+ 1:ioffset+ n) = fcn(1:n)
+       DO j = 1, istage-1
+         hc = ros_c((istage-1) * (istage-2) /2+ j) /(direction* h)
+         CALL waxpy(n, hc, k(n* (j- 1) + 1), 1, k(ioffset+ 1), 1)
+       ENDDO
+       IF ((.not. autonomous).and.(ros_gamma(istage).ne.zero))THEN
+         hg = direction* h* ros_gamma(istage)
+         CALL waxpy(n, hg, dfdt, 1, k(ioffset+ 1), 1)
+       ENDIF
+       CALL ros_solve(ghimj, pivot, k(ioffset+ 1))
+
+   END DO stage
+
+
+!~~~>  compute the new solution
+   !slim: CALL wcopy(n, y, 1, ynew, 1)
+   ynew(1:n) = y(1:n)
+   DO j=1, ros_s
+         CALL waxpy(n, ros_m(j), k(n* (j- 1) + 1), 1, ynew, 1)
+   ENDDO
+
+!~~~>  compute the error estimation
+   !slim: CALL wscal(n, zero, yerr, 1)
+   yerr(1:n) = zero
+   DO j=1, ros_s
+        CALL waxpy(n, ros_e(j), k(n* (j- 1) + 1), 1, yerr, 1)
+   ENDDO
+   err = ros_errornorm(y, ynew, yerr, abstol, reltol, vectortol)
+
+!~~~> new step size is bounded by facmin <= hnew/h <= facmax
+   fac  = min(facmax, max(facmin, facsafe/err** (one/ros_elo)))
+   hnew = h* fac
+
+!~~~>  check the error magnitude and adjust step size
+   istatus(nstp) = istatus(nstp) + 1
+   IF ((err <= one).or.(h <= hmin))THEN  !~~~> accept step
+      istatus(nacc) = istatus(nacc) + 1
+      !slim: CALL wcopy(n, ynew, 1, y, 1)
+      y(1:n) = ynew(1:n)
+      t = t + direction* h
+      hnew = max(hmin, min(hnew, hmax))
+      IF (rejectlasth)THEN  ! no step size increase after a rejected step
+         hnew = min(hnew, h)
+      ENDIF
+      rstatus(nhexit) = h
+      rstatus(nhnew) = hnew
+      rstatus(ntexit) = t
+      rejectlasth = .FALSE.
+      rejectmoreh = .FALSE.
+      h = hnew
+      exit untilaccepted ! exit the loop: WHILE step not accepted
+   ELSE           !~~~> reject step
+      IF (rejectmoreh)THEN
+         hnew = h* facrej
+      ENDIF
+      rejectmoreh = rejectlasth
+      rejectlasth = .TRUE.
+      h = hnew
+      IF (istatus(nacc)>= 1) istatus(nrej) = istatus(nrej) + 1
+   ENDIF ! err <= 1
+
+   END DO untilaccepted
+
+   END DO timeloop
+
+!~~~> succesful exit
+   ierr = 1  !~~~> the integration was successful
+
+  END SUBROUTINE ros_integrator
+
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  REAL(kind=dp)FUNCTION ros_errornorm(y, ynew, yerr, &
+                               abstol, reltol, vectortol)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~> computes the "scaled norm" of the error vector yerr
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+! Input arguments
+   REAL(kind=dp), INTENT(IN):: y(n), ynew(n), &
+          yerr(n), abstol(n), reltol(n)
+   LOGICAL, INTENT(IN)::  vectortol
+! Local variables
+   REAL(kind=dp):: err, scale, ymax
+   INTEGER  :: i
+   REAL(kind=dp), PARAMETER :: zero = 0.0_dp
+
+   err = zero
+   DO i=1, n
+     ymax = max(abs(y(i)), abs(ynew(i)))
+     IF (vectortol)THEN
+       scale = abstol(i) + reltol(i) * ymax
+     ELSE
+       scale = abstol(1) + reltol(1) * ymax
+     ENDIF
+     err = err+ (yerr(i) /scale) ** 2
+   ENDDO
+   err  = sqrt(err/n)
+
+   ros_errornorm = max(err, 1.0d-10)
+
+  END FUNCTION ros_errornorm
+
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  SUBROUTINE ros_funtimederivative(t, roundoff, y, &
+                fcn0, dfdt)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~> the time partial derivative of the FUNCTION by finite differences
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+!~~~> input arguments
+   REAL(kind=dp), INTENT(IN):: t, roundoff, y(n), fcn0(n)
+!~~~> output arguments
+   REAL(kind=dp), INTENT(OUT):: dfdt(n)
+!~~~> local variables
+   REAL(kind=dp):: delta
+   REAL(kind=dp), PARAMETER :: one = 1.0_dp, deltamin = 1.0e-6_dp
+
+   delta = sqrt(roundoff) * max(deltamin, abs(t))
+   CALL funtemplate(t+ delta, y, dfdt)
+   istatus(nfun) = istatus(nfun) + 1
+   CALL waxpy(n, (- one), fcn0, 1, dfdt, 1)
+   CALL wscal(n, (one/delta), dfdt, 1)
+
+  END SUBROUTINE ros_funtimederivative
+
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  SUBROUTINE ros_preparematrix(h, direction, gam, &
+             jac0, ghimj, pivot, singular)
+! --- --- --- --- --- --- --- --- --- --- --- --- ---
+!  Prepares the LHS matrix for stage calculations
+!  1.  Construct Ghimj = 1/(H*ham) - Jac0
+!      "(Gamma H) Inverse Minus Jacobian"
+!  2.  Repeat LU decomposition of Ghimj until successful.
+!       -half the step size if LU decomposition fails and retry
+!       -exit after 5 consecutive fails
+! --- --- --- --- --- --- --- --- --- --- --- --- ---
+
+!~~~> input arguments
+#ifdef full_algebra    
+   REAL(kind=dp), INTENT(IN)::  jac0(n, n)
+#else
+   REAL(kind=dp), INTENT(IN)::  jac0(lu_nonzero)
+#endif   
+   REAL(kind=dp), INTENT(IN)::  gam
+   INTEGER, INTENT(IN)::  direction
+!~~~> output arguments
+#ifdef full_algebra    
+   REAL(kind=dp), INTENT(OUT):: ghimj(n, n)
+#else
+   REAL(kind=dp), INTENT(OUT):: ghimj(lu_nonzero)
+#endif   
+   LOGICAL, INTENT(OUT)::  singular
+   INTEGER, INTENT(OUT)::  pivot(n)
+!~~~> inout arguments
+   REAL(kind=dp), INTENT(INOUT):: h   ! step size is decreased when lu fails
+!~~~> local variables
+   INTEGER  :: i, ising, nconsecutive
+   REAL(kind=dp):: ghinv
+   REAL(kind=dp), PARAMETER :: one  = 1.0_dp, half = 0.5_dp
+
+   nconsecutive = 0
+   singular = .TRUE.
+
+   DO WHILE (singular)
+
+!~~~>    construct ghimj = 1/(h* gam) - jac0
+#ifdef full_algebra    
+     !slim: CALL wcopy(n* n, jac0, 1, ghimj, 1)
+     !slim: CALL wscal(n* n, (- one), ghimj, 1)
+     ghimj = - jac0
+     ghinv = one/(direction* h* gam)
+     DO i=1, n
+       ghimj(i, i) = ghimj(i, i) + ghinv
+     ENDDO
+#else
+     !slim: CALL wcopy(lu_nonzero, jac0, 1, ghimj, 1)
+     !slim: CALL wscal(lu_nonzero, (- one), ghimj, 1)
+     ghimj(1:lu_nonzero) = - jac0(1:lu_nonzero)
+     ghinv = one/(direction* h* gam)
+     DO i=1, n
+       ghimj(lu_diag(i)) = ghimj(lu_diag(i)) + ghinv
+     ENDDO
+#endif   
+!~~~>    compute lu decomposition
+     CALL ros_decomp( ghimj, pivot, ising)
+     IF (ising == 0)THEN
+!~~~>    IF successful done
+        singular = .FALSE.
+     ELSE ! ising .ne. 0
+!~~~>    IF unsuccessful half the step size; IF 5 consecutive fails THEN RETURN
+        istatus(nsng) = istatus(nsng) + 1
+        nconsecutive = nconsecutive+1
+        singular = .TRUE.
+        PRINT*,'Warning: LU Decomposition returned ISING = ',ISING
+        IF (nconsecutive <= 5)THEN ! less than 5 consecutive failed decompositions
+           h = h* half
+        ELSE  ! more than 5 consecutive failed decompositions
+           RETURN
+        ENDIF  ! nconsecutive
+      ENDIF    ! ising
+
+   END DO ! WHILE singular
+
+  END SUBROUTINE ros_preparematrix
+
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  SUBROUTINE ros_decomp( a, pivot, ising)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!  Template for the LU decomposition
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~> inout variables
+#ifdef full_algebra    
+   REAL(kind=dp), INTENT(INOUT):: a(n, n)
+#else   
+   REAL(kind=dp), INTENT(INOUT):: a(lu_nonzero)
+#endif
+!~~~> output variables
+   INTEGER, INTENT(OUT):: pivot(n), ising
+
+#ifdef full_algebra    
+   CALL  dgetrf( n, n, a, n, pivot, ising)
+#else   
+   CALL kppdecomp(a, ising)
+   pivot(1) = 1
+#endif
+   istatus(ndec) = istatus(ndec) + 1
+
+  END SUBROUTINE ros_decomp
+
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  SUBROUTINE ros_solve( a, pivot, b)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!  Template for the forward/backward substitution (using pre-computed LU decomposition)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~> input variables
+#ifdef full_algebra    
+   REAL(kind=dp), INTENT(IN):: a(n, n)
+   INTEGER :: ising
+#else   
+   REAL(kind=dp), INTENT(IN):: a(lu_nonzero)
+#endif
+   INTEGER, INTENT(IN):: pivot(n)
+!~~~> inout variables
+   REAL(kind=dp), INTENT(INOUT):: b(n)
+
+#ifdef full_algebra    
+   CALL  DGETRS( 'N',N ,1,A,N,Pivot,b,N,ISING)
+   IF (info < 0)THEN
+      print* , "error in dgetrs. ising=", ising
+   ENDIF  
+#else   
+   CALL kppsolve( a, b)
+#endif
+
+   istatus(nsol) = istatus(nsol) + 1
+
+  END SUBROUTINE ros_solve
+
+
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  SUBROUTINE ros2
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! --- AN L-STABLE METHOD,2 stages,order 2
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+   double precision g
+
+    g = 1.0_dp + 1.0_dp/sqrt(2.0_dp)
+    rosmethod = rs2
+!~~~> name of the method
+    ros_Name = 'ROS-2'
+!~~~> number of stages
+    ros_s = 2
+
+!~~~> the coefficient matrices a and c are strictly lower triangular.
+!   The lower triangular (subdiagonal) elements are stored in row-wise order:
+!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
+!   The general mapping formula is:
+!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
+!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
+
+    ros_a(1) = (1.0_dp) /g
+    ros_c(1) = (- 2.0_dp) /g
+!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
+!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
+    ros_newf(1) = .TRUE.
+    ros_newf(2) = .TRUE.
+!~~~> m_i = coefficients for new step solution
+    ros_m(1) = (3.0_dp) /(2.0_dp* g)
+    ros_m(2) = (1.0_dp) /(2.0_dp* g)
+! E_i = Coefficients for error estimator
+    ros_e(1) = 1.0_dp/(2.0_dp* g)
+    ros_e(2) = 1.0_dp/(2.0_dp* g)
+!~~~> ros_elo = estimator of local order - the minimum between the
+!    main and the embedded scheme orders plus one
+    ros_elo = 2.0_dp
+!~~~> y_stage_i ~ y( t + h* alpha_i)
+    ros_alpha(1) = 0.0_dp
+    ros_alpha(2) = 1.0_dp
+!~~~> gamma_i = \sum_j  gamma_{i, j}
+    ros_gamma(1) = g
+    ros_gamma(2) = -g
+
+ END SUBROUTINE ros2
+
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  SUBROUTINE ros3
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! --- AN L-STABLE METHOD,3 stages,order 3,2 function evaluations
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+   rosmethod = rs3
+!~~~> name of the method
+   ros_Name = 'ROS-3'
+!~~~> number of stages
+   ros_s = 3
+
+!~~~> the coefficient matrices a and c are strictly lower triangular.
+!   The lower triangular (subdiagonal) elements are stored in row-wise order:
+!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
+!   The general mapping formula is:
+!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
+!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
+
+   ros_a(1) = 1.0_dp
+   ros_a(2) = 1.0_dp
+   ros_a(3) = 0.0_dp
+
+   ros_c(1) = - 0.10156171083877702091975600115545e+01_dp
+   ros_c(2) =  0.40759956452537699824805835358067e+01_dp
+   ros_c(3) =  0.92076794298330791242156818474003e+01_dp
+!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
+!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
+   ros_newf(1) = .TRUE.
+   ros_newf(2) = .TRUE.
+   ros_newf(3) = .FALSE.
+!~~~> m_i = coefficients for new step solution
+   ros_m(1) =  0.1e+01_dp
+   ros_m(2) =  0.61697947043828245592553615689730e+01_dp
+   ros_m(3) = - 0.42772256543218573326238373806514_dp
+! E_i = Coefficients for error estimator
+   ros_e(1) =  0.5_dp
+   ros_e(2) = - 0.29079558716805469821718236208017e+01_dp
+   ros_e(3) =  0.22354069897811569627360909276199_dp
+!~~~> ros_elo = estimator of local order - the minimum between the
+!    main and the embedded scheme orders plus 1
+   ros_elo = 3.0_dp
+!~~~> y_stage_i ~ y( t + h* alpha_i)
+   ros_alpha(1) = 0.0_dp
+   ros_alpha(2) = 0.43586652150845899941601945119356_dp
+   ros_alpha(3) = 0.43586652150845899941601945119356_dp
+!~~~> gamma_i = \sum_j  gamma_{i, j}
+   ros_gamma(1) = 0.43586652150845899941601945119356_dp
+   ros_gamma(2) = 0.24291996454816804366592249683314_dp
+   ros_gamma(3) = 0.21851380027664058511513169485832e+01_dp
+
+  END SUBROUTINE ros3
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  SUBROUTINE ros4
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!     L-STABLE ROSENBROCK METHOD OF ORDER 4,WITH 4 STAGES
+!     L-STABLE EMBEDDED ROSENBROCK METHOD OF ORDER 3
+!
+!      E. HAIRER AND G. WANNER,SOLVING ORDINARY DIFFERENTIAL
+!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
+!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
+!      SPRINGER-VERLAG (1990)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+   rosmethod = rs4
+!~~~> name of the method
+   ros_Name = 'ROS-4'
+!~~~> number of stages
+   ros_s = 4
+
+!~~~> the coefficient matrices a and c are strictly lower triangular.
+!   The lower triangular (subdiagonal) elements are stored in row-wise order:
+!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
+!   The general mapping formula is:
+!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
+!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
+
+   ros_a(1) = 0.2000000000000000e+01_dp
+   ros_a(2) = 0.1867943637803922e+01_dp
+   ros_a(3) = 0.2344449711399156_dp
+   ros_a(4) = ros_a(2)
+   ros_a(5) = ros_a(3)
+   ros_a(6) = 0.0_dp
+
+   ros_c(1) = -0.7137615036412310e+01_dp
+   ros_c(2) = 0.2580708087951457e+01_dp
+   ros_c(3) = 0.6515950076447975_dp
+   ros_c(4) = -0.2137148994382534e+01_dp
+   ros_c(5) = -0.3214669691237626_dp
+   ros_c(6) = -0.6949742501781779_dp
+!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
+!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
+   ros_newf(1) = .TRUE.
+   ros_newf(2) = .TRUE.
+   ros_newf(3) = .TRUE.
+   ros_newf(4) = .FALSE.
+!~~~> m_i = coefficients for new step solution
+   ros_m(1) = 0.2255570073418735e+01_dp
+   ros_m(2) = 0.2870493262186792_dp
+   ros_m(3) = 0.4353179431840180_dp
+   ros_m(4) = 0.1093502252409163e+01_dp
+!~~~> e_i  = coefficients for error estimator
+   ros_e(1) = -0.2815431932141155_dp
+   ros_e(2) = -0.7276199124938920e-01_dp
+   ros_e(3) = -0.1082196201495311_dp
+   ros_e(4) = -0.1093502252409163e+01_dp
+!~~~> ros_elo  = estimator of local order - the minimum between the
+!    main and the embedded scheme orders plus 1
+   ros_elo  = 4.0_dp
+!~~~> y_stage_i ~ y( t + h* alpha_i)
+   ros_alpha(1) = 0.0_dp
+   ros_alpha(2) = 0.1145640000000000e+01_dp
+   ros_alpha(3) = 0.6552168638155900_dp
+   ros_alpha(4) = ros_alpha(3)
+!~~~> gamma_i = \sum_j  gamma_{i, j}
+   ros_gamma(1) = 0.5728200000000000_dp
+   ros_gamma(2) = -0.1769193891319233e+01_dp
+   ros_gamma(3) = 0.7592633437920482_dp
+   ros_gamma(4) = -0.1049021087100450_dp
+
+  END SUBROUTINE ros4
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  SUBROUTINE rodas3
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! --- A STIFFLY-STABLE METHOD,4 stages,order 3
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+   rosmethod = rd3
+!~~~> name of the method
+   ros_Name = 'RODAS-3'
+!~~~> number of stages
+   ros_s = 4
+
+!~~~> the coefficient matrices a and c are strictly lower triangular.
+!   The lower triangular (subdiagonal) elements are stored in row-wise order:
+!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
+!   The general mapping formula is:
+!       A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
+!       C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
+
+   ros_a(1) = 0.0_dp
+   ros_a(2) = 2.0_dp
+   ros_a(3) = 0.0_dp
+   ros_a(4) = 2.0_dp
+   ros_a(5) = 0.0_dp
+   ros_a(6) = 1.0_dp
+
+   ros_c(1) = 4.0_dp
+   ros_c(2) = 1.0_dp
+   ros_c(3) = -1.0_dp
+   ros_c(4) = 1.0_dp
+   ros_c(5) = -1.0_dp
+   ros_c(6) = -(8.0_dp/3.0_dp)
+
+!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
+!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
+   ros_newf(1) = .TRUE.
+   ros_newf(2) = .FALSE.
+   ros_newf(3) = .TRUE.
+   ros_newf(4) = .TRUE.
+!~~~> m_i = coefficients for new step solution
+   ros_m(1) = 2.0_dp
+   ros_m(2) = 0.0_dp
+   ros_m(3) = 1.0_dp
+   ros_m(4) = 1.0_dp
+!~~~> e_i  = coefficients for error estimator
+   ros_e(1) = 0.0_dp
+   ros_e(2) = 0.0_dp
+   ros_e(3) = 0.0_dp
+   ros_e(4) = 1.0_dp
+!~~~> ros_elo  = estimator of local order - the minimum between the
+!    main and the embedded scheme orders plus 1
+   ros_elo  = 3.0_dp
+!~~~> y_stage_i ~ y( t + h* alpha_i)
+   ros_alpha(1) = 0.0_dp
+   ros_alpha(2) = 0.0_dp
+   ros_alpha(3) = 1.0_dp
+   ros_alpha(4) = 1.0_dp
+!~~~> gamma_i = \sum_j  gamma_{i, j}
+   ros_gamma(1) = 0.5_dp
+   ros_gamma(2) = 1.5_dp
+   ros_gamma(3) = 0.0_dp
+   ros_gamma(4) = 0.0_dp
+
+  END SUBROUTINE rodas3
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  SUBROUTINE rodas4
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!     STIFFLY-STABLE ROSENBROCK METHOD OF ORDER 4,WITH 6 STAGES
+!
+!      E. HAIRER AND G. WANNER,SOLVING ORDINARY DIFFERENTIAL
+!      EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
+!      SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
+!      SPRINGER-VERLAG (1996)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+    rosmethod = rd4
+!~~~> name of the method
+    ros_Name = 'RODAS-4'
+!~~~> number of stages
+    ros_s = 6
+
+!~~~> y_stage_i ~ y( t + h* alpha_i)
+    ros_alpha(1) = 0.000_dp
+    ros_alpha(2) = 0.386_dp
+    ros_alpha(3) = 0.210_dp
+    ros_alpha(4) = 0.630_dp
+    ros_alpha(5) = 1.000_dp
+    ros_alpha(6) = 1.000_dp
+
+!~~~> gamma_i = \sum_j  gamma_{i, j}
+    ros_gamma(1) = 0.2500000000000000_dp
+    ros_gamma(2) = -0.1043000000000000_dp
+    ros_gamma(3) = 0.1035000000000000_dp
+    ros_gamma(4) = -0.3620000000000023e-01_dp
+    ros_gamma(5) = 0.0_dp
+    ros_gamma(6) = 0.0_dp
+
+!~~~> the coefficient matrices a and c are strictly lower triangular.
+!   The lower triangular (subdiagonal) elements are stored in row-wise order:
+!   A(2,1) = ros_A(1),A(3,1) =ros_A(2),A(3,2) =ros_A(3),etc.
+!   The general mapping formula is:  A(i,j) = ros_A( (i-1)*(i-2)/2 + j)
+!                  C(i,j) = ros_C( (i-1)*(i-2)/2 + j)
+
+    ros_a(1) = 0.1544000000000000e+01_dp
+    ros_a(2) = 0.9466785280815826_dp
+    ros_a(3) = 0.2557011698983284_dp
+    ros_a(4) = 0.3314825187068521e+01_dp
+    ros_a(5) = 0.2896124015972201e+01_dp
+    ros_a(6) = 0.9986419139977817_dp
+    ros_a(7) = 0.1221224509226641e+01_dp
+    ros_a(8) = 0.6019134481288629e+01_dp
+    ros_a(9) = 0.1253708332932087e+02_dp
+    ros_a(10) = -0.6878860361058950_dp
+    ros_a(11) = ros_a(7)
+    ros_a(12) = ros_a(8)
+    ros_a(13) = ros_a(9)
+    ros_a(14) = ros_a(10)
+    ros_a(15) = 1.0_dp
+
+    ros_c(1) = -0.5668800000000000e+01_dp
+    ros_c(2) = -0.2430093356833875e+01_dp
+    ros_c(3) = -0.2063599157091915_dp
+    ros_c(4) = -0.1073529058151375_dp
+    ros_c(5) = -0.9594562251023355e+01_dp
+    ros_c(6) = -0.2047028614809616e+02_dp
+    ros_c(7) = 0.7496443313967647e+01_dp
+    ros_c(8) = -0.1024680431464352e+02_dp
+    ros_c(9) = -0.3399990352819905e+02_dp
+    ros_c(10) = 0.1170890893206160e+02_dp
+    ros_c(11) = 0.8083246795921522e+01_dp
+    ros_c(12) = -0.7981132988064893e+01_dp
+    ros_c(13) = -0.3152159432874371e+02_dp
+    ros_c(14) = 0.1631930543123136e+02_dp
+    ros_c(15) = -0.6058818238834054e+01_dp
+
+!~~~> m_i = coefficients for new step solution
+    ros_m(1) = ros_a(7)
+    ros_m(2) = ros_a(8)
+    ros_m(3) = ros_a(9)
+    ros_m(4) = ros_a(10)
+    ros_m(5) = 1.0_dp
+    ros_m(6) = 1.0_dp
+
+!~~~> e_i  = coefficients for error estimator
+    ros_e(1) = 0.0_dp
+    ros_e(2) = 0.0_dp
+    ros_e(3) = 0.0_dp
+    ros_e(4) = 0.0_dp
+    ros_e(5) = 0.0_dp
+    ros_e(6) = 1.0_dp
+
+!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
+!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
+    ros_newf(1) = .TRUE.
+    ros_newf(2) = .TRUE.
+    ros_newf(3) = .TRUE.
+    ros_newf(4) = .TRUE.
+    ros_newf(5) = .TRUE.
+    ros_newf(6) = .TRUE.
+
+!~~~> ros_elo  = estimator of local order - the minimum between the
+!        main and the embedded scheme orders plus 1
+    ros_elo = 4.0_dp
+
+  END SUBROUTINE rodas4
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  SUBROUTINE rang3
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! STIFFLY-STABLE W METHOD OF ORDER 3,WITH 4 STAGES
+!
+! J. RANG and L. ANGERMANN
+! NEW ROSENBROCK W-METHODS OF ORDER 3
+! FOR PARTIAL DIFFERENTIAL ALGEBRAIC
+!        EQUATIONS OF INDEX 1                                             
+! BIT Numerical Mathematics (2005) 45: 761-787
+!  DOI: 10.1007/s10543-005-0035-y
+! Table 4.1-4.2
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+
+    rosmethod = rg3
+!~~~> name of the method
+    ros_Name = 'RANG-3'
+!~~~> number of stages
+    ros_s = 4
+
+    ros_a(1) = 5.09052051067020d+00;
+    ros_a(2) = 5.09052051067020d+00;
+    ros_a(3) = 0.0d0;
+    ros_a(4) = 4.97628111010787d+00;
+    ros_a(5) = 2.77268164715849d-02;
+    ros_a(6) = 2.29428036027904d-01;
+
+    ros_c(1) = - 1.16790812312283d+01;
+    ros_c(2) = - 1.64057326467367d+01;
+    ros_c(3) = - 2.77268164715850d-01;
+    ros_c(4) = - 8.38103960500476d+00;
+    ros_c(5) = - 8.48328409199343d-01;
+    ros_c(6) =  2.87009860433106d-01;
+
+    ros_m(1) =  5.22582761233094d+00;
+    ros_m(2) = - 5.56971148154165d-01;
+    ros_m(3) =  3.57979469353645d-01;
+    ros_m(4) =  1.72337398521064d+00;
+
+    ros_e(1) = - 5.16845212784040d+00;
+    ros_e(2) = - 1.26351942603842d+00;
+    ros_e(3) = - 1.11022302462516d-16;
+    ros_e(4) =  2.22044604925031d-16;
+
+    ros_alpha(1) = 0.0d00;
+    ros_alpha(2) = 2.21878746765329d+00;
+    ros_alpha(3) = 2.21878746765329d+00;
+    ros_alpha(4) = 1.55392337535788d+00;
+
+    ros_gamma(1) =  4.35866521508459d-01;
+    ros_gamma(2) = - 1.78292094614483d+00;
+    ros_gamma(3) = - 2.46541900496934d+00;
+    ros_gamma(4) = - 8.05529997906370d-01;
+
+
+!~~~> does the stage i require a new FUNCTION evaluation (ros_newf(i) =true)
+!   or does it re-use the function evaluation from stage i-1 (ros_NewF(i) =FALSE)
+    ros_newf(1) = .TRUE.
+    ros_newf(2) = .TRUE.
+    ros_newf(3) = .TRUE.
+    ros_newf(4) = .TRUE.
+
+!~~~> ros_elo  = estimator of local order - the minimum between the
+!        main and the embedded scheme orders plus 1
+    ros_elo = 3.0_dp
+
+  END SUBROUTINE rang3
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!   End of the set of internal Rosenbrock subroutines
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+END SUBROUTINE rosenbrock
+ 
+SUBROUTINE funtemplate( t, y, ydot)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!  Template for the ODE function call.
+!  Updates the rate coefficients (and possibly the fixed species) at each call
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~> input variables
+   REAL(kind=dp):: t, y(nvar)
+!~~~> output variables
+   REAL(kind=dp):: ydot(nvar)
+!~~~> local variables
+   REAL(kind=dp):: told
+
+   told = time
+   time = t
+   CALL fun( y, fix, rconst, ydot)
+   time = told
+
+END SUBROUTINE funtemplate
+ 
+SUBROUTINE jactemplate( t, y, jcb)
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!  Template for the ODE Jacobian call.
+!  Updates the rate coefficients (and possibly the fixed species) at each call
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!~~~> input variables
+    REAL(kind=dp):: t, y(nvar)
+!~~~> output variables
+#ifdef full_algebra    
+    REAL(kind=dp):: jv(lu_nonzero), jcb(nvar, nvar)
+#else
+    REAL(kind=dp):: jcb(lu_nonzero)
+#endif   
+!~~~> local variables
+    REAL(kind=dp):: told
+#ifdef full_algebra    
+    INTEGER :: i, j
+#endif   
+
+    told = time
+    time = t
+#ifdef full_algebra    
+    CALL jac_sp(y, fix, rconst, jv)
+    DO j=1, nvar
+      DO i=1, nvar
+         jcb(i, j) = 0.0_dp
+      ENDDO
+    ENDDO
+    DO i=1, lu_nonzero
+       jcb(lu_irow(i), lu_icol(i)) = jv(i)
+    ENDDO
+#else
+    CALL jac_sp( y, fix, rconst, jcb)
+#endif   
+    time = told
+
+END SUBROUTINE jactemplate
+ 
+  SUBROUTINE kppdecomp( jvs, ier)                                    
+  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
+  !        sparse lu factorization                                    
+  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
+  ! loop expansion generated by kp4                                   
+                                                                      
+    INTEGER  :: ier                                                   
+    REAL(kind=dp):: jvs(lu_nonzero), w(nvar), a             
+    INTEGER  :: k, kk, j, jj                                          
+                                                                      
+    a = 0.                                                            
+    ier = 0                                                           
+                                                                      
+!   i = 1
+!   i = 2
+!   i = 3
+!   i = 4
+!   i = 5
+    jvs(12) = (jvs(12)) / jvs(7) 
+    jvs(14) = jvs(14) - jvs(8) * jvs(12)
+!   i = 6
+    jvs(16) = (jvs(16)) / jvs(7) 
+    jvs(17) = (jvs(17)) / jvs(9) 
+    a = 0.0; a = a  - jvs(10) * jvs(17)
+    jvs(18) = (jvs(18) + a) / jvs(13) 
+    jvs(19) = jvs(19) - jvs(8) * jvs(16) - jvs(14) * jvs(18)
+    jvs(22) = jvs(22) - jvs(11) * jvs(17) - jvs(15) * jvs(18)
+!   i = 7
+    jvs(23) = (jvs(23)) / jvs(9) 
+    a = 0.0; a = a  - jvs(10) * jvs(23)
+    jvs(24) = (jvs(24) + a) / jvs(13) 
+    a = 0.0; a = a  - jvs(14) * jvs(24)
+    jvs(25) = (jvs(25) + a) / jvs(19) 
+    jvs(26) = jvs(26) - jvs(20) * jvs(25)
+    jvs(27) = jvs(27) - jvs(21) * jvs(25)
+    jvs(28) = jvs(28) - jvs(11) * jvs(23) - jvs(15) * jvs(24) - jvs(22) * jvs(25)
+!   i = 8
+    jvs(29) = (jvs(29)) / jvs(26) 
+    jvs(30) = jvs(30) - jvs(27) * jvs(29)
+    jvs(31) = jvs(31) - jvs(28) * jvs(29)
+!   i = 9
+    jvs(32) = (jvs(32)) / jvs(9) 
+    a = 0.0; a = a  - jvs(10) * jvs(32)
+    jvs(33) = (jvs(33) + a) / jvs(13) 
+    a = 0.0; a = a  - jvs(14) * jvs(33)
+    jvs(34) = (jvs(34) + a) / jvs(19) 
+    a = 0.0; a = a  - jvs(20) * jvs(34)
+    jvs(35) = (jvs(35) + a) / jvs(26) 
+    a = 0.0; a = a  - jvs(21) * jvs(34) - jvs(27) * jvs(35)
+    jvs(36) = (jvs(36) + a) / jvs(30) 
+    jvs(37) = jvs(37) - jvs(11) * jvs(32) - jvs(15) * jvs(33) - jvs(22) * jvs(34) - jvs(28) * jvs(35)&
+          - jvs(31) * jvs(36)
+    RETURN                                                            
+                                                                      
+  END SUBROUTINE kppdecomp                                            
+ 
+SUBROUTINE chem_gasphase_integrate (time_step_len, conc, tempi, qvapi, fakti, photo, ierrf, xnacc, xnrej, istatus, l_debug, pe, &
+                     icntrl_i, rcntrl_i) 
+                                                                    
+  IMPLICIT NONE                                                     
+                                                                    
+  REAL(dp), INTENT(IN)                  :: time_step_len           
+  REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: conc                    
+  REAL(dp), DIMENSION(:, :), INTENT(IN)   :: photo                   
+  REAL(dp), DIMENSION(:), INTENT(IN)     :: tempi                   
+  REAL(dp), DIMENSION(:), INTENT(IN)     :: qvapi                   
+  REAL(dp), DIMENSION(:), INTENT(IN)     :: fakti                   
+  INTEGER, INTENT(OUT), OPTIONAL        :: ierrf(:)               
+  INTEGER, INTENT(OUT), OPTIONAL        :: xnacc(:)               
+  INTEGER, INTENT(OUT), OPTIONAL        :: xnrej(:)               
+  INTEGER, INTENT(INOUT), OPTIONAL      :: istatus(:)             
+  INTEGER, INTENT(IN), OPTIONAL         :: pe                      
+  LOGICAL, INTENT(IN), OPTIONAL         :: l_debug                 
+  INTEGER, DIMENSION(nkppctrl), INTENT(IN), OPTIONAL  :: icntrl_i         
+  REAL(dp), DIMENSION(nkppctrl), INTENT(IN), OPTIONAL  :: rcntrl_i         
+                                                                    
+  INTEGER                                 :: k   ! loop variable     
+  REAL(dp)                               :: dt                      
+  INTEGER, DIMENSION(20)                :: istatus_u               
+  INTEGER                                :: ierr_u                  
+  INTEGER                                :: istatf                  
+  INTEGER                                :: vl_dim_lo               
+                                                                    
+                                                                    
+  IF (PRESENT (istatus)) istatus = 0                             
+  IF (PRESENT (icntrl_i)) icntrl  = icntrl_i                      
+  IF (PRESENT (rcntrl_i)) rcntrl  = rcntrl_i                      
+                                                                    
+  vl_glo = size(tempi, 1)                                           
+                                                                    
+  vl_dim_lo = vl_dim                                                
+  DO k=1, vl_glo, vl_dim_lo                                           
+    is = k                                                          
+    ie = min(k+ vl_dim_lo-1, vl_glo)                                 
+    vl = ie-is+ 1                                                    
+                                                                    
+    c(:) = conc(is, :)                                            
+                                                                    
+    temp = tempi(is)                                             
+                                                                    
+    qvap = qvapi(is)                                             
+                                                                    
+    fakt = fakti(is)                                             
+                                                                    
+    CALL initialize                                                 
+                                                                    
+    phot(:) = photo(is, :)                                            
+                                                                    
+    CALL update_rconst                                              
+                                                                    
+    dt = time_step_len                                              
+                                                                    
+    ! integrate from t=0 to t=dt                                    
+    CALL integrate(0._dp, dt, icntrl, rcntrl, istatus_u = istatus_u, ierr_u=ierr_u)
+                                                                    
+                                                                    
+   IF (PRESENT(l_debug) .AND. PRESENT(pe)) THEN                       
+      IF (l_debug) CALL error_output(conc(is, :), ierr_u, pe)          
+   ENDIF                                                              
+                                                                      
+    conc(is, :) = c(:)                                                
+                                                                    
+    ! RETURN diagnostic information                                 
+                                                                    
+    IF (PRESENT(ierrf))   ierrf(is) = ierr_u                      
+    IF (PRESENT(xnacc))   xnacc(is) = istatus_u(4)               
+    IF (PRESENT(xnrej))   xnrej(is) = istatus_u(5)               
+                                                                    
+    IF (PRESENT (istatus)) THEN                                   
+      istatus(1:8) = istatus(1:8) + istatus_u(1:8)                
+    ENDIF                                                          
+                                                                    
+  END DO                                                            
+ 
+                                                                    
+! Deallocate input arrays                                           
+                                                                    
+                                                                    
+  data_loaded = .FALSE.                                             
+                                                                    
+  RETURN                                                            
+END SUBROUTINE chem_gasphase_integrate                                        
+
+END MODULE chem_gasphase_mod
+ 
Index: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_simple/chem_gasphase_mod.kpp
===================================================================
--- palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_simple/chem_gasphase_mod.kpp	(revision 2696)
+++ palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_simple/chem_gasphase_mod.kpp	(revision 3298)
@@ -1,3 +1,7 @@
 //chem_gasphase_mod.kpp
+//
+// Current revision
+// forkel 20.04.2018: added F90_INIT for initialization of fix compounds
+//
 //
 //Former revisions
@@ -11,7 +15,4 @@
 #HESSIAN    on
 #STOICMAT   on
-#INLINE F90_GLOBAL
-  !   Declaration of global variable declarations for photolysis will come from INLINE F90_DATA
-#ENDINLINE
 //
 // *******************************************************************************************
Index: palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_simple/simple.eqn
===================================================================
--- palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_simple/simple.eqn	(revision 2696)
+++ palm/trunk/UTIL/chemistry/gasphase_preproc/mechanisms/def_simple/simple.eqn	(revision 3298)
@@ -1,10 +1,11 @@
 {simple.eqn
-Former revisions
+Current revision
 ----------------
- $Id: smog.eqn 2459 2017-09-13 14:10:33Z forkel $
+ 201711   Created simple.eqn with 6 equations forkel
+ 20180316 Added equation no. 7                forkel
 }
 #EQUATIONS
 
-{ very simple mechanism - do NOT show to a chemist! }
+{ home made very simple mechanism - please do not show to a chemist! }
 
   { 1.}  NO2 + hv = NO + O3 :                   phot(j_no2)  ;
@@ -14,2 +15,3 @@
   { 5.}  RO2 + NO = NO2 + RCHO + HO2 :       arr2(4.2E-12_dp, -180.0_dp, temp);
   { 6.}  HO2 + NO = NO2 + OH :               arr2(3.7E-12_dp, -240.0_dp, temp) ;
+  { 7.}  NO2 + OH = HNO3     :               arr2(1.15E-11_dp, 0.0_dp, temp) ;
