!
  !---------------------------------------------------------------------------------
  subroutine dynbloch2wan ( nmodes, nq, xk, dynq, nrr, irvec, wslen, rdw )
  !---------------------------------------------------------------------------------
  !
  !  From the Dynamical Matrix in Bloch representation (coarse mesh), 
  !  find the corresponding matrix in Wannier representation 
  !
  !  input  : 
  !
  !  output : 
  !
  !  NOTA BENE: it seems to be very important that the matrix is kept real
  !  physically these are truely the interatomic force constants.
  !  If you use a complex matrix instead, you may get some spurious 
  !  oscillations when you interpolate the phonon dispersions.
  !
  !  Note also that the acoustic sum rule for the q=0 case has been imposed
  !  already in readmat_shuffle.f90
  !
  !
  !  Feliciano Giustino, UCB
  !
  !--------------------------------------------------------------------------------
  !
#include "f_defs.h"
  USE kinds, only : DP
  use pwcom, only : at, bg, celldm
  use control_flags, ONLY : iverbosity
#ifdef __PARA
  use para
  USE io_global, ONLY : ionode_id
  USE mp_global, ONLY : mpime
  USE mp, ONLY : mp_barrier
#endif
  implicit none
  !
  !  input variables
  !
  integer :: nmodes, nq, nrr, irvec (3, nrr)
  ! number of branches
  ! number of qpoints
  ! number of WS points and coordinates
  complex(kind=DP) :: dynq (nmodes, nmodes, nq)
  ! dynamical matrix in bloch representation (Cartesian coordinates)
  real(kind=DP) :: xk (3, nq), wslen (nrr) 
  ! kpoint coordinates (cartesian in units of 2piba)
  ! WS vectors length (alat units)
  !
  !  output variables
  !
!@  real(kind=DP) :: rdw (nmodes, nmodes, nrr)  
complex(kind=DP) :: rdw (nmodes, nmodes, nrr)  
  !
  ! work variables 
  !
  real(kind=DP), parameter :: bohr2ang = 0.5291772108, twopi = 6.28318530717959
  complex(kind=DP), parameter :: ci = (0.d0,1.d0), czero = (0.d0, 0.d0)
  integer :: ik, ikk, ikq, ibnd, jbnd, imode0, ipol, nsize, rest, &
             ik0, ind, ix, jx, ir, mbnd, pbnd, i
  real(kind=DP) :: rdotk, tmp
  complex(kind=DP) :: cfac, ctmp
  !
  !----------------------------------------------------------
  !  Fourier transform to go into Wannier basis
  !----------------------------------------------------------
  !
  !  D (R) = (1/nk) sum_k e^{-ikR} D (k)
  !  rdw (nmodes, nmodes, ir) is D (R)
  !
  ! bring xk in crystal coordinates
  !
  call cryst_to_cart (nq, xk, at, -1)
  !
  rdw ( :, :, :) = 0.d0
  ! 
  do ir = 1, nrr
    !
    do ik = 1, nq
       !
       rdotk = twopi * dot_product( xk ( :, ik), float(irvec( :, ir) ))
       cfac = exp( -ci*rdotk ) / float(nq)
 rdw ( :, :, ir ) = rdw ( :, :, ir ) +  cfac * dynq ( :, :, ik ) 
!@     rdw ( :, :, ir ) = rdw ( :, :, ir ) + real ( cfac * dynq ( :, :, ik ) ) 
       !                                     ^^^^
       !                                   note this
       !
    enddo
    !
  enddo
  !
  ! bring xk back into cart coord
  !
  call cryst_to_cart (nq, xk, bg, 1)
  !
!@  if (iverbosity.eq.1) then
    !
    !  check spatial decay of dynamical matrix in Wannier basis
    !  the unit in r-space is angstrom, and I am plotting
    !  the matrix for the first mode only
    !
#ifdef __PARA
    if (mpime.eq.ionode_id) then
#endif
!@    write(302, '(/3x,a/)') 'Spatial decay of Dynamical matrix in Wannier basis'
      do ir = 1, nrr
        !
        tmp =  maxval ( abs( rdw (:,:,ir)) )
        write(302, *) wslen(ir) * celldm (1) * bohr2ang, tmp
        !
      enddo
#ifdef __PARA
    endif
    call mp_barrier()
#endif
!@  endif
  !
  return
  end subroutine dynbloch2wan
  !-----------------------------------------------------