flow_mpi_initialize Subroutine

   subroutine flow_mpi_initialize(mesh, flow, comm_parent, max_gather_components)
      type(mesh_t), intent(in) :: mesh
      type(flow_mpi_t), intent(inout) :: flow
      type(MPI_Comm), intent(in) :: comm_parent
      integer, intent(in), optional :: max_gather_components

      integer :: ierr, base, rem, c
      integer :: max_components

      call flow_mpi_finalize(flow)

      call MPI_Comm_dup(comm_parent, flow%comm, ierr)
      call check_mpi(ierr, 'MPI_Comm_dup flow')

      call MPI_Comm_rank(flow%comm, flow%rank, ierr)
      call check_mpi(ierr, 'MPI_Comm_rank flow')

      call MPI_Comm_size(flow%comm, flow%nprocs, ierr)
      call check_mpi(ierr, 'MPI_Comm_size flow')

      ! Contiguous cell range calculation:
      ! The total ncells are split into roughly equal blocks. If ncells is not 
      ! divisible by nprocs, the first 'rem' ranks get an extra cell to ensure 
      ! all cells are covered.
      base = mesh%ncells / flow%nprocs
      rem = mod(mesh%ncells, flow%nprocs)

      if (flow%rank < rem) then
         flow%nlocal = base + 1
         flow%first_cell = flow%rank * (base + 1) + 1
      else
         flow%nlocal = base
         flow%first_cell = rem * (base + 1) + (flow%rank - rem) * base + 1
      end if

      ! The global cell indices owned by this rank are [first_cell, last_cell].
      flow%last_cell = flow%first_cell + flow%nlocal - 1

      allocate(flow%owned(mesh%ncells))
      flow%owned = .false.

      do c = flow%first_cell, flow%last_cell
         if (c >= 1 .and. c <= mesh%ncells) flow%owned(c) = .true.
      end do

      max_components = 4
      if (present(max_gather_components)) max_components = max(max_components, max_gather_components)
      flow%cell_halo_max_components = max(1, max_components)

      call setup_owned_gather(mesh, flow, max_gather_components)
      call setup_cell_owners(mesh, flow)
      call setup_owned_faces(mesh, flow)
      call setup_cell_halo(mesh, flow)
      call setup_face_halo(mesh, flow)
   end subroutine flow_mpi_initialize