communicate_threshold_based_SIs Subroutine

    subroutine communicate_threshold_based_SIs(si_buf, ref_buf, refs_found, all_refs_found)
        use semi_stoch_procs, only: return_largest_indices
        implicit none
        integer(n_int), intent(out) :: si_buf(0:NIfTot, maxNRefs)
        integer(n_int), intent(in) :: ref_buf(0:NIfTot, maxNRefs)
        integer, intent(in) :: refs_found
        integer, intent(out) :: all_refs_found
        integer(n_int), allocatable :: mpi_buf(:, :)
        integer(MPIArg) :: refs_found_per_proc(0:nProcessors - 1), refs_displs(0:nProcessors - 1)
        integer(MPIArg) :: mpi_refs_found
        integer :: ierr, i
        integer :: largest_inds(maxNRefs)
        real(dp), allocatable :: buf_signs(:)
        real(dp) :: tmp_sgn(lenof_sign)

        si_buf = 0
        ! here, we gather the potential SIs found by the procs and gather them
        ! Communicate the refs_found info
        mpi_refs_found = int(refs_found, MPIArg)
        call MPIAllGather(mpi_refs_found, refs_found_per_proc, ierr)
        ! total number of SI candiates
        all_refs_found = sum(refs_found_per_proc)
        refs_displs(0) = 0
        do i = 1, nProcessors - 1
            refs_displs(i) = refs_displs(i - 1) + refs_found_per_proc(i - 1)
        end do
        ! Store them on all processors
        allocate(mpi_buf(0:NIfTot, all_refs_found), stat=ierr)
        call MPIAllGatherV(ref_buf(0:NIfTot, 1:refs_found), mpi_buf, refs_found_per_proc, refs_displs)

        ! now, if we have have more than the maximum number of potential SIs, take the
        ! most populated only
        if (all_refs_found > maxNRefs) then
            ! get the indices of the largest elements in the communicated buffer
            ! first, extract the signs
            allocate(buf_signs(all_refs_found), stat=ierr)
            do i = 1, all_refs_found
                call extract_sign(mpi_buf(:, i), tmp_sgn)
                buf_signs(i) = sum(abs(tmp_sgn))
            end do
            ! then get the indices
            call return_largest_indices(maxNRefs, all_refs_found, buf_signs, largest_inds)
            ! and then copy those elements to the output buffer
            do i = 1, maxNRefs
                si_buf(0:NIfTot, i) = mpi_buf(0:NIfTot, largest_inds(i))
            end do
            if (iProcIndex == root .and. NoTypeN > 1) &
                write(stdout, '(A,I5,A,I5,A)') "In total ", all_refs_found, &
                " SIs were found, which is more than the maximum number of ", &
                maxNRefs, " - truncating"
            ! make it look to the outside as though maxNRefs were found
            all_refs_found = maxNRefs
        else
            ! else, take all elements
            si_buf(0:NIfTot, 1:all_refs_found) = mpi_buf(0:NIfTot, 1:all_refs_found)
        end if

        deallocate(mpi_buf)
    end subroutine communicate_threshold_based_SIs