diff --git a/photochem/cython/PhotochemWrk.pyx b/photochem/cython/PhotochemWrk.pyx
index 166c8fa2..763133fd 100644
--- a/photochem/cython/PhotochemWrk.pyx
+++ b/photochem/cython/PhotochemWrk.pyx
@@ -147,6 +147,39 @@ cdef class PhotochemWrk:
       wrk_pxd.photochemwrk_rx_rates_get(self._ptr, &dim1, &dim2, <double *>arr.data)
       return arr
       
+  property transport_rates:
+    """ndarray[double,dim=2], shape (nz,nrT). Reaction rate constants in various units
+    involving molecules cm^3 and s. These rates include 3rd body contributions.
+    """
+    def __get__(self):
+      cdef int dim1, dim2
+      wrk_pxd.photochemwrk_transport_rates_get_size(self._ptr, &dim1, &dim2)
+      cdef ndarray arr = np.empty((dim1, dim2), np.double, order="F")
+      wrk_pxd.photochemwrk_transport_rates_get(self._ptr, &dim1, &dim2, <double *>arr.data)
+      return arr
+
+  property rainout_rates:
+    """ndarray[double,dim=2], shape (nz,nrT). Reaction rate constants in various units
+    involving molecules cm^3 and s. These rates include 3rd body contributions.
+    """
+    def __get__(self):
+      cdef int dim1, dim2
+      wrk_pxd.photochemwrk_rainout_rates_get_size(self._ptr, &dim1, &dim2)
+      cdef ndarray arr = np.empty((dim1, dim2), np.double, order="F")
+      wrk_pxd.photochemwrk_rainout_rates_get(self._ptr, &dim1, &dim2, <double *>arr.data)
+      return arr
+
+  property distributed_fluxes:
+    """ndarray[double,dim=2], shape (nz,nrT). Reaction rate constants in various units
+    involving molecules cm^3 and s. These rates include 3rd body contributions.
+    """
+    def __get__(self):
+      cdef int dim1, dim2
+      wrk_pxd.photochemwrk_distributed_fluxes_get_size(self._ptr, &dim1, &dim2)
+      cdef ndarray arr = np.empty((dim1, dim2), np.double, order="F")
+      wrk_pxd.photochemwrk_distributed_fluxes_get(self._ptr, &dim1, &dim2, <double *>arr.data)
+      return arr
+
   property mubar:
     """ndarray[double,dim=1], shape (nz). The mean molar mass of each atmospheric layer
     (g/mol)
@@ -201,3 +234,4 @@ cdef class PhotochemWrk:
       wrk_pxd.photochemwrk_surf_radiance_get(self._ptr, &dim1, <double *>arr.data)
       return arr
   
+
diff --git a/photochem/cython/PhotochemWrk_pxd.pxd b/photochem/cython/PhotochemWrk_pxd.pxd
index edf3e704..3aadb366 100644
--- a/photochem/cython/PhotochemWrk_pxd.pxd
+++ b/photochem/cython/PhotochemWrk_pxd.pxd
@@ -41,6 +41,15 @@ cdef extern void photochemwrk_densities_get(PhotochemWrk *ptr, int *dim1, int *d
 cdef extern void photochemwrk_rx_rates_get_size(PhotochemWrk *ptr, int *dim1, int *dim2)
 cdef extern void photochemwrk_rx_rates_get(PhotochemWrk *ptr, int *dim1, int *dim2, double *arr)
 
+cdef extern void photochemwrk_transport_rates_get_size(PhotochemWrk *ptr, int *dim1, int *dim2)
+cdef extern void photochemwrk_transport_rates_get(PhotochemWrk *ptr, int *dim1, int *dim2, double *arr)
+
+cdef extern void photochemwrk_rainout_rates_get_size(PhotochemWrk *ptr, int *dim1, int *dim2)
+cdef extern void photochemwrk_rainout_rates_get(PhotochemWrk *ptr, int *dim1, int *dim2, double *arr)
+
+cdef extern void photochemwrk_distributed_fluxes_get_size(PhotochemWrk *ptr, int *dim1, int *dim2)
+cdef extern void photochemwrk_distributed_fluxes_get(PhotochemWrk *ptr, int *dim1, int *dim2, double *arr)
+
 cdef extern void photochemwrk_mubar_get_size(PhotochemWrk *ptr, int *dim1)
 cdef extern void photochemwrk_mubar_get(PhotochemWrk *ptr, int *dim1, double *arr)
 
@@ -58,4 +67,4 @@ cdef extern void photochemwrk_surf_radiance_get(PhotochemWrk *ptr, int *dim1, do
 
 # PhotochemWrkEvo
 cdef extern void photochemwrkevo_pressure_hydro_get_size(PhotochemWrkEvo *ptr, int *dim1)
-cdef extern void photochemwrkevo_pressure_hydro_get(PhotochemWrkEvo *ptr, int *dim1, double *arr)
\ No newline at end of file
+cdef extern void photochemwrkevo_pressure_hydro_get(PhotochemWrkEvo *ptr, int *dim1, double *arr)
diff --git a/photochem/fortran/PhotochemWrk_wrapper.f90 b/photochem/fortran/PhotochemWrk_wrapper.f90
index e8bdf2dd..94731a09 100644
--- a/photochem/fortran/PhotochemWrk_wrapper.f90
+++ b/photochem/fortran/PhotochemWrk_wrapper.f90
@@ -184,6 +184,60 @@ subroutine photochemwrk_rx_rates_get(ptr, dim1, dim2, arr) bind(c)
     arr = wrk%rx_rates
   end subroutine
 
+  subroutine photochemwrk_transport_rates_get_size(ptr, dim1, dim2) bind(c)
+    type(c_ptr), value, intent(in) :: ptr
+    integer(c_int), intent(out) :: dim1, dim2
+    type(PhotochemWrk), pointer :: wrk
+    call c_f_pointer(ptr, wrk)
+    dim1 = size(wrk%transport_rates,1)
+    dim2 = size(wrk%transport_rates,2)
+  end subroutine
+  
+  subroutine photochemwrk_transport_rates_get(ptr, dim1, dim2, arr) bind(c)
+    type(c_ptr), value, intent(in) :: ptr
+    integer(c_int), intent(in) :: dim1, dim2
+    real(c_double), intent(out) :: arr(dim1, dim2)
+    type(PhotochemWrk), pointer :: wrk
+    call c_f_pointer(ptr, wrk)
+    arr = wrk%transport_rates
+  end subroutine
+
+  subroutine photochemwrk_rainout_rates_get_size(ptr, dim1, dim2) bind(c)
+    type(c_ptr), value, intent(in) :: ptr
+    integer(c_int), intent(out) :: dim1, dim2
+    type(PhotochemWrk), pointer :: wrk
+    call c_f_pointer(ptr, wrk)
+    dim1 = size(wrk%rainout_rates,1)
+    dim2 = size(wrk%rainout_rates,2)
+  end subroutine
+  
+  subroutine photochemwrk_rainout_rates_get(ptr, dim1, dim2, arr) bind(c)
+    type(c_ptr), value, intent(in) :: ptr
+    integer(c_int), intent(in) :: dim1, dim2
+    real(c_double), intent(out) :: arr(dim1, dim2)
+    type(PhotochemWrk), pointer :: wrk
+    call c_f_pointer(ptr, wrk)
+    arr = wrk%rainout_rates
+  end subroutine
+
+  subroutine photochemwrk_distributed_fluxes_get_size(ptr, dim1, dim2) bind(c)
+    type(c_ptr), value, intent(in) :: ptr
+    integer(c_int), intent(out) :: dim1, dim2
+    type(PhotochemWrk), pointer :: wrk
+    call c_f_pointer(ptr, wrk)
+    dim1 = size(wrk%distributed_fluxes,1)
+    dim2 = size(wrk%distributed_fluxes,2)
+  end subroutine
+  
+  subroutine photochemwrk_distributed_fluxes_get(ptr, dim1, dim2, arr) bind(c)
+    type(c_ptr), value, intent(in) :: ptr
+    integer(c_int), intent(in) :: dim1, dim2
+    real(c_double), intent(out) :: arr(dim1, dim2)
+    type(PhotochemWrk), pointer :: wrk
+    call c_f_pointer(ptr, wrk)
+    arr = wrk%distributed_fluxes
+  end subroutine
+
   subroutine photochemwrk_mubar_get_size(ptr, dim1) bind(c)
     type(c_ptr), value, intent(in) :: ptr
     integer(c_int), intent(out) :: dim1
@@ -289,4 +343,4 @@ subroutine photochemwrkevo_pressure_hydro_get(ptr, dim1, arr) bind(c)
     type(PhotochemWrkEvo), pointer :: wrk
     call c_f_pointer(ptr, wrk)
     arr = wrk%pressure_hydro
-  end subroutine
\ No newline at end of file
+  end subroutine
diff --git a/src/evoatmosphere/photochem_evoatmosphere_rhs.f90 b/src/evoatmosphere/photochem_evoatmosphere_rhs.f90
index c03e77da..b68aba6f 100644
--- a/src/evoatmosphere/photochem_evoatmosphere_rhs.f90
+++ b/src/evoatmosphere/photochem_evoatmosphere_rhs.f90
@@ -795,33 +795,63 @@ module subroutine rhs_evo_gas(self, neqs, tn, usol_flat, rhs, err)
     do j = 2,var%nz-1
       do i = 1,dat%nq
         k = i + (j-1)*dat%nq
-        rhs(k) = rhs(k) + wrk%DU(i,j)*wrk%usol(i,j+1) + wrk%ADU(i,j)*wrk%usol(i,j+1) &
-                        + wrk%DD(i,j)*wrk%usol(i,j) + wrk%ADD(i,j)*wrk%usol(i,j) &
-                        + wrk%DL(i,j)*wrk%usol(i,j-1) + wrk%ADL(i,j)*wrk%usol(i,j-1)
+  
+        wrk%transport_rates(i, j) = wrk%DU(i,j)*wrk%usol(i,j+1) + wrk%ADU(i,j)*wrk%usol(i,j+1) &
+        + wrk%DD(i,j)*wrk%usol(i,j) + wrk%ADD(i,j)*wrk%usol(i,j) &
+        + wrk%DL(i,j)*wrk%usol(i,j-1) + wrk%ADL(i,j)*wrk%usol(i,j-1)
+
+        rhs(k) = rhs(k) + wrk%transport_rates(i, j)
+
       enddo
     enddo
+
+    
     
     ! Lower boundary
     do i = 1,dat%nq
       if (var%lowerboundcond(i) == VelocityBC .or. &
           var%lowerboundcond(i) == VelocityDistributedFluxBC) then
-        rhs(i) = rhs(i) + 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) &
-                        - wrk%lower_vdep_copy(i)*wrk%usol(i,1)/var%dz(1)
+
+        wrk%transport_rates(i, 1) = &
+            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) &
+          - wrk%lower_vdep_copy(i)*wrk%usol(i,1)/var%dz(1)
+
+        rhs(i) = rhs(i) + wrk%transport_rates(i, 1) 
       elseif (var%lowerboundcond(i) == DensityBC .or. &
               var%lowerboundcond(i) == PressureBC) then
+
+        ! surface flux is molecules required to sustain the lower boundary
+    ! chemical production + diffusion production = total change in lower cell    
+
+      wrk%transport_rates(i, 1)  = &
+          !diffusion 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) &
+        ! chemical production 
+        + rhs(i)
+                            
         rhs(i) = 0.0_dp
       elseif (var%lowerboundcond(i) == FluxBC) then
-        rhs(i) = rhs(i) + 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%lower_flux(i)/var%dz(1)
+
+        wrk%transport_rates(i, 1) = &
+            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%lower_flux(i)/var%dz(1)
+
+        rhs(i) = rhs(i) + wrk%transport_rates(i, 1)
+
       ! Moses (2001) boundary condition for gas giants
       ! A deposition velocity controled by how quickly gases
       ! turbulantly mix vertically
       elseif (var%lowerboundcond(i) == MosesBC) then
-        rhs(i) = rhs(i) + 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%edd(1)/wrk%scale_height(1))*wrk%usol(i,1)/var%dz(1)
+        
+        wrk%transport_rates(i, 1) = &
+            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%edd(1)/wrk%scale_height(1))*wrk%usol(i,1)/var%dz(1)
+
+        rhs(i) = rhs(i) + wrk%transport_rates(i, 1)
       endif
     enddo
 
@@ -829,13 +859,21 @@ module subroutine rhs_evo_gas(self, neqs, tn, usol_flat, rhs, err)
     do i = 1,dat%nq
       k = i + (var%nz-1)*dat%nq
       if (var%upperboundcond(i) == VelocityBC) then
-        rhs(k) = rhs(k) + 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) &
-                        - wrk%upper_veff_copy(i)*wrk%usol(i,var%nz)/var%dz(var%nz)    
+
+        wrk%transport_rates(i, var%nz) = &
+            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) &
+          - wrk%upper_veff_copy(i)*wrk%usol(i,var%nz)/var%dz(var%nz)  
+
+        rhs(k) = rhs(k) + wrk%transport_rates(i, var%nz) 
+
       elseif (var%upperboundcond(i) == FluxBC) then
-        rhs(k) = rhs(k) + 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%upper_flux(i)/var%dz(var%nz)
+
+        wrk%transport_rates(i, var%nz)  = & 
+            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%upper_flux(i)/var%dz(var%nz)
+        rhs(k) = rhs(k) + wrk%transport_rates(i, var%nz)
       endif
     enddo
 
@@ -856,7 +894,8 @@ module subroutine rhs_evo_gas(self, neqs, tn, usol_flat, rhs, err)
         ztop1 = var%z(jdisth) + 0.5e0_dp*var%dz(jdisth)
         do j = 2,jdisth
           k = i + (j-1)*dat%nq
-          rhs(k) = rhs(k) + 2.0_dp*var%lower_flux(i)*(ztop1-var%z(j))/(ztop**2.0_dp)
+          wrk%distributed_fluxes(i, j) = 2.0_dp*var%lower_flux(i)*(ztop1-var%z(j))/(ztop**2.0_dp)
+          rhs(k) = rhs(k) + wrk%distributed_fluxes(i, j)
         enddo
         endif
       endif
@@ -1236,3 +1275,4 @@ module subroutine production_and_loss(self, species, usol, pl, err)
 end submodule
 
 
+
diff --git a/src/photochem_types.f90 b/src/photochem_types.f90
index f913a663..baaac328 100644
--- a/src/photochem_types.f90
+++ b/src/photochem_types.f90
@@ -564,6 +564,8 @@ subroutine time_dependent_rate_fcn(tn, nz, rate)
     real(dp), allocatable :: densities(:,:) !! (nsp+1,nz)
     real(dp), allocatable :: density(:) !! (nz)
     real(dp), allocatable :: rx_rates(:,:) !! (nz,nrT)
+    real(dp), allocatable :: transport_rates(:,:) !! (nq,nz)
+    real(dp), allocatable :: distributed_fluxes(:,:) !! (nq,nz)
     real(dp), allocatable :: mubar(:) !! (nz)
     real(dp), allocatable :: pressure(:) !! (nz)
     real(dp), allocatable :: H2O_rh(:) !! (nz)
@@ -651,6 +653,8 @@ subroutine init_PhotochemWrk(self, nsp, np, nq, nz, nrT, kj, nw)
       deallocate(self%H2O_sat_mix)
       deallocate(self%densities)
       deallocate(self%rx_rates)
+      deallocate(self%transport_rates)
+      deallocate(self%distributed_fluxes)
       deallocate(self%prates)
       deallocate(self%surf_radiance)
       deallocate(self%amean_grd)
@@ -685,6 +689,8 @@ subroutine init_PhotochemWrk(self, nsp, np, nq, nz, nrT, kj, nw)
     allocate(self%H2O_sat_mix(nz))
     allocate(self%densities(nsp+1,nz))
     allocate(self%rx_rates(nz,nrT))
+    allocate(self%transport_rates(nq, nz))
+    allocate(self%distributed_fluxes(nq, nz))
     allocate(self%prates(nz,kj))
     allocate(self%surf_radiance(nw))
     allocate(self%amean_grd(nz,nw))
@@ -772,3 +778,4 @@ subroutine SundialsDataFinalizer_final(self)
 end module
 
 
+