set_press_temp_eddy for EvoAtmosphere

photochem/cython/EvoAtmosphere.pyx

@@ -71,7 +71,7 @@ cdef class EvoAtmosphere:
     integration.
     """
     def __get__(self):
-      wrk = PhotochemWrk(alloc = False)
+      wrk = PhotochemWrkEvo(alloc = False)
       wrk._ptr = self._wrk_ptr
       return wrk
 
@@ -288,6 +288,48 @@ cdef class EvoAtmosphere:
                                        &trop_alt_, &trop_alt_present, err)
     if len(err.strip()) > 0:
       raise PhotoException(err.decode("utf-8").strip())
+
+  def set_press_temp_edd(self, ndarray[double, ndim=1] P, ndarray[double, ndim=1] T, ndarray[double, ndim=1] edd, trop_p = None, hydro_pressure = None):
+    """Given an input P, T, and edd, the code will find the temperature and eddy diffusion profile
+    on the current altitude-grid that matches the inputs.
+    
+    Parameters
+    ----------
+    P : ndarray[double,ndim=1]
+        Pressure (dynes/cm^2)
+    T : ndarray[double,ndim=1]
+        Temperature (K)
+    edd : ndarray[double,ndim=1]
+        Eddy diffusion (cm^2/s)
+    trop_p : float, optional
+        Tropopause pressure (dynes/cm^2). Only necessary if rainout == True, 
+        or fix_water_in_trop == True.
+    hydro_pressure : bool, optional
+        If True, then use hydrostatic pressure. If False then use the
+        actual pressure in the atmosphere. Default is True.
+    """
+    cdef char err[ERR_LEN+1]
+    cdef int P_dim1 = P.size
+    cdef int T_dim1 = T.size
+    cdef int edd_dim1 = edd.size
+
+    cdef double trop_p_ = 0.0
+    cdef bool trop_p_present = False
+    if trop_p != None:
+      trop_p_present = True
+      trop_p_ = trop_p
+
+    cdef bool hydro_pressure_ = True
+    cdef bool hydro_pressure_present = False
+    if hydro_pressure != None:
+      hydro_pressure_present = True
+      hydro_pressure_ = hydro_pressure
+
+    ea_pxd.evoatmosphere_set_press_temp_edd_wrapper(&self._ptr, &P_dim1, <double *>P.data, &T_dim1, <double *>T.data, &edd_dim1, <double *>edd.data,
+                                               &trop_p_, &trop_p_present, 
+                                               &hydro_pressure_, &hydro_pressure_present, err)
+    if len(err.strip()) > 0:
+      raise PhotoException(err.decode("utf-8").strip())
 
   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 

photochem/cython/EvoAtmosphere_pxd.pxd

@@ -25,6 +25,9 @@ cdef extern void evoatmosphere_set_upper_bc_wrapper(void *ptr, char *species,
 cdef extern void evoatmosphere_set_rate_fcn_wrapper(void *ptr, char *species_c, time_dependent_rate_fcn fcn, char *err)
 cdef extern void evoatmosphere_set_temperature_wrapper(void *ptr, int *nz, double *temperature, 
                                                     double *trop_alt, bool *trop_alt_present, char *err)
+cdef extern void evoatmosphere_set_press_temp_edd_wrapper(void *ptr, int *P_dim1, double *P, int *T_dim1, double *T, int *edd_dim1, double *edd,
+                                                      double *trop_p, bool *trop_p_present, 
+                                                      bool *hydro_pressure, bool *hydro_pressure_present, char *err)
 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/PhotochemWrk.pyx

@@ -1,5 +1,18 @@
 cimport PhotochemWrk_pxd as wrk_pxd
 
+cdef class PhotochemWrkEvo(PhotochemWrk):
+
+  property pressure_hydro:
+    """ndarray[double,dim=1], shape (nz). The hydrostatic pressure at the center of each 
+    atmospheric layer (dynes/cm^2).
+    """
+    def __get__(self):
+      cdef int dim1
+      wrk_pxd.photochemwrkevo_pressure_hydro_get_size(&self._ptr, &dim1)
+      cdef ndarray arr = np.empty(dim1, np.double)
+      wrk_pxd.photochemwrkevo_pressure_hydro_get(&self._ptr, &dim1, <double *>arr.data)
+      return arr
+
 cdef class PhotochemWrk:
   """This class contains data that changes during each step 
   when integrating the photochemical model

photochem/cython/PhotochemWrk_pxd.pxd

@@ -34,4 +34,8 @@ cdef extern void photochemwrk_optical_depth_get_size(void *ptr, int *dim1, int *
 cdef extern void photochemwrk_optical_depth_get(void *ptr, int *dim1, int *dim2, double *arr)
 
 cdef extern void photochemwrk_surf_radiance_get_size(void *ptr, int *dim1)
-cdef extern void photochemwrk_surf_radiance_get(void *ptr, int *dim1, double *arr)
+cdef extern void photochemwrk_surf_radiance_get(void *ptr, int *dim1, double *arr)
+
+# PhotochemWrkEvo
+cdef extern void photochemwrkevo_pressure_hydro_get_size(void *ptr, int *dim1)
+cdef extern void photochemwrkevo_pressure_hydro_get(void *ptr, int *dim1, double *arr)

photochem/fortran/EvoAtmosphere_wrapper.f90

@@ -245,6 +245,47 @@ subroutine evoatmosphere_set_temperature_wrapper(ptr, nz, temperature, &
 
   end subroutine
 
+  subroutine evoatmosphere_set_press_temp_edd_wrapper(ptr, P_dim1, P, T_dim1, T, edd_dim1, edd, &
+                                                      trop_p, trop_p_present, &
+                                                      hydro_pressure, hydro_pressure_present, err) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    integer(c_int), intent(in) :: P_dim1
+    real(c_double), intent(in) :: P(P_dim1)
+    integer(c_int), intent(in) :: T_dim1
+    real(c_double), intent(in) :: T(T_dim1)
+    integer(c_int), intent(in) :: edd_dim1
+    real(c_double), intent(in) :: edd(edd_dim1)
+    real(c_double), intent(in) :: trop_p
+    logical(c_bool), intent(in) :: trop_p_present
+    logical(c_bool), intent(in) :: hydro_pressure
+    logical(c_bool), intent(in) :: hydro_pressure_present
+    character(kind=c_char), intent(out) :: err(err_len+1)
+
+    character(:), allocatable :: err_f
+    logical :: hydro_pressure_f
+    type(EvoAtmosphere), pointer :: pc
+
+    call c_f_pointer(ptr, pc)
+
+    hydro_pressure_f = hydro_pressure
+
+    if (trop_p_present .and. hydro_pressure_present) then
+      call pc%set_press_temp_edd(P, T, edd, trop_p=trop_p, hydro_pressure=hydro_pressure_f, err=err_f)
+    elseif (trop_p_present .and. .not.hydro_pressure_present) then
+      call pc%set_press_temp_edd(P, T, edd, trop_p=trop_p, err=err_f)
+    elseif (.not.trop_p_present .and. hydro_pressure_present) then
+      call pc%set_press_temp_edd(P, T, edd, hydro_pressure=hydro_pressure_f, err=err_f)
+    else
+      call pc%set_press_temp_edd(P, T, edd, err=err_f)
+    endif
+    err(1) = c_null_char
+    if (allocated(err_f)) then
+      call copy_string_ftoc(err_f, err)
+    endif
+
+  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/PhotochemWrk_wrapper.f90

@@ -203,3 +203,22 @@ subroutine photochemwrk_surf_radiance_get(ptr, dim1, arr) bind(c)
     call c_f_pointer(ptr, wrk)
     arr = wrk%surf_radiance
   end subroutine
+
+! PhotochemWrkEvo
+
+  subroutine photochemwrkevo_pressure_hydro_get_size(ptr, dim1) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    integer(c_int), intent(out) :: dim1
+    type(PhotochemWrkEvo), pointer :: wrk
+    call c_f_pointer(ptr, wrk)
+    dim1 = size(wrk%pressure_hydro,1)
+  end subroutine
+
+  subroutine photochemwrkevo_pressure_hydro_get(ptr, dim1, arr) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    integer(c_int), intent(in) :: dim1
+    real(c_double), intent(out) :: arr(dim1)
+    type(PhotochemWrkEvo), pointer :: wrk
+    call c_f_pointer(ptr, wrk)
+    arr = wrk%pressure_hydro
+  end subroutine

photochem/fortran/photochem_c_api.f90

@@ -4,7 +4,7 @@ module photochem_c_api
   use photochem, only: EvoAtmosphere
   use photochem_types, only: PhotochemData
   use photochem_types, only: PhotochemVars
-  use photochem_types, only: PhotochemWrk
+  use photochem_types, only: PhotochemWrk, PhotochemWrkEvo
   use photochem_types, only: AtomConservation, ProductionLoss
   use photochem, only: err_len
   use photochem_const, only: s_str_len, m_str_len

src/evoatmosphere/photochem_evoatmosphere.f90

@@ -61,6 +61,7 @@ function temp_dependent_albedo_fcn(T_surf) result(albedo)
     procedure :: set_upper_bc
     procedure :: set_rate_fcn
     procedure :: set_temperature
+    procedure :: set_press_temp_edd
     procedure :: rebin_update_vertical_grid
     procedure :: regrid_prep_atmosphere
 
@@ -280,6 +281,20 @@ module subroutine set_temperature(self, temperature, trop_alt, err)
       character(:), allocatable, intent(out) :: err
     end subroutine
 
+    !> Given an input P, T, and edd, the code will find the temperature and eddy diffusion profile
+    !> on the current altitude-grid that matches the inputs.
+    module subroutine set_press_temp_edd(self, P, T, edd, trop_p, hydro_pressure, err)
+      class(EvoAtmosphere), target, intent(inout) :: self
+      real(dp), intent(in) :: P(:) !! Pressure (dynes/cm^2)
+      real(dp), intent(in) :: T(:) !! Temperature (K)
+      real(dp), intent(in) :: edd(:) !! Eddy diffusion (cm^2/s)
+      real(dp), optional, intent(in) :: trop_p !! Tropopause pressure (dynes/cm^2)
+      !> If .true., then use hydrostatic pressure. If .false. then use the
+      !> actual pressure in the atmosphere. Default is .true..
+      logical, optional, intent(in) :: hydro_pressure
+      character(:), allocatable, intent(out) :: err
+    end subroutine
+
     module subroutine rebin_update_vertical_grid(self, usol_old, top_atmos, usol_new, err)
       class(EvoAtmosphere), target, intent(inout) :: self
       real(dp), intent(in) :: usol_old(:,:)