!
    !-----------------------------------------------------------------------
    subroutine loadumat ( nbnd, nbndsub, nks, nkstot, xxq, cu )
    !-----------------------------------------------------------------------
    !
    !   wannier interpolation of e-p vertex:
    !   load rotation matrix on coarse mesh and distribute
    !
    !   Felciiano Giustino, UCB
    !
    !-----------------------------------------------------------------------
    !
#include "f_defs.h"
    USE kinds, only : DP
    USE io_global, ONLY : stdout
    use klist, only : kmap
    use el_phon, only : nkbl
    use control_ph, only : lgamma
    use phcom, only : filukk
    use control_flags, ONLY : iverbosity
#ifdef __PARA
    use para
    USE io_global, ONLY : ionode_id
    USE mp_global, ONLY : mpime
#endif
    implicit none
    !
    !  input variables
    ! 
    integer :: nbnd, nbndsub, nks, nkstot, iunukk
    ! number of bands
    ! number of bands in the optimal subspace
    ! number of kpoints 
    ! total number of kpoints across pools
    real(kind=DP) :: xxq(3)
    ! the qpoint for folding of U
    !
    !  output variables
    !   
    complex(kind=DP) :: cu ( nbnd, nbndsub, nks)
    !
    integer, parameter :: iunukk  = 71               
    ! unit with rotation matrix U(k) from wannier code
    !
    ! work variables 
    !
    integer :: ik, ikk, ikq, ibnd, jbnd, nsize, ik0, ind, nbndfirst
    real(kind=DP) :: aux( 2*nbnd*nbndsub, nkstot) 
    complex(kind=DP) :: cu_big ( nbnd, nbndsub, nkstot)

    !
#ifdef __PARA
    if (mpime.eq.ionode_id) then
#endif
      !
      ! first proc read rotation matrix (coarse mesh) from file
      !
      open ( unit = iunukk, file = filukk, status = 'old', form = 'formatted')
      !
      do ik = 1, nkbl
        !
        if (lgamma) then
          ikk = ik 
          ikq = ik 
        else
          ikk = 2 * ik - 1
          ikq = ikk + 1
        endif
        !
        do ibnd = 1, nbnd
!@ TaC: nbndfirst = 2, not tested
!@        if (ibnd.lt.nbndfirst) cu_big (ibnd, :, ikk) = czero
!@
          do jbnd = 1, nbndsub
            read(iunukk, *) cu_big (ibnd, jbnd, ikk)
          enddo
        enddo
      enddo
      !
      close ( iunukk )
      !
      !  generate U(k+q) through the map 
      !
      if (.not.lgamma) then
         !
         ! here we create the map k+q --> k
         ! (first proc has a copy of all kpoints)
         !
         call createkmap2 ( xxq )
         !
         !  and we generate the matrix for the q-displaced mesh
         !
         do ik = 1, nkbl
           ikq = 2*ik
           ikk = 2*kmap(ik)-1
           cu_big (:, :, ikq) = cu_big (:, :, ikk )  
         enddo
         !
      endif
      !
#ifdef __PARA
    endif
    !
    !  scatter the rotation matrix acroos the pools
    !
    ! first pack the array...
    ! 
    nsize = 2 * nbnd * nbndsub
    do ik = 1, nkstot
      ind = 0
      do ibnd = 1, nbnd
        do jbnd = 1, nbndsub
          ind = ind + 1
          aux( ind, ik) =  real ( cu_big (ibnd, jbnd, ik) ) 
          ind = ind + 1
          aux( ind, ik) = aimag ( cu_big (ibnd, jbnd, ik) ) 
          if (ind.gt.nsize) call errore ('ephwann','problem with packing of cu',1)
        enddo
      enddo
    enddo
    !
    ! ... then scatter
    !
    call poolscatter (nsize, nkstot, aux, nks, aux )
    !
    ! ... finally unpack
    !
    do ik = 1, nks
      ind = 0
      do ibnd = 1, nbnd
        do jbnd = 1, nbndsub
          ind = ind + 1
          cu_big (ibnd, jbnd, ik) = cmplx ( aux( 2*ind-1, ik), aux( 2*ind, ik) )
        enddo
      enddo
    enddo
    !
#endif
    !
    cu = cu_big (:, :, 1:nks)
    !
    return
    end subroutine loadumat
    !----------------------------------------------------------------------
    !