! 
  !---------------------------------------------------------------------------
  subroutine ephwan2blochp ( nmodes, xxq, irvec, ndegen, nrr_q, cuf, epmatw, epmatf, nbnd, nrr_k )
  !---------------------------------------------------------------------------
  !
  ! even though this is for phonons, I use the same notations
  ! adopted for the electronic case (nmodes->nmodes etc)
  !
  !
#include "f_defs.h"
  USE kinds, only : DP
  implicit none
  !
  !  input variables
  !
  integer :: nmodes, nrr_q, irvec ( 3, nrr_q), ndegen (nrr_q), nbnd, nrr_k
  ! number of bands (possibly in tyhe optimal subspace)
  ! number of WS points
  ! coordinates of WS points
  ! degeneracy of WS points
  ! n of bands
  ! n of electronic WS points
  complex(kind=DP) :: epmatw ( nbnd, nbnd, nrr_k, nmodes, nrr_q), cuf (nmodes, nmodes)
  ! e-p matrix in Wanner representation
  ! rotation matrix U(k)
  real(kind=DP) :: xxq(3)
  ! kpoint for the interpolation (WARNING: this must be in crystal coord!)
  !
  !  output variables
  !
  complex(kind=DP) :: epmatf (nbnd, nbnd, nrr_k, nmodes)
  ! e-p matrix in Bloch representation, fine grid
  !
  ! work variables 
  !
  integer :: ibnd, jbnd, ir, ire
  real(kind=DP), parameter :: twopi = 6.28318530717959
  real(kind=DP) :: rdotk
  complex(kind=DP) :: cfac, eptmp( nbnd, nbnd, nrr_k, nmodes)
  complex(kind=DP), parameter :: ci = (0.d0,1.d0), &
     czero = (0.d0, 0.d0), cone = (1.d0, 0.d0)
  !
  !----------------------------------------------------------
  !  STEP 3: inverse Fourier transform of g to fine k mesh
  !----------------------------------------------------------
  !
  !  g~ (k') = sum_R 1/ndegen(R) e^{-ik'R} g (R)
  !
  !  g~(k') is epmatf (nmodes, nmodes, ik )
  !  every pool works with its own subset of k points on the fine grid
  !
  eptmp = czero
  !
  do ir = 1, nrr_q
     !
     ! note xxq is assumed to be already in cryst coord
     !
     rdotk = twopi * dot_product ( xxq, float(irvec(:, ir)) )
     cfac = exp( ci*rdotk ) / float( ndegen(ir) )
     !
     do ibnd = 1, nbnd
      do jbnd = 1, nbnd
       do ire = 1, nrr_k
        !
        eptmp (ibnd,jbnd,ire,:) = eptmp (ibnd,jbnd,ire,:) + cfac * epmatw ( ibnd, jbnd, ire, :, ir)
        !
       enddo
      enddo
     enddo
     !
  enddo
  !
  !----------------------------------------------------------
  !  STEP 4: un-rotate to Bloch space, fine grid
  !----------------------------------------------------------
  !
  ! epmatf(j) = sum_i eptmp(i) * uf(i,j)
  !
  do ibnd = 1, nbnd
   do jbnd = 1, nbnd
    do ire = 1, nrr_k
     !
     call zgemv ('t', nmodes, nmodes, cone, cuf, nmodes, eptmp(ibnd,jbnd,ire,:), &
        1, czero, epmatf(ibnd,jbnd,ire,:), 1 )
     !
    enddo
   enddo
  enddo
  !
  return
  end subroutine ephwan2blochp
  !