added gas_fluxes to EvoAtmosphere

photochem/cython/Atmosphere.pyx

@@ -184,8 +184,9 @@ cdef class Atmosphere:
     """
     cdef ndarray surf_fluxes = np.empty(self.dat.nq, np.double)
     cdef ndarray top_fluxes = np.empty(self.dat.nq, np.double)
+    cdef int nq = self.dat.nq
     cdef char err[ERR_LEN+1]
-    a_pxd.atmosphere_gas_fluxes_wrapper(&self._ptr, <double *>surf_fluxes.data, <double *>top_fluxes.data, err)
+    a_pxd.atmosphere_gas_fluxes_wrapper(&self._ptr, &nq, <double *>surf_fluxes.data, <double *>top_fluxes.data, err)
     if len(err.strip()) > 0:
       raise PhotoException(err.decode("utf-8").strip())
     surface = {}

photochem/cython/Atmosphere_pxd.pxd

@@ -19,7 +19,7 @@ cdef extern void atmosphere_photochemical_equilibrium_wrapper(void *ptr, bool *s
 
 cdef extern void atmosphere_out2atmosphere_txt_wrapper(void *ptr, char *filename, bool *overwrite, bool *clip, char *err)                         
 cdef extern void atmosphere_out2in_wrapper(void *ptr, char *err)
-cdef extern void atmosphere_gas_fluxes_wrapper(void *ptr, double *surf_fluxes, double *top_fluxes, char *err)
+cdef extern void atmosphere_gas_fluxes_wrapper(void *ptr, int *nq, double *surf_fluxes, double *top_fluxes, char *err)
 cdef extern void atmosphere_set_lower_bc_wrapper(void *ptr, char *species, char *bc_type, 
                                                     double *vdep, double *mix, double *flux, double *height, bool *missing, char *err)
 cdef extern void atmosphere_set_upper_bc_wrapper(void *ptr, char *species, 

photochem/cython/EvoAtmosphere.pyx

@@ -94,6 +94,34 @@ cdef class EvoAtmosphere:
     ea_pxd.evoatmosphere_out2atmosphere_txt_wrapper(&self._ptr, filename_c, &overwrite, &clip, err)  
     if len(err.strip()) > 0:
       raise PhotoException(err.decode("utf-8").strip())
+
+  def gas_fluxes(self):
+    """Computes gas fluxes at model boundaries in order to maintain
+    current atmospheric concentrations. Uses the densities stored in
+    self.wrk.usol.
+
+    Returns
+    -------
+    tuple
+        First element are the surface fluxes, and the second are top-of-atmosphere
+        fluxes.
+    """
+    cdef ndarray surf_fluxes = np.empty(self.dat.nq, np.double)
+    cdef ndarray top_fluxes = np.empty(self.dat.nq, np.double)
+    cdef int nq = self.dat.nq
+    cdef char err[ERR_LEN+1]
+    ea_pxd.evoatmosphere_gas_fluxes_wrapper(&self._ptr, &nq, <double *>surf_fluxes.data, <double *>top_fluxes.data, err)
+    if len(err.strip()) > 0:
+      raise PhotoException(err.decode("utf-8").strip())
+    surface = {}
+    names = self.dat.species_names
+    for i in range(self.dat.nq):
+      surface[names[i]] = surf_fluxes[i]
+    top = {}
+    names = self.dat.species_names
+    for i in range(self.dat.nq):
+      top[names[i]] = top_fluxes[i]
+    return surface, top
 
   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

@@ -16,7 +16,7 @@ cdef extern void evoatmosphere_create_wrapper(void *ptr, char *mechanism_file,
                                             void *wrk_ptr, char *err);
 
 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_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/Atmosphere_wrapper.f90

@@ -146,17 +146,18 @@ subroutine atmosphere_out2in_wrapper(ptr, err) bind(c)
 
   end subroutine
 
-  subroutine atmosphere_gas_fluxes_wrapper(ptr, surf_fluxes, top_fluxes, err) bind(c)
+  subroutine atmosphere_gas_fluxes_wrapper(ptr, nq, surf_fluxes, top_fluxes, err) bind(c)
     type(c_ptr), intent(in) :: ptr
-    real(c_double), intent(out) :: surf_fluxes(*)
-    real(c_double), intent(out) :: top_fluxes(*)
+    integer(c_int), intent(in) :: nq
+    real(c_double), intent(out) :: surf_fluxes(nq)
+    real(c_double), intent(out) :: top_fluxes(nq)
     character(kind=c_char), intent(out) :: err(err_len+1)
 
     character(:), allocatable :: err_f
     type(Atmosphere), pointer :: pc
 
     call c_f_pointer(ptr, pc)
-    call pc%gas_fluxes(surf_fluxes(1:pc%dat%nq),top_fluxes(1:pc%dat%nq),err_f)
+    call pc%gas_fluxes(surf_fluxes,top_fluxes,err_f)
     err(1) = c_null_char
     if (allocated(err_f)) then
       call copy_string_ftoc(err_f, err)

photochem/fortran/EvoAtmosphere_wrapper.f90

@@ -100,6 +100,24 @@ subroutine evoatmosphere_out2atmosphere_txt_wrapper(ptr, filename, overwrite, cl
 
   end subroutine
 
+  subroutine evoatmosphere_gas_fluxes_wrapper(ptr, nq, surf_fluxes, top_fluxes, err) bind(c)
+    type(c_ptr), intent(in) :: ptr
+    integer(c_int), intent(in) :: nq
+    real(c_double), intent(out) :: surf_fluxes(nq)
+    real(c_double), intent(out) :: top_fluxes(nq)
+    character(kind=c_char), intent(out) :: err(err_len+1)
+
+    character(:), allocatable :: err_f
+    type(EvoAtmosphere), pointer :: pc
+
+    call c_f_pointer(ptr, pc)
+    call pc%gas_fluxes(surf_fluxes,top_fluxes,err_f)
+    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/evoatmosphere/photochem_evoatmosphere.f90

@@ -56,6 +56,7 @@ function temp_dependent_albedo_fcn(T_surf) result(albedo)
 
     !!! photochem_evoatmosphere_utils.f90 !!!
     procedure :: out2atmosphere_txt
+    procedure :: gas_fluxes
     procedure :: rebin_update_vertical_grid
     procedure :: regrid_prep_atmosphere
 
@@ -221,6 +222,16 @@ module subroutine out2atmosphere_txt(self, filename, overwrite, clip, err)
       character(:), allocatable, intent(out) :: err
     end subroutine
 
+    !> Computes gas fluxes at model boundaries in order to maintain
+    !> current atmospheric concentrations. Uses the densities stored in
+    !> self%wrk%usol.
+    module subroutine gas_fluxes(self, surf_fluxes, top_fluxes, err)
+      class(EvoAtmosphere), target, intent(inout) :: self
+      real(dp), intent(out) :: surf_fluxes(:) !! surface fluxes (molecules/cm^2/s)
+      real(dp), intent(out) :: top_fluxes(:) !! top-of-atmosphere fluxes (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(:,:)

src/evoatmosphere/photochem_evoatmosphere_utils.f90

@@ -31,6 +31,56 @@ module subroutine out2atmosphere_txt(self, filename, overwrite, clip, err)
 
   end subroutine
 
+  module subroutine gas_fluxes(self, surf_fluxes, top_fluxes, err)
+    class(EvoAtmosphere), target, intent(inout) :: self
+    real(dp), intent(out) :: surf_fluxes(:)
+    real(dp), intent(out) :: top_fluxes(:)
+    character(:), allocatable, intent(out) :: err
+
+    real(dp) :: rhs(self%var%neqs)  
+    real(dp) :: diffusive_production
+    real(dp) :: chemical_production
+    type(PhotochemData), pointer :: dat
+    type(PhotochemVars), pointer :: var
+    type(PhotochemWrkEvo), pointer :: wrk
+
+    integer :: i
+
+    dat => self%dat
+    var => self%var
+    wrk => self%wrk
+
+    if (size(surf_fluxes) /= dat%nq .or. size(top_fluxes) /= dat%nq) then
+      err = "Input fluxes to gas_fluxes has the wrong dimensions"
+      return
+    endif
+
+    call self%right_hand_side_chem(wrk%usol, rhs, err)
+    if (allocated(err)) return
+
+    ! surface flux is molecules required to sustain the lower boundary
+    ! chemical production + diffusion production = total change in lower cell    
+    do i = 1,dat%nq
+      diffusive_production = (wrk%DU(i,1)*wrk%usol(i,2) + wrk%ADU(i,1)*wrk%usol(i,2) &
+                            + wrk%DD(i,1)*wrk%usol(i,1) + wrk%ADD(i,1)*wrk%usol(i,1)) &
+                              *var%dz(1)
+      chemical_production = rhs(i)*var%dz(1)
+      surf_fluxes(i) = -(diffusive_production + chemical_production)
+    enddo
+
+    ! fluxes going into or out of the top of the atmosphere.
+    do i = 1,dat%nq
+      diffusive_production = &
+         (wrk%DD(i,var%nz)*wrk%usol(i,var%nz) + wrk%ADD(i,var%nz)*wrk%usol(i,var%nz) &
+        + wrk%DL(i,var%nz)*wrk%usol(i,var%nz-1) + wrk%ADL(i,var%nz)*wrk%usol(i,var%nz-1)) &
+          *var%dz(var%nz)
+
+      chemical_production = rhs(i + (var%nz-1)*dat%nq)*var%dz(var%nz)
+      top_fluxes(i) = diffusive_production + chemical_production
+    enddo
+
+  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(:,:)