ability to change planet mass radius for evoatmosphere

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()

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.

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, 

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'.

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)
 

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

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

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)

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)