!
  !-----------------------------------------------------------------------
  subroutine readmat_shuffle ( q, dyn, w2, iq, nqc)
  !-----------------------------------------------------------------------
  !
#include "f_defs.h"
  USE kinds, only : DP
  use io_files, only : prefix, tmp_dir
  USE cell_base, ONLY : ibrav, celldm, omega
  USE ions_base, ONLY : amass, tau, nat, ntyp => nsp, ityp
  USE el_phon
  !
  implicit none
  !
  ! Input
  !
  integer :: iq, nqc
  !  the current q point
  !  the total number of qpoints in the list
  !
  ! output
  !
  complex(kind=DP) :: dyn (3 * nat, 3 * nat)
  real(kind=DP) :: w2 (3 * nat)
  !
  ! local (control variables)
  !
  integer :: iudyn
  real(kind=DP) :: q (3), eps
  integer :: ntyp_, nat_, ibrav_, ityp_, ios
  real(kind=DP) :: celldm_ (6), amass_, tau_ (3), q_ (3)
  !
  ! local
  !
  real(kind=DP) :: dynr (2, 3, nat, 3, nat), massfac, sumr(2)
  character(len=80) :: line
  character(len=3)  :: atm
  integer :: nt, na, nb, naa, nbb, nu, mu, i, j, ipol, jpol
  !
  character (len=256) :: tempfile
  character (len=3) :: filelab
  ! file label
  !
  eps = 1.d-5
  !
  !  the call to set_ndnmbr is just a trick to get quickly
  !  a file label by exploiting an existing subroutine
  !  (if you look at the sub you will find that the original
  !  purpose was for pools and nodes)
  !
  call set_ndnmbr ( 0, iq, 1, nqc, filelab)
  tempfile = trim(tmp_dir) // trim(prefix) // '.dyn_q' // filelab
  !
  iudyn = 81
  open (unit = iudyn, file = tempfile, status = 'unknown', err = 100, iostat = ios)
100  call errore ('openfilq', 'opening file'//tempfile, abs (ios) )
  rewind (iudyn)
  !
  !
  !
  read (iudyn, '(a)') line
  read (iudyn, '(a)') line
  read (iudyn, * ) ntyp_, nat_, ibrav_, celldm_
  if (ntyp.ne.ntyp_.or.nat.ne.nat_.or.ibrav_.ne.ibrav.or.abs ( &
     celldm_ (1) - celldm (1) ) .gt.1.0d-5) call errore ('readmat', &
     'inconsistent data', 1)
  do nt = 1, ntyp
     read (iudyn, * ) i, atm, amass_
     if (nt.ne.i.or.abs (amass_ - amass (nt) ) .gt.1.0d-5) call errore ( &
          'readmat', 'inconsistent data', 1 + nt)
  enddo
  do na = 1, nat
     read (iudyn, * ) i, ityp_, tau_
     if (na.ne.i.or.ityp_.ne.ityp (na) ) call errore ('readmat', &
        'inconsistent data', 10 + na)
  enddo
  !
  !
  read (iudyn, '(///a)') line
  read (line (11:80), * ) (q_ (i), i = 1, 3)
  if ( abs(q_ (1)-q(1)).gt.eps .or. &
       abs(q_ (2)-q(2)).gt.eps .or. &
       abs(q_ (3)-q(3)).gt.eps ) call errore &
     ('readmat', 'wrong qpoint', 1)
  read (iudyn, '(a)') line
  !
  !
  do na = 1, nat
     do nb = 1, nat
        read (iudyn, * ) naa, nbb
        if (na.ne.naa.or.nb.ne.nbb) call errore &
             ('readmat', 'error reading file', nb)
        read (iudyn, * ) &
            ((dynr (1,i,na,j,nb), dynr (2,i,na,j,nb), j = 1, 3), i = 1, 3)
     enddo
  enddo
  !
  if ( abs(q(1)).lt.eps .and. abs(q(2)).lt.eps .and. abs(q(3)).lt.eps ) then
    !
    ! in the case q=0 we impose the acoustic sum rule
    ! [Gonze and Lee, PRB 55, 10361 (1998), Eq. (45)]
    !
    write(6,'(8x,a)') 'Imposing acoustic sum rule on the dynamical matrix'
    !
    do na = 1,nat
     do ipol = 1,3
      do jpol = ipol,3
         !
         sumr(1) = sum ( dynr (1,ipol,na,jpol,:) )
         sumr(2) = sum ( dynr (2,ipol,na,jpol,:) )
         !
         dynr (:,ipol,na,jpol,na) = dynr (:,ipol,na,jpol,na) - sumr
         !
      end do
     end do
    end do
    !
  endif
  !
  !  fill the two-indices dynamical matrix in cartesian coordinates
  !
  do na = 1,nat
   do nb = 1,nat
      do ipol = 1,3
       do jpol = 1,3
         !
         mu = (na-1)*3+ipol
         nu = (nb-1)*3+jpol
         dynq ( mu, nu, iq) = dcmplx ( dynr (1,ipol,na,jpol,nb), dynr (2,ipol,na,jpol,nb) )
         !
       end do
      end do
   end do
  end do
  !
  close(iudyn)
  !
  return
  end subroutine readmat_shuffle
  !