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

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