setup_face_halo Subroutine

   subroutine setup_face_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, face_owner_rank, my_index, first, last
      integer :: total_recv, total_send, total_requests

      allocate(need(mesh%nfaces, 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)
               face_owner_rank = flow%cell_owner(mesh%faces(f)%owner)
               if (face_owner_rank /= q - 1) need(f, q) = .true.
            end do
         end do
      end do

      my_index = flow%rank + 1
      do f = 1, mesh%nfaces
         if (.not. need(f, my_index)) cycle
         face_owner_rank = flow%cell_owner(mesh%faces(f)%owner)
         recv_counts_all(face_owner_rank + 1) = recv_counts_all(face_owner_rank + 1) + 1
      end do

      do q = 1, flow%nprocs
         if (q == my_index) cycle
         do f = 1, mesh%nfaces
            if (need(f, q) .and. flow%cell_owner(mesh%faces(f)%owner) == 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%nface_recv_ranks, flow%face_recv_ranks, &
                              flow%face_recv_counts, flow%face_recv_displs)
      call pack_rank_metadata(send_counts_all, send_displs_all, &
                              flow%nface_send_ranks, flow%face_send_ranks, &
                              flow%face_send_counts, flow%face_send_displs)

      allocate(flow%face_recv_faces(total_recv))
      allocate(flow%face_send_faces(total_send))
      allocate(flow%face_recvbuf(max(1, total_recv)))
      allocate(flow%face_sendbuf(max(1, total_send)))

      recv_next = recv_displs_all
      do f = 1, mesh%nfaces
         if (.not. need(f, my_index)) cycle
         face_owner_rank = flow%cell_owner(mesh%faces(f)%owner) + 1
         recv_next(face_owner_rank) = recv_next(face_owner_rank) + 1
         flow%face_recv_faces(recv_next(face_owner_rank)) = f
      end do

      send_next = send_displs_all
      do q = 1, flow%nprocs
         if (q == my_index) cycle
         do f = 1, mesh%nfaces
            if (.not. need(f, q)) cycle
            if (flow%cell_owner(mesh%faces(f)%owner) /= flow%rank) cycle
            send_next(q) = send_next(q) + 1
            flow%face_send_faces(send_next(q)) = f
         end do
      end do

      total_requests = flow%nface_recv_ranks + flow%nface_send_ranks
      allocate(flow%face_requests(max(1, total_requests)))
      flow%face_recvbuf = zero
      flow%face_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_face_halo