From 03623d29519079927c250763c73b6d8c02b0bec1 Mon Sep 17 00:00:00 2001 From: Nick Wogan Date: Thu, 9 Mar 2023 08:49:30 -0800 Subject: [PATCH] fast_arbitrary_rate was too slow. Also made it non-constant --- CMakeLists.txt | 2 +- photochem/cython/PhotochemVars.pyx | 8 ++++++++ photochem/cython/PhotochemVars_pxd.pxd | 5 ++++- photochem/fortran/PhotochemVars_wrapper.f90 | 16 ++++++++++++++++ src/atmosphere/photochem_atmosphere_rhs.f90 | 4 ++-- src/atmosphere/photochem_atmosphere_utils.f90 | 3 +-- .../photochem_evoatmosphere_rhs.f90 | 4 ++-- src/photochem_const.f90 | 1 - src/photochem_types.f90 | 2 ++ 9 files changed, 36 insertions(+), 9 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 1c38e77..27f0906 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,5 @@ cmake_minimum_required(VERSION "3.14") -project(Photochem LANGUAGES Fortran C VERSION "0.3.14") +project(Photochem LANGUAGES Fortran C VERSION "0.3.15") include(FortranCInterface) FortranCInterface_VERIFY() diff --git a/photochem/cython/PhotochemVars.pyx b/photochem/cython/PhotochemVars.pyx index cf38a90..d20a548 100644 --- a/photochem/cython/PhotochemVars.pyx +++ b/photochem/cython/PhotochemVars.pyx @@ -144,5 +144,13 @@ cdef class PhotochemVars: def __set__(self, int val): var_pxd.photochemvars_verbose_set(&self._ptr, &val) + property fast_arbitrary_rate: + def __get__(self): + cdef double val + var_pxd.photochemvars_fast_arbitrary_rate_get(&self._ptr, &val) + return val + def __set__(self, double val): + var_pxd.photochemvars_fast_arbitrary_rate_set(&self._ptr, &val) + \ No newline at end of file diff --git a/photochem/cython/PhotochemVars_pxd.pxd b/photochem/cython/PhotochemVars_pxd.pxd index 702a520..1707f31 100644 --- a/photochem/cython/PhotochemVars_pxd.pxd +++ b/photochem/cython/PhotochemVars_pxd.pxd @@ -50,4 +50,7 @@ cdef extern void photochemvars_equilibrium_time_get(void *ptr, double *val) cdef extern void photochemvars_equilibrium_time_set(void *ptr, double *val) cdef extern void photochemvars_verbose_get(void *ptr, int *val) -cdef extern void photochemvars_verbose_set(void *ptr, int *val) \ No newline at end of file +cdef extern void photochemvars_verbose_set(void *ptr, int *val) + +cdef extern void photochemvars_fast_arbitrary_rate_get(void *ptr, double *val) +cdef extern void photochemvars_fast_arbitrary_rate_set(void *ptr, double *val) \ No newline at end of file diff --git a/photochem/fortran/PhotochemVars_wrapper.f90 b/photochem/fortran/PhotochemVars_wrapper.f90 index 5ce8e3b..e9fdf15 100644 --- a/photochem/fortran/PhotochemVars_wrapper.f90 +++ b/photochem/fortran/PhotochemVars_wrapper.f90 @@ -267,4 +267,20 @@ subroutine photochemvars_verbose_set(ptr, val) bind(c) call c_f_pointer(ptr, var) var%verbose = val end subroutine + + subroutine photochemvars_fast_arbitrary_rate_get(ptr, val) bind(c) + type(c_ptr), intent(in) :: ptr + real(c_double), intent(out) :: val + type(PhotochemVars), pointer :: var + call c_f_pointer(ptr, var) + val = var%fast_arbitrary_rate + end subroutine + + subroutine photochemvars_fast_arbitrary_rate_set(ptr, val) bind(c) + type(c_ptr), intent(in) :: ptr + real(c_double), intent(in) :: val + type(PhotochemVars), pointer :: var + call c_f_pointer(ptr, var) + var%fast_arbitrary_rate = val + end subroutine \ No newline at end of file diff --git a/src/atmosphere/photochem_atmosphere_rhs.f90 b/src/atmosphere/photochem_atmosphere_rhs.f90 index 53edb6f..e3cd36d 100644 --- a/src/atmosphere/photochem_atmosphere_rhs.f90 +++ b/src/atmosphere/photochem_atmosphere_rhs.f90 @@ -18,7 +18,7 @@ subroutine dochem(self, neqs, nsp, np, nsl, nq, nz, trop_ind, nrT, usol, density use photochem_enum, only: CondensingParticle use photochem_common, only: chempl, chempl_sl use photochem_eqns, only: damp_condensation_rate - use photochem_const, only: N_avo, pi, small_real, T_crit_H2O, fast_arbitrary_rate + use photochem_const, only: N_avo, pi, small_real, T_crit_H2O class(Atmosphere), target, intent(in) :: self integer, intent(in) :: neqs, nsp, np, nsl, nq, nz, trop_ind, nrT @@ -84,7 +84,7 @@ subroutine dochem(self, neqs, nsp, np, nsl, nq, nz, trop_ind, nrT, usol, density if (dat%fix_water_in_trop) then do j = 1,trop_ind k = dat%LH2O + (j - 1) * dat%nq - rhs(k) = rhs(k) + fast_arbitrary_rate*(H2O_sat_mix(j)*H2O_rh(j) - usol(dat%LH2O,j)) + rhs(k) = rhs(k) + var%fast_arbitrary_rate*(H2O_sat_mix(j)*H2O_rh(j) - usol(dat%LH2O,j)) enddo endif ! H2O condensation diff --git a/src/atmosphere/photochem_atmosphere_utils.f90 b/src/atmosphere/photochem_atmosphere_utils.f90 index 24e900a..a399f23 100644 --- a/src/atmosphere/photochem_atmosphere_utils.f90 +++ b/src/atmosphere/photochem_atmosphere_utils.f90 @@ -100,7 +100,6 @@ module function atom_conservation(self, atom, err) result(con) use photochem_enum, only: VelocityDistributedFluxBC use photochem_eqns, only: damp_condensation_rate use photochem_types, only: AtomConservation - use photochem_const, only: fast_arbitrary_rate class(Atmosphere), target, intent(inout) :: self character(len=*), intent(in) :: atom character(:), allocatable, intent(out) :: err @@ -184,7 +183,7 @@ module function atom_conservation(self, atom, err) result(con) if (dat%fix_water_in_trop) then do j = 1,var%trop_ind - con_evap_rate = fast_arbitrary_rate*(wrk%H2O_sat_mix(j)*wrk%H2O_rh(j) - wrk%usol(dat%LH2O,j)) & + con_evap_rate = var%fast_arbitrary_rate*(wrk%H2O_sat_mix(j)*wrk%H2O_rh(j) - wrk%usol(dat%LH2O,j)) & *wrk%density(j)*var%dz(j)*dat%species_composition(kk,dat%LH2O) if (con_evap_rate > 0.0_dp) then con%in_other = con%in_other + con_evap_rate diff --git a/src/evoatmosphere/photochem_evoatmosphere_rhs.f90 b/src/evoatmosphere/photochem_evoatmosphere_rhs.f90 index 04e283a..89c25c6 100644 --- a/src/evoatmosphere/photochem_evoatmosphere_rhs.f90 +++ b/src/evoatmosphere/photochem_evoatmosphere_rhs.f90 @@ -23,7 +23,7 @@ subroutine dochem(self, usol, rx_rates, & use photochem_enum, only: CondensingParticle use photochem_common, only: chempl, chempl_sl use photochem_eqns, only: damp_condensation_rate - use photochem_const, only: N_avo, pi, small_real, T_crit_H2O, fast_arbitrary_rate + use photochem_const, only: N_avo, pi, small_real, T_crit_H2O class(EvoAtmosphere), target, intent(in) :: self real(dp), intent(in) :: usol(:,:) @@ -89,7 +89,7 @@ subroutine dochem(self, usol, rx_rates, & if (dat%fix_water_in_trop) then do j = 1,var%trop_ind k = dat%LH2O + (j - 1) * dat%nq - rhs(k) = rhs(k) + fast_arbitrary_rate*(density(j)*H2O_sat_mix(j)*H2O_rh(j) - usol(dat%LH2O,j)) + rhs(k) = rhs(k) + var%fast_arbitrary_rate*(density(j)*H2O_sat_mix(j)*H2O_rh(j) - usol(dat%LH2O,j)) enddo endif if (dat%water_cond) then diff --git a/src/photochem_const.f90 b/src/photochem_const.f90 index 783e618..80366e7 100644 --- a/src/photochem_const.f90 +++ b/src/photochem_const.f90 @@ -22,5 +22,4 @@ module photochem_const real(dp), parameter :: small_real = tiny(1.0_dp)**0.25_dp real(dp), parameter :: ln_small_real = log(small_real) - real(dp), parameter :: fast_arbitrary_rate = 1.0e-5_dp ! arbitrary rate that is fast (1/s) end module diff --git a/src/photochem_types.f90 b/src/photochem_types.f90 index efa6106..4124278 100644 --- a/src/photochem_types.f90 +++ b/src/photochem_types.f90 @@ -441,6 +441,8 @@ subroutine time_dependent_rate_fcn(tn, nz, rate) integer(c_int) :: max_order = 5 real(dp) :: epsj = 1.d-9 ! perturbation for jacobian calculation integer :: verbose = 1 ! 0 == no printing. 1 == some printing. 2 == bunch of printing. + !> arbitrary rate that is fast (1/s). Used for keeping H2O at saturation in troposphere + real(dp) :: fast_arbitrary_rate = 1.0e-2_dp end type type :: SundialsData