From ef95e44518db167fde92b14489cac3fb4f8dea56 Mon Sep 17 00:00:00 2001 From: Nick Wogan Date: Tue, 16 Apr 2024 12:04:11 -0700 Subject: [PATCH] ability to change planet mass radius for evoatmosphere --- CMakeLists.txt | 2 +- photochem/cython/EvoAtmosphere.pyx | 12 ++++++++++++ photochem/cython/EvoAtmosphere_pxd.pxd | 1 + photochem/cython/PhotochemData.pyx | 14 ++++++++++++++ photochem/cython/PhotochemData_pxd.pxd | 4 ++++ photochem/fortran/EvoAtmosphere_wrapper.f90 | 10 ++++++++++ photochem/fortran/PhotochemData_wrapper.f90 | 16 ++++++++++++++++ src/evoatmosphere/photochem_evoatmosphere.f90 | 7 +++++++ .../photochem_evoatmosphere_utils.f90 | 10 ++++++++++ 9 files changed, 75 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index ab3499a..b3d7599 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,6 +1,6 @@ cmake_minimum_required(VERSION "3.14") cmake_policy(SET CMP0148 OLD) -project(Photochem LANGUAGES Fortran C VERSION "0.5.3") +project(Photochem LANGUAGES Fortran C VERSION "0.5.4") include(FortranCInterface) FortranCInterface_VERIFY() diff --git a/photochem/cython/EvoAtmosphere.pyx b/photochem/cython/EvoAtmosphere.pyx index 53bbdc2..17d5843 100644 --- a/photochem/cython/EvoAtmosphere.pyx +++ b/photochem/cython/EvoAtmosphere.pyx @@ -385,6 +385,18 @@ cdef class EvoAtmosphere: if len(err.strip()) > 0: raise PhotoException(err.decode("utf-8").strip()) + def update_planet_mass_radius(self, double planet_mass, double planet_radius): + """Updates planet mass and radius. The routine recomputes gravity vs. altitude. + + Parameters + ---------- + planet_mass : float + In grams + planet_radius : float + In cm + """ + ea_pxd.evoatmosphere_update_planet_mass_radius_wrapper(&self._ptr, &planet_mass, &planet_radius) + def regrid_prep_atmosphere(self, ndarray[double, ndim=2] usol, double top_atmos): """This subroutine calculates re-grids the model so that the top of the model domain is at `top_atmos` the computes reaction rates, photolysis rates, etc. diff --git a/photochem/cython/EvoAtmosphere_pxd.pxd b/photochem/cython/EvoAtmosphere_pxd.pxd index 5906b0d..d3fc808 100644 --- a/photochem/cython/EvoAtmosphere_pxd.pxd +++ b/photochem/cython/EvoAtmosphere_pxd.pxd @@ -32,6 +32,7 @@ cdef extern void evoatmosphere_set_press_temp_edd_wrapper(void *ptr, int *P_dim1 bool *hydro_pressure, bool *hydro_pressure_present, char *err) cdef extern void evoatmosphere_update_vertical_grid_wrapper(void *ptr, double *toa_alt, bool *toa_alt_present, double *toa_pressure, bool *toa_pressure_present, char *err) +cdef extern void evoatmosphere_update_planet_mass_radius_wrapper(void *ptr, double *planet_mass, double *planet_radius) cdef extern void evoatmosphere_regrid_prep_atmosphere_wrapper(void *ptr, int *nq, int *nz, double *usol, double *top_atmos, char *err) cdef extern void evoatmosphere_evolve_wrapper(void *ptr, char *filename, diff --git a/photochem/cython/PhotochemData.pyx b/photochem/cython/PhotochemData.pyx index fbfe967..fd8eb26 100644 --- a/photochem/cython/PhotochemData.pyx +++ b/photochem/cython/PhotochemData.pyx @@ -69,6 +69,20 @@ cdef class PhotochemData: dat_pxd.photochemdata_nw_get(&self._ptr, &val) return val + property planet_mass: + "Planet mass (g)." + def __get__(self): + cdef double val + dat_pxd.photochemdata_planet_mass_get(&self._ptr, &val) + return val + + property planet_radius: + "Planet radius (cm)." + def __get__(self): + cdef double val + dat_pxd.photochemdata_planet_radius_get(&self._ptr, &val) + return val + property species_names: """List, shape (nsp+2). A list of the species in the model (particles and gases). The last two elements are 'hv' and 'M'. diff --git a/photochem/cython/PhotochemData_pxd.pxd b/photochem/cython/PhotochemData_pxd.pxd index da5335b..4c430d3 100644 --- a/photochem/cython/PhotochemData_pxd.pxd +++ b/photochem/cython/PhotochemData_pxd.pxd @@ -13,6 +13,10 @@ cdef extern void photochemdata_nsl_get(void *ptr, int *nq) cdef extern void photochemdata_nll_get(void *ptr, int *nq) cdef extern void photochemdata_nw_get(void *ptr, int *nq) +cdef extern void photochemdata_planet_mass_get(void *ptr, double *val) + +cdef extern void photochemdata_planet_radius_get(void *ptr, double *val) + cdef extern void photochemdata_species_names_get_size(void *ptr, int *dim1) cdef extern void photochemdata_species_names_get(void *ptr, int *dim1, char* species_names) diff --git a/photochem/fortran/EvoAtmosphere_wrapper.f90 b/photochem/fortran/EvoAtmosphere_wrapper.f90 index 25f2e35..c2c253b 100644 --- a/photochem/fortran/EvoAtmosphere_wrapper.f90 +++ b/photochem/fortran/EvoAtmosphere_wrapper.f90 @@ -334,6 +334,16 @@ subroutine evoatmosphere_update_vertical_grid_wrapper(ptr, TOA_alt, TOA_alt_pres end subroutine + subroutine evoatmosphere_update_planet_mass_radius_wrapper(ptr, planet_mass, planet_radius) bind(c) + type(c_ptr), intent(in) :: ptr + real(c_double), intent(in) :: planet_mass, planet_radius + type(EvoAtmosphere), pointer :: pc + + call c_f_pointer(ptr, pc) + call pc%update_planet_mass_radius(planet_mass, planet_radius) + + end subroutine + subroutine evoatmosphere_regrid_prep_atmosphere_wrapper(ptr, nq, nz, usol, top_atmos, err) bind(c) type(c_ptr), intent(in) :: ptr integer(c_int), intent(in) :: nq, nz diff --git a/photochem/fortran/PhotochemData_wrapper.f90 b/photochem/fortran/PhotochemData_wrapper.f90 index 6140958..7e4369d 100644 --- a/photochem/fortran/PhotochemData_wrapper.f90 +++ b/photochem/fortran/PhotochemData_wrapper.f90 @@ -76,6 +76,22 @@ subroutine photochemdata_nw_get(ptr, val) bind(c) call c_f_pointer(ptr, dat) val = dat%nw end subroutine + + subroutine photochemdata_planet_mass_get(ptr, val) bind(c) + type(c_ptr), intent(in) :: ptr + real(c_double), intent(out) :: val + type(PhotochemData), pointer :: dat + call c_f_pointer(ptr, dat) + val = dat%planet_mass + end subroutine + + subroutine photochemdata_planet_radius_get(ptr, val) bind(c) + type(c_ptr), intent(in) :: ptr + real(c_double), intent(out) :: val + type(PhotochemData), pointer :: dat + call c_f_pointer(ptr, dat) + val = dat%planet_radius + end subroutine subroutine photochemdata_species_names_get_size(ptr, dim1) bind(c) type(c_ptr), intent(in) :: ptr diff --git a/src/evoatmosphere/photochem_evoatmosphere.f90 b/src/evoatmosphere/photochem_evoatmosphere.f90 index efaec68..2a8c962 100644 --- a/src/evoatmosphere/photochem_evoatmosphere.f90 +++ b/src/evoatmosphere/photochem_evoatmosphere.f90 @@ -64,6 +64,7 @@ function temp_dependent_albedo_fcn(T_surf) result(albedo) procedure :: set_temperature procedure :: set_press_temp_edd procedure :: update_vertical_grid + procedure :: update_planet_mass_radius procedure :: rebin_update_vertical_grid procedure :: regrid_prep_atmosphere @@ -314,6 +315,12 @@ module subroutine update_vertical_grid(self, TOA_alt, TOA_pressure, err) character(:), allocatable, intent(out) :: err end subroutine + !> Updates planet mass and radius. The routine recomputes gravity vs. altitude. + module subroutine update_planet_mass_radius(self, planet_mass, planet_radius) + class(EvoAtmosphere), target, intent(inout) :: self + real(dp), intent(in) :: planet_mass, planet_radius + end subroutine + ! Both below are needed only for `evolve` routine and probably should not be used most of the time. module subroutine rebin_update_vertical_grid(self, usol_old, top_atmos, usol_new, err) diff --git a/src/evoatmosphere/photochem_evoatmosphere_utils.f90 b/src/evoatmosphere/photochem_evoatmosphere_utils.f90 index af6a800..85e26b6 100644 --- a/src/evoatmosphere/photochem_evoatmosphere_utils.f90 +++ b/src/evoatmosphere/photochem_evoatmosphere_utils.f90 @@ -793,6 +793,16 @@ module subroutine update_vertical_grid(self, TOA_alt, TOA_pressure, err) end subroutine + module subroutine update_planet_mass_radius(self, planet_mass, planet_radius) + use photochem_eqns, only: gravity + class(EvoAtmosphere), target, intent(inout) :: self + real(dp), intent(in) :: planet_mass, planet_radius + self%dat%planet_mass = planet_mass + self%dat%planet_radius = planet_radius + call gravity(self%dat%planet_radius, self%dat%planet_mass, & + self%var%nz, self%var%z, self%var%grav) + end subroutine + ! Below is mostly stuff needed in `evolve` routine module subroutine rebin_update_vertical_grid(self, usol_old, top_atmos, usol_new, err)