MODULE chem_gasphase_mod
!------------------------------------------------------------------------------!
!
! ******Module chem_gasphase_mod is automatically generated by kpp4palm ******
!
! *********Please do NOT change this Code *********
!
!------------------------------------------------------------------------------!
! 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 .
!
! Copyright 1997-2018 Leibniz Universitaet Hannover
!--------------------------------------------------------------------------------!
!
! Current revisions:
! ------------------
!
!
! Former revisions:
! -----------------
! $Id: module_header 2460 2017-09-13 14:47:48Z forkel $
! Removed preprocessor directive __chem
!
! 2460 2017-09-13 14:47:48Z forkel
!
!
! 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 :: phot
PUBLIC :: rconst
PUBLIC :: nvar
PUBLIC :: nphot
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 = 0
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 : Fri Dec 8 11:54:15 2017
! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
! Equation file : chem_gasphase_mod.kpp
! Output root filename : chem_gasphase_mod
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! NSPEC - Number of chemical species
INTEGER,PARAMETER :: nspec = 2
! NVAR - Number of Variable species
INTEGER,PARAMETER :: nvar = 2
! NVARACT - Number of Active species
INTEGER,PARAMETER :: nvaract = 2
! NFIX - Number of Fixed species
INTEGER,PARAMETER :: nfix = 1
! NREACT - Number of reactions
INTEGER,PARAMETER :: nreact = 2
! NVARST - Starting of variables in conc. vect.
INTEGER,PARAMETER :: nvarst = 1
! NFIXST - Starting of fixed in conc. vect.
INTEGER,PARAMETER :: nfixst = 3
! NONZERO - Number of nonzero entries in Jacobian
INTEGER,PARAMETER :: nonzero = 2
! LU_NONZERO - Number of nonzero entries in LU factoriz. of Jacobian
INTEGER,PARAMETER :: lu_nonzero = 2
! CNVAR - (NVAR+1) Number of elements in compressed row format
INTEGER,PARAMETER :: cnvar = 3
! CNEQN - (NREACT+1) Number stoicm elements in compressed col format
INTEGER,PARAMETER :: cneqn = 3
! NHESS - Length of Sparse Hessian
INTEGER,PARAMETER :: nhess = 1
! NLOOKAT - Number of species to look at
INTEGER,PARAMETER :: nlookat = 0
! NMONITOR - Number of species to monitor
INTEGER,PARAMETER :: nmonitor = 0
! 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_pm10 = 1
INTEGER,PARAMETER,PUBLIC :: ind_pm25 = 2
! Index declaration for fixed species in C
! C(ind_spc)
! Index declaration for fixed species in FIX
! FIX(indf_spc) = C(ind_spc) = C(NVAR+indf_spc)
! NJVRP - Length of sparse Jacobian JVRP
INTEGER,PARAMETER :: njvrp = 2
! NSTOICM - Length of Sparse Stoichiometric Matrix
INTEGER,PARAMETER :: nstoicm = 1
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! 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 : Fri Dec 8 11:54:15 2017
! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/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))
! RCONST - Rate constants (global)
REAL(kind=dp):: rconst(nreact)
! TIME - Current integration time
REAL(kind=dp):: time
! SUN - Sunlight intensity between [0,1]
REAL(kind=dp):: sun
! TEMP - Temperature
REAL(dp),dimension(:),allocatable :: temp
! RTOLS - (scalar) Relative tolerance
REAL(kind=dp):: rtols
! TSTART - Integration start time
REAL(kind=dp):: tstart
! TEND - Integration end time
REAL(kind=dp):: tend
! DT - Integration step
REAL(kind=dp):: dt
! 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
! STEPMAX - Upper bound for integration step
REAL(kind=dp):: stepmax
! CFACTOR - Conversion factor for concentration units
REAL(kind=dp):: cfactor
! DDMTYPE - DDM sensitivity w.r.t.: 0=init.val.,1=params
INTEGER :: ddmtype
! INLINED global variable declarations
! declaration of global variable declarations for photolysis will come from
! 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 : Fri Dec 8 11:54:15 2017
! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
! Equation file : chem_gasphase_mod.kpp
! Output root filename : chem_gasphase_mod
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Sparse Jacobian Data
INTEGER,PARAMETER,DIMENSION(2):: lu_irow = (/ &
1, 2 /)
INTEGER,PARAMETER,DIMENSION(2):: lu_icol = (/ &
1, 2 /)
INTEGER,PARAMETER,DIMENSION(3):: lu_crow = (/ &
1, 2, 3 /)
INTEGER,PARAMETER,DIMENSION(3):: lu_diag = (/ &
1, 2, 3 /)
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! 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 : Fri Dec 8 11:54:15 2017
! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/GASPHASE_PREPROC/tmp_kpp4palm
! Equation file : chem_gasphase_mod.kpp
! Output root filename : chem_gasphase_mod
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
CHARACTER(len=15),PARAMETER,DIMENSION(2):: spc_names = (/ &
'PM10 ','PM25 ' /)
INTEGER,DIMENSION(1):: lookat
INTEGER,DIMENSION(1):: monitor
CHARACTER(len=15),DIMENSION(1):: smass
CHARACTER(len=100),PARAMETER,DIMENSION(2):: eqn_names = (/ &
'PM10 --> PM10 ',&
'PM25 --> PM25 ' /)
! 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 = 1
! phot photolysis frequencies
REAL(kind=dp):: phot(nphot)
INTEGER,PARAMETER,PUBLIC :: j_no2 = 1
CHARACTER(len=15),PARAMETER,DIMENSION(nphot):: phot_names = (/ &
'J_NO2 '/)
! 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 : Fri Dec 8 11:54:15 2017
! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/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 : Fri Dec 8 11:54:15 2017
! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/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 : Fri Dec 8 11:54:15 2017
! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/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 : Fri Dec 8 11:54:15 2017
! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/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 : Fri Dec 8 11:54:15 2017
! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/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 : Fri Dec 8 11:54:15 2017
! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/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 : Fri Dec 8 11:54:15 2017
! Working directory : /home/forkel-r/palmstuff/work/chemistry20171117/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 kppdecomp
MODULE PROCEDURE kppdecomp
END INTERFACE kppdecomp
INTERFACE wlamch
MODULE PROCEDURE wlamch
END INTERFACE wlamch
INTERFACE wlamch_add
MODULE PROCEDURE wlamch_add
END INTERFACE wlamch_add
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 not working INTERFACE wcopy
!interface not working MODULE PROCEDURE wcopy
!interface not working END INTERFACE wcopy
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 chem_gasphase_integrate
MODULE PROCEDURE chem_gasphase_integrate
END INTERFACE chem_gasphase_integrate
INTERFACE fill_temp
MODULE PROCEDURE fill_temp
END INTERFACE fill_temp
PUBLIC fill_temp
CONTAINS
SUBROUTINE initialize()
INTEGER :: i
REAL(kind=dp):: x
cfactor = 1.000000e+00_dp
x = (0.)*cfactor
DO i = 1,nvar
var(i) = x
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
! Aggregate function
vdot(1) = 0
vdot(2) = 0
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(2) = x(2)/ jvs(2)
x(1) = x(1)/ jvs(1)
END SUBROUTINE kppsolve
SUBROUTINE kppdecomp( jvs,ier)
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Sparse LU factorization
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
INTEGER :: ier
REAL(kind=dp):: jvs(lu_nonzero),w(nvar),a
INTEGER :: k,kk,j,jj
a = 0. ! mz_rs_20050606
ier = 0
DO k=1,nvar
! mz_rs_20050606: don't check if real value == 0
! IF(jvs( lu_diag(k)).eq. 0.)THEN
IF(abs(jvs(lu_diag(k)))< tiny(a))THEN
ier = k
RETURN
ENDIF
DO kk = lu_crow(k),lu_crow(k+ 1)- 1
w( lu_icol(kk)) = jvs(kk)
ENDDO
DO kk = lu_crow(k),lu_diag(k)- 1
j = lu_icol(kk)
a = - w(j)/ jvs( lu_diag(j))
w(j) = - a
DO jj = lu_diag(j)+ 1,lu_crow(j+ 1)- 1
w( lu_icol(jj)) = w( lu_icol(jj))+ a*jvs(jj)
ENDDO
ENDDO
DO kk = lu_crow(k),lu_crow(k+ 1)- 1
jvs(kk) = w( lu_icol(kk))
ENDDO
ENDDO
END SUBROUTINE kppdecomp
REAL(kind=dp)FUNCTION wlamch( c)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! returns epsilon machine
! after LAPACK
! replace this by the function from the optimized LAPACK implementation:
! CALL SLAMCH('E') or CALL DLAMCH('E')
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! USE chem_gasphase_mod_Precision
CHARACTER :: c
INTEGER :: i
REAL(kind=dp),SAVE :: eps
REAL(kind=dp) :: suma
REAL(kind=dp),PARAMETER :: one=1.0_dp, half=0.5_dp
LOGICAL,SAVE :: first=.true.
IF (first)THEN
first = .false.
eps = half**(16)
DO i = 17,80
eps = eps*half
CALL wlamch_add(one,eps,suma)
IF (suma.le.one)goto 10
ENDDO
PRINT*,'ERROR IN WLAMCH. EPS < ',Eps
RETURN
10 eps = eps*2
i = i- 1
ENDIF
wlamch = eps
END FUNCTION wlamch
SUBROUTINE wlamch_add( a,b,suma)
! USE chem_gasphase_mod_Precision
REAL(kind=dp)a,b,suma
suma = a + b
END SUBROUTINE wlamch_add
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(2)
! B(1) = dA(1)/dV(1)
b(1) = rct(1)
! B(2) = dA(2)/dV(2)
b(2) = rct(2)
! Construct the Jacobian terms from B's
! JVS(1) = Jac_FULL(1,1)
jvs(1) = 0
! JVS(2) = Jac_FULL(2,2)
jvs(2) = 0
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 :: j,k
k = is
! Begin INLINED RCONST
! End INLINED RCONST
rconst(1) = (1.0_dp)
rconst(2) = (1.0_dp)
END SUBROUTINE update_rconst
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,iou,modstr)
! i/o
INTEGER, INTENT(out):: status
INTEGER, INTENT(in) :: iou ! LOGICAL i/o unit
CHARACTER(len=*),INTENT(in) :: modstr ! read .nml
! 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 wcopy(n,x,incx,y,incy)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! copies a vector,x,to a vector,y: y <- x
! only for incX=incY=1
! after BLAS
! replace this by the function from the optimized BLAS implementation:
! CALL SCOPY(N,X,1,Y,1) or CALL DCOPY(N,X,1,Y,1)
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! USE chem_gasphase_mod_Precision
INTEGER :: i,incx,incy,m,mp1,n
REAL(kind=dp):: x(n),y(n)
IF (n.le.0)RETURN
m = mod(n,8)
IF( m .ne. 0)THEN
DO i = 1,m
y(i) = x(i)
ENDDO
IF( n .lt. 8)RETURN
ENDIF
mp1 = m+ 1
DO i = mp1,n,8
y(i) = x(i)
y(i + 1) = x(i + 1)
y(i + 2) = x(i + 2)
y(i + 3) = x(i + 3)
y(i + 4) = x(i + 4)
y(i + 5) = x(i + 5)
y(i + 6) = x(i + 6)
y(i + 7) = x(i + 7)
ENDDO
END SUBROUTINE wcopy
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 = WLAMCH('E')
!~~~> 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 chem_gasphase_integrate (time_step_len,conc,tempk,photo,ierrf,xnacc,xnrej,istatus,l_debug,pe)
IMPLICIT NONE
REAL(dp), INTENT(IN) :: time_step_len
REAL(dp), DIMENSION(:,:), INTENT(INOUT) :: conc
REAL(dp), DIMENSION(:,:), INTENT(INOUT) :: photo
REAL(dp), DIMENSION(:), INTENT(IN) :: tempk
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 :: k ! loop variable
REAL(dp) :: dt
INTEGER, DIMENSION(20) :: istatus_u
INTEGER :: ierr_u
if (present (istatus)) istatus = 0
DO k=1,vl_glo,vl_dim
is = k
ie = min(k+vl_dim-1,vl_glo)
vl = ie-is+1
c(:) = conc(is,:)
temp = tempk(is)
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
if (ALLOCATED(temp)) DEALLOCATE(temp)
data_loaded = .false.
RETURN
END SUBROUTINE chem_gasphase_integrate
SUBROUTINE fill_temp(status,array)
INTEGER, INTENT(OUT) :: status
REAL(dp), INTENT(IN), DIMENSION(:) :: array
status = 0
IF (.not. ALLOCATED(temp)) &
ALLOCATE(temp(size(array)))
IF (data_loaded .AND. (vl_glo /= size(array,1))) THEN
status = 1
RETURN
END IF
vl_glo = size(array,1)
temp = array
data_loaded = .TRUE.
RETURN
END SUBROUTINE fill_temp
END MODULE chem_gasphase_mod