MODULE chem_gasphase_mod ! Mechanism: passive ! !------------------------------------------------------------------------------! ! ! ******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 . ! ! Copyright 1997-2018 Leibniz Universitaet Hannover !--------------------------------------------------------------------------------! ! ! Current revisions: ! ------------------ ! ! ! Former revisions: ! ----------------- ! $Id: chem_gasphase_mod.f90 3287 2018-09-28 10:19:58Z forkel $ ! forkel June 2018: qvap,fakt added ! forkel June 2018: reset case in Initialize,Integrate,Update_rconst ! ! ! 3287 2018-09-28 10:19:58Z 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:34:57 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 = 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 ! 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 : Tue Sep 25 18:34:57 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)) ! 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:34:57 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(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 : Tue Sep 25 18:34:57 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(2):: spc_names = (/ & 'PM10 ','PM25 ' /) 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 : Tue Sep 25 18:34:57 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:34:57 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:34:57 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:34:57 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:34:57 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:34:57 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:34:57 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 ! 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 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 :: k k = is ! Begin INLINED RCONST ! End INLINED RCONST rconst(1) = (1.0_dp) rconst(2) = (1.0_dp) 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 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