setup_cell_halo Subroutine

   subroutine setup_cell_halo(mesh, flow)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(inout) :: flow
      logical, allocatable :: need(:,:)
      integer, allocatable :: recv_counts_all(:), recv_displs_all(:), recv_next(:)
      integer, allocatable :: send_counts_all(:), send_displs_all(:), send_next(:)
      integer :: q, c, lf, f, nb, owner, my_index, first, last
      integer :: total_recv, total_send, total_requests

      allocate(need(mesh%ncells, flow%nprocs))
      allocate(recv_counts_all(flow%nprocs), recv_displs_all(flow%nprocs), recv_next(flow%nprocs))
      allocate(send_counts_all(flow%nprocs), send_displs_all(flow%nprocs), send_next(flow%nprocs))

      need = .false.
      recv_counts_all = 0
      send_counts_all = 0

      do q = 1, flow%nprocs
         first = flow%gather_firsts(q)
         last = first + flow%gather_counts(q) - 1
         do c = first, last
            do lf = 1, mesh%ncell_faces(c)
               f = mesh%cell_faces(lf, c)
               nb = mesh_neighbor_for_cell(mesh, f, c)
               if (nb <= 0) cycle
               owner = flow%cell_owner(nb)
               if (owner /= q - 1) need(nb, q) = .true.
            end do
         end do
      end do

      my_index = flow%rank + 1
      do c = 1, mesh%ncells
         if (.not. need(c, my_index)) cycle
         owner = flow%cell_owner(c)
         recv_counts_all(owner + 1) = recv_counts_all(owner + 1) + 1
      end do

      do q = 1, flow%nprocs
         if (q == my_index) cycle
         do c = 1, mesh%ncells
            if (need(c, q) .and. flow%cell_owner(c) == flow%rank) then
               send_counts_all(q) = send_counts_all(q) + 1
            end if
         end do
      end do

      call prefix_counts(recv_counts_all, recv_displs_all)
      call prefix_counts(send_counts_all, send_displs_all)

      total_recv = sum(recv_counts_all)
      total_send = sum(send_counts_all)

      call pack_rank_metadata(recv_counts_all, recv_displs_all, &
                              flow%ncell_recv_ranks, flow%cell_recv_ranks, &
                              flow%cell_recv_counts, flow%cell_recv_displs)
      call pack_rank_metadata(send_counts_all, send_displs_all, &
                              flow%ncell_send_ranks, flow%cell_send_ranks, &
                              flow%cell_send_counts, flow%cell_send_displs)

      allocate(flow%cell_recv_cells(total_recv))
      allocate(flow%cell_send_cells(total_send))
      allocate(flow%ghost_cells(total_recv))
      allocate(flow%cell_recvbuf(max(1, total_recv * flow%cell_halo_max_components)))
      allocate(flow%cell_sendbuf(max(1, total_send * flow%cell_halo_max_components)))

      recv_next = recv_displs_all
      do c = 1, mesh%ncells
         if (.not. need(c, my_index)) cycle
         owner = flow%cell_owner(c) + 1
         recv_next(owner) = recv_next(owner) + 1
         flow%cell_recv_cells(recv_next(owner)) = c
      end do
      flow%ghost_cells = flow%cell_recv_cells

      send_next = send_displs_all
      do q = 1, flow%nprocs
         if (q == my_index) cycle
         do c = 1, mesh%ncells
            if (.not. need(c, q)) cycle
            if (flow%cell_owner(c) /= flow%rank) cycle
            send_next(q) = send_next(q) + 1
            flow%cell_send_cells(send_next(q)) = c
         end do
      end do

      total_requests = flow%ncell_recv_ranks + flow%ncell_send_ranks
      allocate(flow%cell_requests(max(1, total_requests)))
      flow%cell_recvbuf = zero
      flow%cell_sendbuf = zero

      deallocate(need)
      deallocate(recv_counts_all, recv_displs_all, recv_next)
      deallocate(send_counts_all, send_displs_all, send_next)
   end subroutine setup_cell_halo