!
  !-----------------------------------------------------------------------
  subroutine loadkmesh
  !-----------------------------------------------------------------------
  !
  !  load fine k mesh and distribute among pools
  !
  !
#include "f_defs.h"
  !
  USE kinds, only : DP
  USE io_global, ONLY : stdout
  use pwcom
  use phcom
  use el_phon
  USE para, ONLY : kunit
#ifdef __PARA
  use para
  USE io_global, ONLY : ionode_id
  USE mp_global, ONLY : nproc, my_pool_id, nproc_pool,      &
     intra_image_comm, inter_pool_comm, me_pool, root_pool, &
     mpime, intra_pool_comm
  USE mp, ONLY : mp_barrier, mp_bcast
#endif
  implicit none
  !
  integer :: ik, ikk, ikq, ibnd, jbnd, imode0, ipol, nsize, rest, &
             ik0, ind, ix, jx, ir, mbnd, irkk, irkq, irk, pbnd, i, iq
  real(kind=DP), parameter :: eps = 1.0e-4
  !
  integer, parameter :: &
    iunqf   = 77,     &! unit with fine q point mesh (cartesian coord)
    iunkkf  = 70,     &! unit with fine k point mesh (crys coord)
    iunukk  = 71       ! unit with rotation matrix U(k) from wannier code
  !
  ! work variables 
  !
  real(kind=DP) :: tmp, xxq(3)
  real(kind=DP), allocatable :: xkf_ (:,:), wkf_ (:), aux(:,:)
  !

  !
  ! read fine qmesh - no pool division at the moment
  !
  open ( unit = iunqf, file = filqf, status = 'old', form = 'formatted')
  read (iunqf,*) nxqf
  allocate ( xqf (3, nxqf), wqf(nxqf) )
  do iq = 1, nxqf
    read (iunqf, *) xqf (:, iq), wqf(iq)
  enddo
  !
  ! bring xqf in crystal coordinates
  !
  call cryst_to_cart (nxqf, xqf, at, -1)
  !
!  if (nqc.ne.nq1*nq2*nq3) call errore &
!    ('ephwann_shuffle','nqc .ne. nq1*nq2*nq3', nqc)
  !
  close ( unit = iunqf )
  write( stdout, '(/5x,"Size of q point mesh for interpolation: ",i10)' ) nxqf
  !


  !
#ifdef __PARA
  if (mpime.eq.ionode_id) then
#endif
    !
    ! read fine mesh
    !
    open ( unit = iunkkf, file = filkf, status = 'old', form = 'formatted')
    read(iunkkf, *) nkstotf 
    write( stdout, '(5x,"Size of k point mesh for interpolation: ",i10)' ) nkstotf 
    !
#ifdef __PARA
  endif
  !
  call mp_bcast (nkstotf, ionode_id, inter_pool_comm)
  call mp_bcast (nkstotf, root_pool, intra_pool_comm)
  !
#endif
  !
  !  temporary storage of k points and weights
  !
  if (lgamma) then
    allocate( xkf_ (3, nkstotf), wkf_ (nkstotf) )
  else
  allocate( xkf_ (3, 2*nkstotf), wkf_ (2*nkstotf) )
  endif
  !
#ifdef __PARA
  if (mpime.eq.ionode_id) then
#endif
    if (lgamma) then
      !
      do ik = 1, nkstotf
        read (iunkkf, *) (xkf_ (ipol, ik), ipol=1,3), wkf_ (ik)
      enddo
      !
      ! bring the points to crystal coordinates
      !
      call cryst_to_cart ( nkstotf, xkf_ ,at,-1)
      !
    else
      !
      ! bring q to crystal coordinates
      !
      xxq = xq
      call cryst_to_cart (1,xxq,at,-1)
      !
      ! read k and generate k+q in crystal coordinates
      !
      do ik = 1, nkstotf
        !
        ikk = 2 * ik - 1
        ikq = ikk + 1
        !
        read (iunkkf, *) (xkf_ (ipol, ikk ), ipol=1,3), wkf_ (ikk)
        !
        !  bring the k point to crystal coordinates
        call cryst_to_cart ( 1, xkf_ (:,ikk), at, -1)
        !
        xkf_ (:, ikq) = xkf_ (:, ikk) + xxq(:)
        wkf_ ( ikq ) = 0.d0
        !
      enddo
      !
      ! redefine nkstotf to include the k+q points
      !
      nkstotf = 2 * nkstotf
      !
    endif
    close ( iunkkf )
    !
    !  check the sum of the kpoint weights: 
    !  must be 2 for spin-unpolarized calculations
    !
    if (abs(sum( wkf_ )-2.d0).gt.eps) call errore &
         ('ephwann','wrong weights in kpoint grid',1)  
    !
#ifdef __PARA
  endif
  !
  call mp_bcast (nkstotf, ionode_id, inter_pool_comm)
  call mp_bcast (nkstotf, root_pool, intra_pool_comm)
  !
  !  scatter the k points of the fine mesh across the pools
  !
  nksf = kunit * ( nkstotf / kunit / npool )
  rest = ( nkstotf - nksf * npool ) / kunit
  if ( mypool .le. rest ) nksf = nksf + kunit
  nksqf = nksf / kunit 
  !
  nsize = 3
  call poolscatter (nsize, nkstotf, xkf_ , nksf, xkf_ )
  nsize = 1
  call poolscatter (nsize, nkstotf, wkf_ , nksf, wkf_ )
  !
  write( stdout, '(5x,"Max number of k points per pool:",7x,i10)' ) nksf 
  !
  nkbl  = nkstot / kunit
#else
  if (lgamma) then
    kunit = 1
  else
    kunit = 2
  endif
  nkbl  = nkstot / kunit
  nksqf = nkbl
#endif
  !
  !  now copy the kpoints and weights into smaller arrays
  !  and deallocate the bigger ones
  !
  allocate ( xkf ( 3, nksf ), wkf( nksf ) )
  xkf = xkf_ ( :, 1:nksf )
  wkf = wkf_ ( 1:nksf )
  deallocate ( xkf_ , wkf_ )
  !
  return
  end subroutine loadkmesh
  !------------------------------------------------------------