!
  !--------------------------------------------------------
  subroutine ktokpmq ( xk, xq, sign, ipool, nkq, nkq_abs)
  !--------------------------------------------------------
  !
  !   For a given k point in cart coord, find the index 
  !   of the corresponding (k + sign*q) point
  !
  !   In the parallel case, determine also the pool number
  !   nkq is the in-pool index, nkq_abs is the absolute
  !   index
  !
  !--------------------------------------------------------
  !
  USE kinds, only : DP
  use io_global, only : stdout
  USE control_flags, ONLY : iverbosity
  use pwcom, only : nk1, nk2, nk3, nkstot, at
#ifdef __PARA
  use para
  USE mp_global, ONLY : nproc, my_pool_id, nproc_pool,    &
     inter_pool_comm, me_pool, root_pool, mpime
  USE mp, ONLY : mp_barrier, mp_bcast
#endif
  implicit none
  !
  real(kind=DP) :: xk (3), xq (3), xkq(3)
  ! input: coordinates of k points and q points
  ! output: coordinates of k+q point
  integer :: sign, ipool, nkq, nkq_abs
  ! input: +1 for searching k+q, -1 for k-q
  ! output: in the parallel case, the pool hosting the k+-q point    
  ! output: the index of k+sign*q
  ! output: the absolute index of k+sign*q (in the full k grid)
  !
  ! work variables
  !
  real(kind=DP) :: xxk (3), xxq (3)
  integer :: nkl, nkbl, nkr, iks, ik, i, j, k, n, jpool
  real(kind=DP) :: xx, yy, zz, eps
  logical :: in_the_list
  !
  ! loosy tolerance, no problem since we use integer comparisons
  eps = 1.d-5
  if (abs(sign).ne.1) call errore('ktokpmq','sign must be +1 or -1',1)
  !
  ! bring k and q in crystal coordinates
  !
  xxk = xk
  xxq = xq
  call cryst_to_cart (1, xxk, at, -1)
  call cryst_to_cart (1, xxq, at, -1)
  !
  !  check that k is actually on a uniform mesh centered at gamma
  !
  xx = xxk(1)*nk1
  yy = xxk(2)*nk2
  zz = xxk(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('ktokpmq','is this a uniform k-mesh?',1)
  !
  !  now add the phonon wavevector and check that k+q falls again on the k grid
  !
  xxk = xxk + float(sign) * xxq
  !
  xx = xxk(1)*nk1
  yy = xxk(2)*nk2
  zz = xxk(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('ktokpmq','k+q does not fall on k-grid',1)
  !
  !  find the index of this k+q in the k-grid
  !
  i = mod ( nint ( xx + 2*nk1), nk1 )
  j = mod ( nint ( yy + 2*nk2), nk2 )
  k = mod ( nint ( zz + 2*nk3), nk3 )
  n = i*nk2*nk3 + j*nk3 + k + 1
  !
  !  take into account the even/odd ordering of the k/k+q points
  n = 2 * n - 1
  !
  !  Now n represents the index of k+sign*q in the original k grid.
  !  In the parallel case we have to find the corresponding pool 
  !  and index in the pool
  !
#ifdef __PARA
  !
  npool = nproc/nproc_pool
  nkbl  = nkstot / kunit
  !
  do jpool = 0, npool-1
    !
    nkl   = kunit * ( nkbl / npool )
    nkr   = ( nkstot - nkl * npool ) / kunit
    !
    !  the reminder goes to the first nkr pools (0...nkr-1)
    !
    if ( jpool < nkr ) nkl = nkl + kunit
    !
    !  the index of the first k point in this pool
    !
    iks = nkl * jpool + 1
    if ( jpool >= nkr ) iks = iks + nkr * kunit
    !
    if (n.ge.iks) then
      ipool = jpool+1
      nkq = n - iks + 1
    endif
    !
  enddo
  !
#else
  !
  nkq = n
  !
#endif
  !
  nkq_abs = n
  !
  end subroutine ktokpmq
  !--------------------------------------------------------