def gather_ps(rank, size, comm, k_allmodels, P21_allmodels, PHII_allmodels,
              first_snap_allmodels, last_snap_allmodels):
    """
    Gathers the power spectra calculated on each processor onto the root rank.
    Each rank calculates the spectra of only a subset of snapshots so here we
    gather the spectra of all models and snapshots onto root rank. 

    Parameters
    ----------

    rank : Integer
        This processor rank.

    size : Integer
        The total number of processors executing the pipeline.

    comm : Class ``mpi4py.MPI.Intracomm``
        The ``mpi4py`` communicator.

    k_allmodels : 3D nested list of floats. Outer length is number of models,
                  next is number of snapshots processed by this processor and
                  final is the number of wavenumber bins.
        Wavenumber the spectra are binned on (units of Mpc/h). Processor
        unique.

    P21_allmodels, PHII_allmodels : 3D nested lists of floats. Dimensions are
                                    identical to ``k_allmodels``.
        The 21cm and HII power spectra for each model at each snapshot.
        Processor unique.

    first_snap_allmodels, last_snap_allmodels : List of integers. Length is
                                                number of models.
        The first and last snapshot that defines the snapshot range that
        ``cifog`` was run on. 

    Returns
    ---------

    For all non-zero ranks, the returns are:
        None, None, None.

    k_master, P21_master, PHII_master : 3D nested lists of floats. Dimensions
                                        are identical to `k_allmodels` except
                                        the 2nd axis length is the snapshot
                                        range for that particular model (not a
                                        subset).
        The wavenumber bins, 21cm and HII power spectra for each model for all
        snapshots. 
    """

    def generate_tag(rank):
        tag = int(rank*100)

        return tag

    # Rank 0 will gather the wavenumber bins/power spectra from all other
    # ranks. 
    if rank == 0:
        k_master = []
        P21_master = []
        PHII_master = []

        # Go through each model. 
        for model_number in range(len(k_allmodels)):

            k_master.append([])
            P21_master.append([])
            PHII_master.append([])

            model_k = k_allmodels[model_number]
            model_P21 = P21_allmodels[model_number]
            model_PHII = PHII_allmodels[model_number]

            num_snaps = last_snap_allmodels[model_number] - \
                        first_snap_allmodels[model_number]
            rank_count = 0
            my_count = 0

            # Then go through each snapshot.
            # In the main data loop (``generate_data()``) the snapshots are
            # scatter sequentially. Hence when we gather, we get snap0 from
            # rank 0, snap1 from rank 1 etc. So we increase rank_count for each
            # snapshot and then reset it when we reach `size`.
            for snap_idx in range(num_snaps):

                if rank_count == 0:
                    this_k = model_k[my_count] 
                    this_P21 = model_P21[my_count] 
                    this_PHII = model_PHII[my_count] 
                    my_count += 1
                else:
                    # Each rank will use a unique tag.
                    tag = generate_tag(rank_count) 

                    # Then the tag is offset for each data array. 
                    this_k = comm.recv(source = rank_count,
                                       tag = tag)
                    this_P21 = comm.recv(source = rank_count,
                                         tag = tag+1)
                    this_PHII = comm.recv(source = rank_count,
                                          tag = tag+2)

                # Now we have the data, append it to the master.
                k_master[model_number].append(this_k)
                P21_master[model_number].append(this_P21)
                PHII_master[model_number].append(this_PHII)

                rank_count += 1
                if rank_count == size:
                    rank_count = 0

            # Snapshot Loop.
        # Model Loop.

        return k_master, P21_master, PHII_master

    else:

        # For all other ranks, go through the power spectra it calculated and
        # send it back to the root rank.
        for model_number in range(len(k_allmodels)):
            for idx in range(len(P21_allmodels[model_number])):

                tag = generate_tag(rank) 

                k_this_idx = k_allmodels[model_number][idx]
                P21_this_idx = P21_allmodels[model_number][idx]
                PHII_this_idx = PHII_allmodels[model_number][idx]

                comm.send(k_this_idx, dest = 0, tag = tag)
                comm.send(P21_this_idx, dest = 0, tag = tag+1)
                comm.send(PHII_this_idx, dest = 0, tag = tag+2)

        # Non-zero ranks return junk.
        return None, None, None