boundary condition functions EvoAtmosphere

photochem/cython/EvoAtmosphere.pyx

@@ -122,6 +122,111 @@ cdef class EvoAtmosphere:
     for i in range(self.dat.nq):
       top[names[i]] = top_fluxes[i]
     return surface, top
+
+  def set_lower_bc(self, species = None, bc_type = None, vdep = None, den = None, press = None,
+                            flux = None, height = None):
+    """Sets a lower boundary condition.
+
+    Parameters
+    ----------
+    species : str
+        Species to set boundary condition
+    bc_type : str
+        Boundary condition type
+    vdep : float, optional
+        Deposition velocity (cm/s)
+    den : float, optional
+        Density (molecules/cm^3)
+    press : float, optional
+        Pressure (dynes/cm^2)
+    flux : float, optional
+        Flux (molecules/cm^2/s)
+    height : float, optional
+        Height in atmosphere (km)
+    """
+    cdef bytes species_b = pystring2cstring(species)
+    cdef char *species_c = species_b
+    cdef bytes bc_type_b = pystring2cstring(bc_type)
+    cdef char *bc_type_c = bc_type_b
+    cdef double vdep_c = 0
+    cdef double den_c = 0
+    cdef double press_c = 0
+    cdef double flux_c = 0
+    cdef double height_c = 0
+    cdef bool missing = False
+    if bc_type == 'vdep':
+      if vdep == None:
+        missing = True
+      else:
+        vdep_c = vdep
+    elif bc_type == 'den':
+      if den == None:
+        missing = True
+      else:
+        den_c = den
+    elif bc_type == 'press':
+      if press == None:
+        missing = True
+      else:
+        press_c = press
+    elif bc_type == 'flux':
+      if flux == None:
+        missing = True
+      else:
+        flux_c = flux
+    elif bc_type == 'vdep + dist flux':
+      if vdep == None or flux == None or height == None:
+        missing = True
+      else:
+        vdep_c = vdep
+        flux_c = flux
+        height_c = height
+    elif bc_type == 'Moses':
+      pass
+
+    cdef char err[ERR_LEN+1]
+    ea_pxd.evoatmosphere_set_lower_bc_wrapper(&self._ptr, species_c, bc_type_c, 
+                                      &vdep_c, &den_c, &press_c, &flux_c, &height_c, &missing, err);
+    if len(err.strip()) > 0:
+      raise PhotoException(err.decode("utf-8").strip())
+
+  def set_upper_bc(self, species = None, bc_type = None, veff = None,flux = None):
+    """Sets upper boundary condition.
+
+    Parameters
+    ----------
+    species : str
+        Species to set boundary condition
+    bc_type : str
+        Boundary condition type
+    veff : float, optional
+        effusion velocity (cm/s)
+    flux : float, optional
+        Flux (molecules/cm^2/s)
+    """
+    cdef bytes species_b = pystring2cstring(species)
+    cdef char *species_c = species_b
+    cdef bytes bc_type_b = pystring2cstring(bc_type)
+    cdef char *bc_type_c = bc_type_b
+    cdef double veff_c = 0
+    cdef double flux_c = 0
+    cdef bool missing = False
+    if bc_type == 'veff':
+      if veff == None:
+        missing = True
+      else:
+        veff_c = veff
+    elif bc_type == 'flux':
+      if flux == None:
+        missing = True
+      else:
+        flux_c = flux
+
+    cdef char err[ERR_LEN+1]
+    ea_pxd.evoatmosphere_set_upper_bc_wrapper(&self._ptr, species_c, bc_type_c, 
+                                      &veff_c, &flux_c, &missing, 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

@@ -17,6 +17,10 @@ cdef extern void evoatmosphere_create_wrapper(void *ptr, char *mechanism_file,
 
 cdef extern void evoatmosphere_out2atmosphere_txt_wrapper(void *ptr, char *filename, bool *overwrite, bool *clip, char *err)   
 cdef extern void evoatmosphere_gas_fluxes_wrapper(void *ptr, int *nq, double *surf_fluxes, double *top_fluxes, char *err)
+cdef extern void evoatmosphere_set_lower_bc_wrapper(void *ptr, char *species, char *bc_type, 
+                                                    double *vdep, double *den, double *press, double *flux, double *height, bool *missing, char *err)
+cdef extern void evoatmosphere_set_upper_bc_wrapper(void *ptr, char *species, 
+                                                    char *bc_type, double *veff, double *flux, bool *missing, 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/fortran/EvoAtmosphere_wrapper.f90

@@ -118,6 +118,75 @@ subroutine evoatmosphere_gas_fluxes_wrapper(ptr, nq, surf_fluxes, top_fluxes, er
     endif
   end subroutine
 
+  subroutine evoatmosphere_set_lower_bc_wrapper(ptr, species, bc_type, vdep, den, press, flux, height, missing, err) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    character(kind=c_char), intent(in) :: species(*)
+    character(kind=c_char), intent(in) :: bc_type(*)
+    real(c_double), intent(in) :: vdep
+    real(c_double), intent(in) :: den
+    real(c_double), intent(in) :: press
+    real(c_double), intent(in) :: flux
+    real(c_double), intent(in) :: height
+    logical(c_bool), intent(in) :: missing
+    character(kind=c_char), intent(out) :: err(err_len+1)
+
+    character(len=:), allocatable :: species_f
+    character(len=:), allocatable :: bc_type_f
+
+    character(:), allocatable :: err_f
+    type(EvoAtmosphere), pointer :: pc
+    call c_f_pointer(ptr, pc)
+
+    allocate(character(len=len_cstring(species))::species_f)
+    allocate(character(len=len_cstring(bc_type))::bc_type_f)
+
+    call copy_string_ctof(species, species_f)
+    call copy_string_ctof(bc_type, bc_type_f)
+
+    if (missing) then
+      call pc%set_lower_bc(species_f, bc_type_f, err=err_f)
+    else
+      call pc%set_lower_bc(species_f, bc_type_f, vdep=vdep, den=den, press=press, flux=flux, height=height, 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_set_upper_bc_wrapper(ptr, species, bc_type, veff, flux, missing, err) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    character(kind=c_char), intent(in) :: species(*)
+    character(kind=c_char), intent(in) :: bc_type(*)
+    real(c_double), intent(in) :: veff
+    real(c_double), intent(in) :: flux
+    logical(c_bool), intent(in) :: missing
+    character(kind=c_char), intent(out) :: err(err_len+1)
+
+    character(len=:), allocatable :: species_f
+    character(len=:), allocatable :: bc_type_f
+
+    character(:), allocatable :: err_f
+    type(EvoAtmosphere), pointer :: pc
+    call c_f_pointer(ptr, pc)
+
+    allocate(character(len=len_cstring(species))::species_f)
+    allocate(character(len=len_cstring(bc_type))::bc_type_f)
+
+    call copy_string_ctof(species, species_f)
+    call copy_string_ctof(bc_type, bc_type_f)
+
+    if (missing) then
+      call pc%set_upper_bc(species_f, bc_type_f, err=err_f)
+    else
+      call pc%set_upper_bc(species_f, bc_type_f, veff=veff, flux=flux, 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

src/atmosphere/photochem_atmosphere_utils.f90

@@ -419,6 +419,12 @@ module subroutine set_lower_bc(self, species, bc_type, vdep, mix, flux, height,
             "' because it is not in the list of species"
       return
     endif
+    if (self%dat%there_are_particles) then
+      if (ind(1) <= self%dat%np .and. bc_type /= 'vdep') then
+        err = 'Particles must have a deposition velocity lower boundary condition.'
+        return
+      endif
+    endif
 
     if (self%dat%fix_water_in_trop) then
       if (species == "H2O") then

src/evoatmosphere/photochem_evoatmosphere.f90

@@ -57,6 +57,8 @@ function temp_dependent_albedo_fcn(T_surf) result(albedo)
     !!! photochem_evoatmosphere_utils.f90 !!!
     procedure :: out2atmosphere_txt
     procedure :: gas_fluxes
+    procedure :: set_lower_bc
+    procedure :: set_upper_bc
     procedure :: rebin_update_vertical_grid
     procedure :: regrid_prep_atmosphere
 
@@ -232,6 +234,29 @@ module subroutine gas_fluxes(self, surf_fluxes, top_fluxes, err)
       character(:), allocatable, intent(out) :: err
     end subroutine
 
+    !> Sets a lower boundary condition.
+    module subroutine set_lower_bc(self, species, bc_type, vdep, den, press, flux, height, err)
+      class(EvoAtmosphere), intent(inout) :: self
+      character(len=*), intent(in) :: species !! Species to set boundary condition
+      character(len=*), intent(in) :: bc_type !! Boundary condition type
+      real(dp), optional, intent(in) :: vdep !! Deposition velocity (cm/s)
+      real(dp), optional, intent(in) :: den !! density (molecules/cm^3)
+      real(dp), optional, intent(in) :: press !! pressure (dynes/cm^2)
+      real(dp), optional, intent(in) :: flux !! Flux (molecules/cm^2/s)
+      real(dp), optional, intent(in) :: height !! Height in atmosphere (km)
+      character(:), allocatable, intent(out) :: err
+    end subroutine
+
+    !> Sets upper boundary condition
+    module subroutine set_upper_bc(self, species, bc_type, veff, flux, err)
+      class(EvoAtmosphere), intent(inout) :: self
+      character(len=*), intent(in) :: species !! Species to set boundary condition
+      character(len=*), intent(in) :: bc_type !! Boundary condition type
+      real(dp), optional, intent(in) :: veff !! effusion velocity (cm/s)
+      real(dp), optional, intent(in) :: flux !! Flux (molecules/cm^2/s)
+      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(:,:)