diff --git a/README.md b/README.md
index bb223817..694710be 100644
--- a/README.md
+++ b/README.md
@@ -35,11 +35,14 @@ Utility functions
 
 * Array functions
 * Assertions
-* Emulated intrinsic functions
+* Emulated intrinsic functions: `findloc`
+* Emulated collective subroutines: `co_sum`, `co_broadcast`
+* User-defined collective subroutines: `co_all`
 * String functions
 
 Classes
 -------
+* Parallel data partitioning and gathering
 * (Co-)Object pattern abstract parent
 
 Prerequisites
@@ -65,13 +68,19 @@ fpm test \
   --flag "-Wall" \
   --flag "-std=f2018" \
   --flag "-DCOMPILER_LACKS_COLLECTIVE_SUBROUTINES" \
-  --flag "-DCOMPILER_LACKS_FINDLOC" 
+  --flag "-DCOMPILER_LACKS_FINDLOC"
 ```
-where the `COMPILER_LACKS_*` flags exercise the Sourcery Library's 
-emulated instrinsic procedures, which are intended for use with 
+where the `COMPILER_LACKS_*` flags exercise the Sourcery Library's
+emulated instrinsic procedures, which are intended for use with
 compiler versions that lack support for the named features.  Delete
 those flags with compilers that support these features.
 
+Build documentation
+-------------------
+```zsh
+ford doc/ford-documentation.md
+```
+
 [GNU Fortran]: https://gcc.gnu.org
 [OpenCoarrays]: https://github.com/sourceryinstitute/opencoarrays
 [fpm]: https://github.com/fortran-lang/fpm
diff --git a/src/data-partition-implementation.F90 b/src/data-partition-implementation.F90
new file mode 100644
index 00000000..fbdd62b5
--- /dev/null
+++ b/src/data-partition-implementation.F90
@@ -0,0 +1,148 @@
+submodule(data_partition_interface) data_partition_implementation
+  use assertions_interface, only : assert, assertions
+  implicit none
+
+contains
+
+  module procedure define_partitions
+
+    if (allocated(first_datum)) deallocate(first_datum)
+    if (allocated(last_datum)) deallocate(last_datum)
+
+    associate( ni => num_images() )
+
+      call assert( ni<=cardinality, "sufficient data for distribution across images")
+
+      allocate(first_datum(ni), last_datum(ni))
+
+      block
+        integer i, image
+        do image=1,ni
+          associate( remainder => mod(cardinality, ni), quotient => cardinality/ni )
+            first_datum(image) = sum([(quotient+overflow(i, remainder), i=1, image-1)]) + 1
+            last_datum(image) = first_datum(image) + quotient + overflow(image, remainder) - 1
+          end associate
+        end do
+      end block
+    end associate
+
+#ifdef FORD
+  end procedure
+#else
+  contains
+#endif
+
+    pure function overflow(im, excess) result(extra_datum)
+      integer, intent(in) :: im, excess
+      integer extra_datum
+      extra_datum= merge(1,0,im<=excess)
+    end function
+
+#ifndef FORD
+  end procedure
+#endif
+
+  module procedure first
+    if (assertions) call assert( allocated(first_datum), "allocated(first_datum)")
+    first_index= first_datum( image_number )
+  end procedure
+
+  module procedure last
+    if (assertions) call assert( allocated(last_datum), "allocated(last_datum)")
+    last_index = last_datum( image_number )
+  end procedure
+
+  module procedure gather_real_1D_array
+
+    if (present(dim)) call assert (dim==1, "dimensioned partitioned == 1")
+
+    associate( me => this_image() )
+      write(6,*) 'gather_real_1D_array(): executing on image', me
+      flush(6)
+      associate( first=>first(me), last=>last(me) )
+        if (.not. present(result_image)) then
+          a(1:first-1)  = 0.
+          a(last+1:)  = 0.
+          call co_sum(a)
+        else
+          block
+            real(real64), allocatable, dimension(:) :: a_lower, a_upper
+            a_lower = a(1:first-1)
+            a_upper = a(last+1:)
+            a(1:first-1)  = 0.
+            a(last+1:)  = 0.
+            call co_sum(a, result_image=result_image)
+            if (result_image /= me) then
+              a(1:first-1) = a_lower
+              a(last+1:) = a_upper
+            end if
+          end block
+        end if
+      end associate
+    end associate
+  end procedure
+
+  module procedure gather_real_2D_array
+
+    integer dim_
+    if (present(dim)) then
+      dim_ = dim
+    else
+      dim_ = 2
+    end if
+
+    associate( me => this_image() )
+      write(6,*) 'gather_real_2D_array(): executing on image', me
+      flush(6)
+      associate( first => first(me), last => last(me) )
+        if (.not. present(result_image)) then
+          select case(dim_)
+            case(1)
+              a(1:first-1, :) = 0.
+              a(last+1:, :) = 0.
+            case(2)
+              a(:, 1:first-1) = 0.
+              a(:, last+1:) = 0.
+            case default
+              error stop "gather_real_2D_array: invalid dim argument"
+          end select
+          call co_sum(a)
+        else
+          block
+            real(real64), allocatable, dimension(:,:) :: a_lower, a_upper
+            select case(dim_)
+              case(1)
+                a_lower = a(1:first-1, :)
+                a_upper = a(last+1:, :)
+                a(1:first-1, :) = 0.
+                a(last+1:, :) = 0.
+              case(2)
+                a_lower = a(:, 1:first-1)
+                a_upper = a(:, last+1:)
+                a(:, 1:first-1) = 0.
+                a(:, last+1:) = 0.
+              case default
+                error stop "gather_real_2D_array: invalid dim argument"
+            end select
+
+            call co_sum(a, result_image=result_image)
+
+            if (result_image /= me) then
+              select case(dim_)
+                case(1)
+                  a(1:first-1, :) = a_lower
+                  a(last+1:, :) = a_upper
+                case(2)
+                  a(:, 1:first-1) = a_lower
+                  a(:, last+1:) = a_upper
+                case default
+                  error stop "gather_real_2D_array: invalid dim argument"
+              end select
+            end if
+          end block
+        end if
+      end associate
+    end associate
+  end procedure
+
+end submodule data_partition_implementation
diff --git a/src/data-partition-interface.f90 b/src/data-partition-interface.f90
new file mode 100644
index 00000000..39ed0cb2
--- /dev/null
+++ b/src/data-partition-interface.f90
@@ -0,0 +1,64 @@
+module data_partition_interface
+  !! distribute data identification numbers across images such that the number of
+  !! items differs by at most 1 between any two images.
+  use iso_fortran_env, only : real64
+  implicit none
+
+  private
+  public :: data_partition
+
+  type data_partition
+    !! encapsulate a description of the data subset the executing image owns
+    private
+  contains
+    procedure, nopass :: define_partitions
+    procedure, nopass :: first
+    procedure, nopass :: last
+    procedure, nopass, private :: gather_real_2D_array, gather_real_1D_array
+    generic :: gather => gather_real_2D_array, gather_real_1D_array
+  end type
+
+  integer, allocatable :: first_datum(:), last_datum(:)
+
+  interface
+
+    module subroutine define_partitions(cardinality)
+      !! define the range of data identification numbers owned by the executing image
+      integer, intent(in) :: cardinality
+    end subroutine
+
+    pure module function first(image_number) result(first_index)
+      !! the result is the first identification number owned by the executing image
+      implicit none
+      integer, intent(in) :: image_number
+      integer first_index
+    end function
+
+    pure module function last(image_number) result(last_index)
+      !! the result is the last identification number owned by the executing image
+      implicit none
+      integer, intent(in) :: image_number
+      integer last_index
+    end function
+
+    !! Gathers are inherently expensive and are best used either
+    !! 1. Near the beginning/end of execution to amortize costs across an entire run or
+    !! 2. Temporarily while developing/debugging code.
+
+    module subroutine gather_real_1D_array( a, result_image, dim )
+      !! Gather the elements of an 1D array distributed along dimension dim onto result_image
+      real(real64), intent(inout) :: a(:)
+      integer, intent(in), optional :: result_image
+      integer, intent(in), optional :: dim
+    end subroutine
+
+    module subroutine gather_real_2D_array( a, result_image, dim )
+      !! Gather the elements of an 2D array distributed along dimension dim onto result_image
+      real(real64), intent(inout) :: a(:,:)
+      integer, intent(in), optional :: result_image
+      integer, intent(in), optional :: dim
+    end subroutine
+
+  end interface
+
+end module data_partition_interface
diff --git a/src/emulated_intrinsics_implementation.F90 b/src/emulated_intrinsics_implementation.F90
index cd7ee5cc..0a5a4a4d 100644
--- a/src/emulated_intrinsics_implementation.F90
+++ b/src/emulated_intrinsics_implementation.F90
@@ -12,13 +12,19 @@
 
   module procedure co_all
       call co_reduce(boolean, both)
+#ifdef FORD
+  end procedure
+#else
   contains
+#endif
     pure function both(lhs,rhs) result(lhs_and_rhs)
       logical, intent(in) :: lhs,rhs
       logical lhs_and_rhs
       lhs_and_rhs = lhs .and. rhs
     end function
+#ifndef FORD
   end procedure
+#endif
 
 #ifdef COMPILER_LACKS_COLLECTIVE_SUBROUTINES
   module procedure co_sum_integer
diff --git a/tests/data_partition_test.f90 b/tests/data_partition_test.f90
new file mode 100644
index 00000000..16f3e798
--- /dev/null
+++ b/tests/data_partition_test.f90
@@ -0,0 +1,203 @@
+module data_partition_test
+   !! author: Damian Rouson
+   !!
+   !! summary: verify data partitioning across images and data gathering
+   use vegetables, only: &
+     result_t, input_t, integer_input_t, test_item_t, & ! types
+     describe, it, assert_equals, assert_that, example  ! functions
+   use data_partition_interface, only : data_partition
+   use iso_fortran_env, only : real64
+   implicit none
+
+   private
+   public :: test_data_partition
+
+   type(data_partition) partition
+   integer, parameter :: num_particles=31, gatherer=1, num_steps=9, dummy=0
+
+contains
+
+  function test_data_partition() result(tests)
+    type(test_item_t) tests
+
+    integer iteration
+
+    call partition%define_partitions( cardinality=num_particles)
+
+    associate( me=>this_image() )
+      associate( my_first=>partition%first(me), my_last=>partition%last(me) )
+        tests = describe( &
+         "data_partition class", &
+         [it( &
+           "partitions data in nearly even blocks", &
+           verify_block_partitioning), &
+          it( &
+           "all data partitioned across all images without data loss", &
+           verify_all_particles_partitioned), &
+          it( &
+           "1D real array gathered on all images", &
+           [example(integer_input_t(dummy)), example(integer_input_t(dummy))], &
+           verify_all_gather_1D_real_array), &
+          it( &
+           "dimension 1 of 2D real array gathered on all images witout dim argument", &
+           [example(integer_input_t(dummy)), example(integer_input_t(dummy))], &
+           verify_all_gather_2D_real_array), &
+          it( &
+           "dimension 1 of 2D real array gathered on all images with dim argument", &
+           [example(integer_input_t(dummy)), example(integer_input_t(dummy))], &
+           verify_all_gather_2D_real_array_dim1), &
+          it( &
+           "dimension 1 of 2D real array gathered onto result_image with dim argument", &
+           [example(integer_input_t(dummy)), example(integer_input_t(dummy))], &
+           verify_gather_2D_real_array_dim1)])
+
+      end associate
+    end associate
+  end function
+
+  function verify_block_partitioning() result(result_)
+    !! Verify that the data is partitioned across images evenly to
+    !! within a difference of one datum between any two images.
+    type(data_partition) partition
+    type(result_t) result_
+    integer my_particles
+
+    associate( me=>this_image() )
+      associate( my_first=>partition%first(me), my_last=>partition%last(me) )
+        my_particles = my_last - my_first + 1
+        associate( ni=>num_images() )
+          associate( quotient=>num_particles/ni, remainder=>mod(num_particles,ni)  )
+            result_ = assert_equals( quotient + merge(1, 0, me<=remainder), my_particles, "block distribution" )
+          end associate
+        end associate
+      end associate
+    end associate
+  end function
+
+  function verify_all_particles_partitioned() result(result_)
+    !! Verify that the number of particles on each image sums to the
+    !! total number of particles distributed.
+    type(data_partition) partition
+    type(result_t) result_
+    integer particles
+
+    associate(me => this_image())
+      associate( my_first=>partition%first(me), my_last=>partition%last(me) )
+        particles = my_last - my_first + 1
+        call co_sum(particles)
+        result_ = assert_equals(num_particles, particles, "all particles distributed" )
+      end associate
+    end associate
+  end function
+
+ function verify_all_gather_1D_real_array(unused) result(result_)
+   type(data_partition) partition
+   class(input_t), intent(in) :: unused
+   type(result_t) result_
+   real(real64) :: particle_scalar(num_particles)
+   real(real64), parameter :: junk=-12345._real64, expected=1._real64
+
+   associate( no_op => unused) ! eliminate unused-variable warning
+   end associate
+
+   associate(me => this_image())
+     associate( first=>partition%first(me), last=>partition%last(me) )
+
+       particle_scalar(first:last) = expected !! values to be gathered
+       particle_scalar(1:first-1)  = junk !! values to be overwritten by the gather
+       particle_scalar(last+1:)  = junk !! values to be overwritten by the gather
+
+       call partition%gather(particle_scalar)
+
+       result_ = assert_that( all(particle_scalar==expected), "real 1D array all-gathered" )
+
+     end associate
+   end associate
+ end function
+
+ function verify_all_gather_2D_real_array(unused) result(result_)
+   class(input_t), intent(in) :: unused
+   type(data_partition) partition
+   type(result_t) result_
+   integer, parameter :: vec_space_dim=3
+   real(real64) particle_vector(vec_space_dim, num_particles)
+   real(real64), parameter :: junk=-12345._real64, expected=1._real64
+
+   associate( no_op => unused) ! eliminate unused-variable warning
+   end associate
+
+   associate(me => this_image())
+     associate( first=>partition%first(me), last=>partition%last(me) )
+
+       particle_vector(:, first:last) = expected !! values to be gathered
+       particle_vector(:, 1:first-1)  = junk !! values to be overwritten by the gather
+       particle_vector(:, last+1:)  = junk !! values to be overwritten by the gather
+
+       call partition%gather(particle_vector)
+
+       result_ = assert_that(all(particle_vector==expected), "real 2D array all-gathered implicitly along dimension 1" )
+
+     end associate
+   end associate
+ end function
+
+ function verify_all_gather_2D_real_array_dim1(unused) result(result_)
+   class(input_t), intent(in) :: unused
+   type(data_partition) partition
+   type(result_t) result_
+   integer, parameter :: vec_space_dim=3
+   real(real64) :: vector_transpose(num_particles, vec_space_dim)
+   real(real64), parameter :: junk=-12345._real64, expected=1._real64
+
+   associate( no_op => unused) ! eliminate unused-variable warning
+   end associate
+
+   associate(me => this_image())
+     associate( first=>partition%first(me), last=>partition%last(me) )
+
+       vector_transpose(first:last, :) = expected !! values to be gathered
+       vector_transpose(1:first-1, :)  = junk !! values to be overwritten by the gather
+       vector_transpose(last+1:, :)  = junk !! values to be overwritten by the gather
+
+       call partition%gather( vector_transpose, dim=1)
+
+       result_ = assert_that(all(vector_transpose==expected), "vector_transpose gathered explicitly along dimension 1" )
+
+      end associate
+    end associate
+ end function
+
+ function verify_gather_2D_real_array_dim1(unused) result(result_)
+   class(input_t), intent(in) :: unused
+   type(data_partition) partition
+   type(result_t) result_
+   integer, parameter :: vec_space_dim=3
+   real(real64) :: vector_transpose(num_particles, vec_space_dim)
+   real(real64), parameter :: junk=-12345._real64, expected=1._real64
+
+   associate( no_op => unused) ! eliminate unused-variable warning
+   end associate
+
+   associate(me => this_image())
+     associate( first=>partition%first(me), last=>partition%last(me) )
+
+       vector_transpose(first:last, :) = expected !! values to be gathered
+       vector_transpose(1:first-1, :)  = junk !! values to be overwritten by the gather
+       vector_transpose(last+1:, :)  = junk !! values to be overwritten by the gather
+
+       call partition%gather( vector_transpose, result_image=gatherer, dim=1)
+
+       if (me==gatherer) then
+         result_ = assert_that(all(vector_transpose==expected), "all( particle_vector==expected)")
+       else
+         result_ = &
+           assert_that(all(vector_transpose(1:first-1,:)==junk), "lower transpose data unchanged)") .and. &
+           assert_that(all(vector_transpose(first:last,:)==expected), "expected transpose data gathered") .and. &
+           assert_that(all(vector_transpose(last+1:,:)==junk), "upper transpose data unchanged)" )
+       end if
+
+     end associate
+   end associate
+ end function
+
+end module data_partition_test
diff --git a/tests/main.f90 b/tests/main.f90
index cdaa6e74..543d06d8 100644
--- a/tests/main.f90
+++ b/tests/main.f90
@@ -6,18 +6,21 @@ program main
 contains
     subroutine run()
         use collective_subroutines_test, only: &
-                collective_subroutines_co_sum => test_co_sum, &
-                collective_subroutines_co_all => test_co_all
+                collective_subroutines_co_all => test_co_all, &
+                collective_subroutines_co_sum => test_co_sum
+        use data_partition_test, only: &
+                data_partition_data_partition => test_data_partition
         use single_image_intrinsics_test, only: &
                 single_image_intrinsics_findloc => test_findloc
         use vegetables, only: test_item_t, test_that, run_tests
 
         type(test_item_t) :: tests
-        type(test_item_t) :: individual_tests(3)
+        type(test_item_t) :: individual_tests(4)
 
-        individual_tests(1) = collective_subroutines_co_sum()
-        individual_tests(2) = collective_subroutines_co_all()
-        individual_tests(3) = single_image_intrinsics_findloc()
+        individual_tests(1) = collective_subroutines_co_all()
+        individual_tests(2) = collective_subroutines_co_sum()
+        individual_tests(3) = data_partition_data_partition()
+        individual_tests(4) = single_image_intrinsics_findloc()
         tests = test_that(individual_tests)
 
         call run_tests(tests)