!
  !--------------------------------------------------------------------------
  subroutine kmap_sym ( s, invs, nsym)
  !--------------------------------------------------------------------------
  !
  ! here we generate the map k -> S^-1(k) for all k points and all symmetries
  !
  !  the first proc keeps a copy of all kpoints 
  !
#include "f_defs.h"
  !
  USE kinds, only : DP
  USE io_global, ONLY : stdout
  use pwcom, only : nkstot, xk, wk, at, nk1, nk2, nk3
  use phcom, only : lgamma
  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 :: ix, ik, ikk, ikq, ik0, ri, rj, rk, i, j, k, ism1
  real(kind=DP), parameter :: eps = 1.0e-4
  real(kind=DP) :: ak(3),  xx, yy, zz
  logical :: in_the_list
  integer :: isym, nksqtot, kmapsym (nkstot, 48), s(3,3,48), nsym, invs(48)

  !
  if (lgamma) then
     kunit = 1
  else
     kunit = 2
  endif
  nksqtot = nkstot / kunit
  !
  do ik = 1, nksqtot
    !
    if (lgamma) then
      ikk = ik
      ikq = ik
    else
      ikk = 2 * ik - 1
      ikq = 2 * ik
    endif
    write ( stdout, '(8x,"k(",i4,") = (",3f12.7,"), wk =",f12.7)') &
         ik, (xk (ix, ikk) , ix = 1, 3) , wk (ikk)
    !
    ! go to crystal coordinates and check that the k's are actually 
    ! on a uniform mesh centered at gamma
    !
    ak = xk(:,ikk) 
    call cryst_to_cart ( 1, ak, at, -1)
    !
    !  check that the k's are actually on a uniform mesh centered at gamma
    !
    xx = ak(1)*nk1
    yy = ak(2)*nk2
    zz = ak(3)*nk3
    in_the_list = abs(xx-nint(xx)).le.eps .and. &
       abs(yy-nint(yy)).le.eps .and. abs(zz-nint(zz)).le.eps
    if (.not.in_the_list) &
       call errore('elphon_shuffle_wrap','is this a uniform k-mesh?',1)
    !
    do isym = 1, nsym
      !
      ism1 = invs(isym)
      !
      ! go to integer crys coord and rotate
      !
      i = nint(xx)
      j = nint(yy)
      k = nint(zz)
      !
      ri = s (1, 1, ism1) * i + s (1, 2, ism1) * j + s (1, 3, ism1) * k
      rj = s (2, 1, ism1) * i + s (2, 2, ism1) * j + s (2, 3, ism1) * k
      rk = s (3, 1, ism1) * i + s (3, 2, ism1) * j + s (3, 3, ism1) * k
      ! 
      !  find the index of this S^-1(k) in the k-grid
      !
      ri = mod ( ri + 2*nk1, nk1 )
      rj = mod ( rj + 2*nk2, nk2 )
      rk = mod ( rk + 2*nk3, nk3 )
      !
      kmapsym ( ik, isym ) = ri*nk2*nk3 + rj*nk3 + rk + 1
      !
      write(6,'(3i6)') ik, isym, kmapsym(ik,isym)
      !
    enddo
    !
  enddo
  !
  return
  end subroutine kmap_sym
  !