!
  !-----------------------------------------------------------------------
  subroutine elphon_shuffle ( iq_irr, nqc_irr, iq, gmapsym, eigv, isym, invs0, xq0, timerev )
  !-----------------------------------------------------------------------
  !
  ! Electron-phonon calculation from data saved in fildvscf
  ! Shuffle2 mode (shuffle on electrons + load all phonon q's)
  !
  ! no ultrasoft yet
  !
  !-----------------------------------------------------------------------
  !
#include "f_defs.h"
  !
#ifdef __PARA
  use para
  use mp, ONLY : mp_barrier
  USE mp_global,   ONLY : nproc, my_pool_id, nproc_pool
#endif
  !
  USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau
  !
  use io_global, only : stdout
  use pwcom
  USE kinds, only : DP
  use phcom
  use el_phon
  implicit none
  complex(kind=DP), parameter :: czero = (0.d0,0.d0), cone = (1.d0,0.d0)
  !
  integer :: irr, imode0, ipert, is, iq, iq_irr, nqc_irr
  ! counter on the representations
  ! counter on the modes
  ! the change of Vscf due to perturbations
  ! the current qpoint in the uniform grid
  ! the current ireducible qpoint
  ! the total number of irreducible qpoints in the list
  complex(kind=DP), pointer :: dvscfin(:,:,:), dvscfins (:,:,:)
  logical :: timerev
  !  true if we are using time reversal

  !
  integer :: tmp_pool_id, nkl, nkr, ik0, iks, ik, ibnd, jbnd, mu, nu
  real(kind=DP) :: tmp
  character(len=30) filepmat
  complex(kind=8) :: eye(nmodes,nmodes) ! debug
  !
  integer :: gmapsym ( ngm, 48 ), isym, invs0 (48)
  ! the correspondence G-->S(G)
  ! the symmetry which generates the current q in the star
  ! the index of the inverse operations
  complex(kind=DP) :: eigv (ngm, 48)
  ! e^{ iGv} for 1...nsym ( v the fractional translation)
  real(kind=DP) :: xq0(3)
  ! the first q in the star (cartesian)

  !
  call start_clock ('elphon')
  !
  ! tshuffle2 implies tshuffle
  !
  if (tshuffle2) tshuffle = .true.
  !
  ik0 = 0
  tmp_pool_id = 0
  !
#ifdef __PARA
  !
  npool = nproc / nproc_pool
  if (npool.gt.1) then
    !
    ! number of kpoint blocks, kpoints per pool and reminder
    nkbl  = nkstot / kunit
    nkl   = kunit * ( nkbl / npool )
    nkr   = ( nkstot - nkl * npool ) / kunit
    ! the reminder goes to the first nkr pools
    if ( my_pool_id < nkr ) nkl = nkl + kunit
    !
    iks = nkl * my_pool_id + 1
    if ( my_pool_id >= nkr ) iks = iks + nkr * kunit
    !
    !  the index of the first k point block in this pool - 1
    !  (I will need the index of ik, not ikk)
    !
    ik0 = ( iks - 1 ) / kunit
    !
  endif
  !
#endif
  !    

  !
  ! read Delta Vscf and calculate electron-phonon coefficients
  !
  imode0 = 0
  do irr = 1, nirr
     allocate (dvscfin ( nrxx , nspin , npert(irr)) )
     !
     !   read the <prefix>.dvscf_q[iq] files
     !
     dvscfin = (0.d0,0.d0)
     if (mypool.eq.1) then
        do ipert = 1, npert (irr)
           call readdvscf ( dvscfin(1,1,ipert), imode0 + ipert, iq_irr, nqc_irr )
        end do
     endif
     call poolreduce ( 2 * npert(irr) * nspin * nrxx, dvscfin )
     !
     if (doublegrid) then
        allocate (dvscfins ( nrxxs , nspin , npert(irr)) )
        do is = 1, nspin
           do ipert = 1, npert(irr)
              call cinterpolate (dvscfin(1,is,ipert),dvscfins(1,is,ipert),-1)
           enddo 
        enddo
     else
        dvscfins => dvscfin
     endif
     !
     call elphel2_shuffle (npert (irr), imode0, dvscfins, gmapsym, eigv, isym, invs0, xq0, timerev)
     !
     imode0 = imode0 + npert (irr)
     if (doublegrid) deallocate (dvscfins)
     deallocate (dvscfin)
  enddo
  !
  !  the output e-p matrix in the pattern representation
  !  must be transformed in the cartesian basis
  !  epmat_{CART} = conjg ( U ) * epmat_{PATTERN}
  !
  !  note it is not U^\dagger ! Have a look to symdyn_munu.f90 
  !  for comparison
  !
  do ibnd = 1, nbnd
    do jbnd = 1, nbnd
      do ik = 1, nksq
        ! 
        call zgemv ('n', nmodes, nmodes, cone, conjg ( u ), nmodes, &
          el_ph_mat (ibnd,jbnd,ik,:), 1, czero, epmatq (ibnd,jbnd,ik,:,iq), 1 )
        !
      enddo
    enddo
  enddo
  !
  call stop_clock ('elphon')
  return
end subroutine elphon_shuffle
!