diff --git a/CMakeLists.txt b/CMakeLists.txt
index f623759f6..74065413d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -265,6 +265,13 @@ add_executable(bin_average_comp src/bin_average_comp.F90)
 
 target_link_libraries(bin_average_comp partmclib)
 
+######################################################################
+# bin_deaverage_comp
+
+add_executable(bin_deaverage_comp src/bin_deaverage_comp.F90)
+
+target_link_libraries(bin_deaverage_comp partmclib)
+
 ######################################################################
 # bin_average_size
 
diff --git a/src/aero_state.F90 b/src/aero_state.F90
index 7c45482e6..165343369 100644
--- a/src/aero_state.F90
+++ b/src/aero_state.F90
@@ -1939,6 +1939,232 @@ subroutine aero_state_bin_average_comp(aero_state, bin_grid, aero_data)
 
   end subroutine aero_state_bin_average_comp
 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  !> Set each aerosol particle to have its original total volume, but
+  !> species volume ratios given by the group species volume ratio
+  !> within each bin. If no groups are specified, each particle contain a
+  !> single species (fully externally mixed). If groups are specified,
+  !> particles contain the volume fractions of the speices within the group.
+  !> This preserves per-particle total volumes.
+  subroutine aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, &
+       groups)
+
+    !> Aerosol state to de-average.
+    type(aero_state_t), intent(inout) :: aero_state
+    !> Bin grid to de-average within.
+    type(bin_grid_t), intent(in) :: bin_grid
+    !> Aerosol data.
+    type(aero_data_t), intent(in) :: aero_data
+    !> Sets of species names to group together.
+    character(len=*), optional :: groups(:,:)
+
+    real(kind=dp) :: species_volume_conc(aero_data_n_spec(aero_data))
+    real(kind=dp) :: total_volume_conc, particle_volume, num_conc
+    integer :: i_bin, i_class, i_entry, i_part, i_spec
+    real(kind=dp) :: factors(aero_data_n_spec(aero_data))
+    integer :: n_part_spec, start_val, i, n_part
+    integer :: species_group_numbers(aero_data_n_spec(aero_data))
+    integer :: n_group, i_group, i_name
+    real(kind=dp), allocatable :: group_fractions(:), &
+        group_volume_conc(:)
+    ! Flag to control algorithm selection - false uses low-variance method.
+    logical, parameter :: do_naive_algorithm = .false.
+    real(kind=dp), allocatable :: cumulative_vals(:)
+    real(kind=dp), allocatable :: factors_2d(:,:)
+    integer, allocatable :: particle_group(:)
+    real(kind=dp) :: x
+    integer, allocatable :: particle_index(:)
+    real(kind=dp) :: v_aerosol, v_bulk, p, v_group, v_particle, v_group_prev
+
+    call aero_state_sort(aero_state, aero_data, bin_grid)
+
+    if (present(groups)) then
+       n_group = size(groups, 1)
+       ! species_group_numbers(i_spec) will give the group number for
+       ! each species
+       species_group_numbers = 0
+       do i_group = 1,n_group
+          do i_name = 1,size(groups, 2)
+             if (len_trim(groups(i_group, i_name)) > 0) then
+                i_spec = aero_data_spec_by_name(aero_data, &
+                     groups(i_group, i_name))
+                call assert_msg(926066862, i_spec > 0, &
+                     "unknown species: " // trim(groups(i_group, i_name)))
+                species_group_numbers(i_spec) = i_group
+             end if
+          end do
+       end do
+       ! Assign left overs to their own groups
+       do i_spec = 1,aero_data_n_spec(aero_data)
+          if (species_group_numbers(i_spec) == 0) then
+             n_group = n_group + 1
+             species_group_numbers(i_spec) = n_group
+          end if
+       end do
+    else
+       ! Assign each species to a group
+       do i_spec = 1,aero_data_n_spec(aero_data)
+          species_group_numbers(i_spec) = i_spec
+       end do
+       n_group = aero_data_n_spec(aero_data)
+    end if
+
+    do i_bin = 1,bin_grid_size(bin_grid)
+       species_volume_conc = 0d0
+       total_volume_conc = 0d0
+       do i_class = 1,size(aero_state%awa%weight, 2)
+          do i_entry = 1,integer_varray_n_entry( &
+               aero_state%aero_sorted%size_class%inverse(i_bin, i_class))
+             i_part = aero_state%aero_sorted%size_class%inverse(i_bin, &
+                  i_class)%entry(i_entry)
+             num_conc = aero_weight_array_num_conc(aero_state%awa, &
+                  aero_state%apa%particle(i_part), aero_data)
+             particle_volume = aero_particle_volume( &
+                  aero_state%apa%particle(i_part))
+             species_volume_conc = species_volume_conc &
+                  + num_conc * aero_state%apa%particle(i_part)%vol
+             total_volume_conc = total_volume_conc + num_conc * particle_volume
+          end do
+       end do
+
+       allocate(group_volume_conc(n_group))
+       group_volume_conc = 0.0d0
+       do i_spec = 1,aero_data_n_spec(aero_data)
+          group_volume_conc(species_group_numbers(i_spec)) &
+             =  group_volume_conc(species_group_numbers(i_spec)) &
+             + species_volume_conc(i_spec)
+       end do
+
+       allocate(factors_2d(n_group,aero_data_n_spec(aero_data)))
+       factors_2d = 0.0d0
+       do i_group = 1,n_group
+          do i_spec = 1,aero_data_n_spec(aero_data)
+             if (species_group_numbers(i_spec) == i_group) then
+                factors_2d(i_group,i_spec) = &
+                     species_volume_conc(i_spec) &
+                     / group_volume_conc(i_group)
+             end if
+          end do
+       end do
+
+       if (do_naive_algorithm) then
+          do i_class = 1,size(aero_state%awa%weight, 2)
+             n_part = integer_varray_n_entry( &
+                  aero_state%aero_sorted%size_class%inverse(i_bin, i_class))
+             if (n_part > 0) then
+                allocate(group_fractions(n_group))
+                group_fractions = 0.0d0
+                do i_group = 1,n_group
+                   group_fractions(i_group) = n_part &
+                        * group_volume_conc(i_group) / total_volume_conc
+                end do
+                allocate(cumulative_vals(n_group+1))
+                cumulative_vals = 0.0d0
+                do i_group = 2,n_group + 1
+                   cumulative_vals(i_group) = cumulative_vals(i_group-1) &
+                        + group_fractions(i_group-1)
+                end do
+                allocate(particle_group(n_part))
+                do i_part = 1,n_part
+                   x = pmc_random() * n_part
+                   particle_group(i_part) = find_1d(n_group+1, &
+                       cumulative_vals, x)
+                end do
+
+                do i_entry = 1,n_part
+                   i_part = aero_state%aero_sorted%size_class%inverse( &
+                        i_bin, i_class)%entry(i_entry)
+                   particle_volume = aero_particle_volume( &
+                        aero_state%apa%particle(i_part))
+                   aero_state%apa%particle(i_part)%vol = 0.0d0
+                   i_group = particle_group(i_entry)
+                   do i_spec = 1,aero_data_n_spec(aero_data)
+                      aero_state%apa%particle(i_part)%vol(i_spec) &
+                           = particle_volume * factors_2d(i_group,i_spec)
+                   end do
+                end do
+                deallocate(particle_group)
+                deallocate(cumulative_vals)
+                deallocate(group_fractions)
+             end if
+          end do
+       else
+          n_part = 0
+          do i_class = 1,size(aero_state%awa%weight, 2)
+             n_part = n_part + integer_varray_n_entry( &
+                  aero_state%aero_sorted%size_class%inverse(i_bin, i_class))
+          end do
+
+          if (n_part > 0) then
+             allocate(particle_index(n_part))
+             i_part = 1
+             do i_class = 1,size(aero_state%awa%weight, 2)
+                do i_entry = 1, integer_varray_n_entry( &
+                     aero_state%aero_sorted%size_class%inverse(i_bin, i_class))
+                   particle_index(i_part) = &
+                        aero_state%aero_sorted%size_class%inverse( &
+                           i_bin, i_class)%entry(i_entry)
+                   i_part = i_part + 1
+                end do
+             end do
+
+             call pmc_rand_shuffle_array(particle_index, n_part)
+
+             i_entry = 1
+             v_aerosol = 0.0d0
+             v_bulk = 0.0d0
+             v_particle = 0.0d0
+             v_group = 0.0d0
+             v_group_prev = 0.0d0
+             i_group = 0
+             do while (i_entry <= n_part)
+                if (v_aerosol + v_particle >= v_bulk) then
+                   i_group = i_group + 1
+                   v_group = group_volume_conc(i_group)
+                   v_group_prev = v_bulk
+                   v_bulk = v_bulk + v_group
+                end if
+                i_part = particle_index(i_entry)
+                particle_volume = aero_particle_volume( &
+                     aero_state%apa%particle(i_part))
+                num_conc = aero_weight_array_num_conc(aero_state%awa, &
+                     aero_state%apa%particle(i_part), aero_data)
+                v_particle = particle_volume * num_conc
+                if (min(v_particle + v_aerosol, sum(group_volume_conc)) &
+                     <= v_bulk) then
+                   aero_state%apa%particle(i_part)%vol = 0.0d0
+                   do i_spec = 1,aero_data_n_spec(aero_data)
+                      aero_state%apa%particle(i_part)%vol(i_spec) = &
+                           particle_volume  * factors_2d(i_group,i_spec)
+                   end do
+                   v_aerosol = v_aerosol + v_particle
+                   i_entry = i_entry + 1
+                   v_particle = 0.0d0
+                else
+                   p = (v_bulk - max(v_aerosol, v_group_prev)) &
+                        / (v_particle)
+                   if (pmc_random() < p) then
+                      aero_state%apa%particle(i_part)%vol = 0.0d0
+                      do i_spec = 1,aero_data_n_spec(aero_data)
+                         aero_state%apa%particle(i_part)%vol(i_spec) = &
+                              particle_volume  * factors_2d(i_group,i_spec)
+                      end do
+                      v_aerosol = v_aerosol + v_particle
+                      i_entry = i_entry + 1
+                      v_particle = 0.0d0
+                   end if
+                end if
+             end do
+             deallocate(particle_index)
+          end if
+       end if
+       deallocate(group_volume_conc)
+       deallocate(factors_2d)
+    end do
+
+  end subroutine aero_state_bin_deaverage_comp
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
   !> Set each aerosol particle to have its original species ratios,
diff --git a/src/bin_average_comp.F90 b/src/bin_average_comp.F90
index 9d3d32902..0183a3fe7 100644
--- a/src/bin_average_comp.F90
+++ b/src/bin_average_comp.F90
@@ -66,6 +66,8 @@ program bin_average_comp
 
   call pmc_mpi_init()
 
+  call pmc_srand(0, pmc_mpi_rank())
+
   call bin_grid_make(bin_grid, BIN_GRID_TYPE_LOG, n_bin, diam2rad(d_min), &
        diam2rad(d_max))
 
diff --git a/src/bin_deaverage_comp.F90 b/src/bin_deaverage_comp.F90
new file mode 100644
index 000000000..b65612d38
--- /dev/null
+++ b/src/bin_deaverage_comp.F90
@@ -0,0 +1,99 @@
+! Copyright (C) 2009-2013 Matthew West
+! Licensed under the GNU General Public License version 2 or (at your
+! option) any later version. See the file COPYING for details.
+
+!> \file
+!> The bin_average_comp program.
+
+!> Read a NetCDF file, average the composition of all particles within
+!> each bin, and write the data out as another NetCDF file.
+program bin_deaverage_comp
+
+  use pmc_aero_state
+  use pmc_gas_data
+  use pmc_gas_state
+  use pmc_env_state
+  use pmc_aero_data
+  use pmc_bin_grid
+  use pmc_output
+  use pmc_rand
+  use netcdf
+
+  character(len=1000) :: in_filename, out_prefix
+  type(bin_grid_t) :: bin_grid
+  type(aero_data_t) :: aero_data
+  type(aero_state_t) :: aero_state
+  type(gas_data_t) :: gas_data
+  type(gas_state_t) :: gas_state
+  type(env_state_t) :: env_state
+  integer :: n_bin, index, i_repeat, output_type
+  real(kind=dp) :: d_min, d_max, time, del_t
+  character(len=1000) :: tmp_str
+  logical :: record_removals, dry_volume, record_optical
+  character(len=PMC_UUID_LEN) :: uuid
+  character(len=AERO_NAME_LEN), allocatable :: external_groups(:,:)
+
+  allocate(external_groups(3, 4)) ! 3 groups, max 4 species per group
+  external_groups(1,:) = ["OC    ", "BC    ", "      ", "      "]
+  external_groups(2,:) = ["API1  ", "API2  ", "LIM1  ", "LIM2  "]
+  external_groups(3,:) = ["SO4   ", "NO3   ", "NH4   ", "      "]
+
+  ! process commandline arguments
+  if (command_argument_count() .ne. 6) then
+     write(6,*) 'Usage: bin_average_comp <d_min> <d_max> <n_bin> ' &
+          // '<"wet" or "dry"> <input_filename> <output_prefix>'
+     write(6,*) ''
+     write(6,*) '  d_min: minimum bin diameter (m)'
+     write(6,*) '  d_max: maximum bin diameter (m)'
+     write(6,*) '  n_bin: number of bins'
+     write(6,*) '  wet/dry: average wet or dry sizes'
+     write(6,*) '  input_filename: like scenario_0001_00000001.nc'
+     write(6,*) '  output_prefix: like scenario_comp_average'
+     stop 2
+  endif
+  call get_command_argument(1, tmp_str)
+  d_min = string_to_real(tmp_str)
+  call get_command_argument(2, tmp_str)
+  d_max = string_to_real(tmp_str)
+  call get_command_argument(3, tmp_str)
+  n_bin = string_to_integer(tmp_str)
+  call get_command_argument(4, tmp_str)
+  if (trim(tmp_str) == "wet") then
+     dry_volume = .false.
+  elseif (trim(tmp_str) == "dry") then
+     dry_volume = .true.
+  else
+     write(6,*) 'Argument 4 must be "wet" or "dry", not ' &
+          // trim(tmp_str)
+     stop 1
+  end if
+  call get_command_argument(5, in_filename)
+  call get_command_argument(6, out_prefix)
+
+  call pmc_mpi_init()
+
+  call pmc_srand(0, pmc_mpi_rank())
+
+  call bin_grid_make(bin_grid, BIN_GRID_TYPE_LOG, n_bin, diam2rad(d_min), &
+       diam2rad(d_max))
+
+  call input_state(in_filename, index, time, del_t, i_repeat, uuid, &
+       aero_data, aero_state, gas_data, gas_state, env_state)
+
+  if (dry_volume) then
+     call aero_state_make_dry(aero_state, aero_data)
+  end if
+
+  call aero_state_bin_deaverage_comp(aero_state, bin_grid, aero_data, &
+       external_groups)
+
+  output_type = OUTPUT_TYPE_SINGLE
+  record_removals = .false.
+  record_optical = .true.
+  call output_state(out_prefix, output_type, aero_data, aero_state, &
+       gas_data, gas_state, env_state, index, time, del_t, i_repeat, &
+       record_removals, record_optical, uuid)
+
+  call pmc_mpi_finalize()
+
+end program bin_deaverage_comp
diff --git a/src/rand.F90 b/src/rand.F90
index e7640b75c..1e4472ef9 100644
--- a/src/rand.F90
+++ b/src/rand.F90
@@ -206,6 +206,53 @@ end function pmc_rand_int_gsl
 
   end function pmc_rand_int
 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  !> Shuffles randomly an array of integer values.
+  subroutine pmc_rand_shuffle_array(array, n_values)
+    !> Array of values.
+    integer, intent(inout) :: array(n_values)
+    !> Number of values.
+    integer, intent(in) :: n_values
+
+    integer :: temp
+    integer :: i, j
+    real(kind=dp) :: u
+
+#ifdef PMC_USE_GSL
+    integer(kind=c_int) :: n_c
+    integer(kind=c_int) :: array_c(n_values)
+#endif
+
+#ifdef PMC_USE_GSL
+#ifndef DOXYGEN_SKIP_DOC
+    interface
+       integer(kind=c_int) function pmc_rand_shuffle_gsl(array, n) bind(c)
+         use iso_c_binding
+         integer(kind=c_int), value :: n
+         integer(kind=c_int) :: array(n)
+       end function pmc_rand_shuffle_gsl
+    end interface
+#endif
+#endif
+
+#ifdef PMC_USE_GSL
+    n_c = int(n_values, kind=c_int)
+    array_c = int(array, kind=c_int)
+    call rand_check_gsl(388234845, pmc_rand_shuffle_gsl(array_c, n_c))
+    array = int(array_c)
+#else
+    do i=1,n_values-1
+         u = pmc_random()
+         j = i + floor((n_values-i+1)*u)
+         temp=array(j)
+         array(j)=array(i)
+         array(i)=temp
+    end do
+#endif
+
+  end subroutine pmc_rand_shuffle_array
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
   !> Round val to \c floor(val) or \c ceiling(val) with probability
diff --git a/src/rand_gsl.c b/src/rand_gsl.c
index db200abbf..c14bcd62a 100644
--- a/src/rand_gsl.c
+++ b/src/rand_gsl.c
@@ -150,3 +150,12 @@ int pmc_rand_binomial_gsl(int n, double p, int *harvest)
         *harvest = gsl_ran_binomial(pmc_rand_gsl_rng, p, u);
         return PMC_RAND_GSL_SUCCESS;
 }
+
+int pmc_rand_shuffle_gsl(int *array, int n)
+{
+        if (!pmc_rand_gsl_rng) {
+                return PMC_RAND_GSL_NOT_INIT;
+        }
+        gsl_ran_shuffle(pmc_rand_gsl_rng, array, n, sizeof (int));
+        return PMC_RAND_GSL_SUCCESS;
+}