!
!-----------------------------------------------------------------------
subroutine elphel (npe, imode0, dvscfins)
  !-----------------------------------------------------------------------
  !
  !      Calculation of the electron-phonon matrix elements el_ph_mat
  !         <\psi(k+q)|dV_{SCF}/du^q_{i a}|\psi(k)>
  !      Original routine written by Francesco Mauri
  !
#include "f_defs.h"
  use pwcom
  USE wavefunctions_module,  ONLY: evc
  USE kinds, only : DP
  USE io_files, ONLY: iunigk
  use phcom
  use el_phon
  implicit none
  !
  integer :: npe, imode0
  complex(kind=DP) :: dvscfins (nrxxs, nspin, npe)
  ! LOCAL variables
  integer :: ik, ikk, ikq, ipert, mode, nrec, ibnd, jbnd, ir, ig, &
       ios
  complex(kind=DP) , allocatable :: aux1 (:), elphmat (:,:,:)
  complex(kind=DP) :: ZDOTC
  !
  allocate (aux1    ( nrxxs))    
  allocate (elphmat ( nbnd , nbnd , npe))    
  !
  !  Start the loops over the k-points
  !
  if (nksq.gt.1) rewind (unit = iunigk)
  do ik = 1, nksq
     if (nksq.gt.1) then
        read (iunigk, err = 100, iostat = ios) npw, igk
100     call errore ('elphel', 'reading igk', abs (ios) )
     endif
     !
     !  ik = counter of k-points with vector k
     !  ikk= index of k-point with vector k
     !  ikq= index of k-point with vector k+q
     !       k and k+q are alternated if q!=0, are the same if q=0
     !
     if (lgamma) then
        ikk = ik
        ikq = ik
        npwq = npw
     else
        ikk = 2 * ik - 1
        ikq = ikk + 1
     endif
     if (lsda) current_spin = isk (ikk)
     if (.not.lgamma.and.nksq.gt.1) then
        read (iunigk, err = 200, iostat = ios) npwq, igkq
200     call errore ('elphel', 'reading igkq', abs (ios) )
     endif
     !
     call init_us_2 (npwq, igkq, xk (1, ikq), vkb)
     !
     ! read unperturbed wavefuctions psi(k) and psi(k+q)
     !
     if (nksq.gt.1) then
        if (lgamma) then
           call davcio (evc, lrwfc, iuwfc, ikk, - 1)
        else
           call davcio (evc, lrwfc, iuwfc, ikk, - 1)
           call davcio (evq, lrwfc, iuwfc, ikq, - 1)
        endif
     endif
     !
     do ipert = 1, npe
        nrec = (ipert - 1) * nksq + ik
        !
        !  dvbare_q*psi_kpoint is read from file (if available) or recalculated
        !
        if (trans) then
           call davcio (dvpsi, lrbar, iubar, nrec, - 1)
        else
           mode = imode0 + ipert
           ! TODO : .false. or .true. ???
           call dvqpsi_us (ik, mode, u (1, mode), .false. )
        endif
        !
        ! calculate dvscf_q*psi_k
        !
        do ibnd = 1, nbnd
           aux1(:) = (0.d0, 0.d0)
           do ig = 1, npw
              aux1 (nls (igk (ig) ) ) = evc (ig, ibnd)
           enddo
           call cft3s (aux1, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, + 2)
           do ir = 1, nrxxs
              aux1 (ir) = aux1 (ir) * dvscfins (ir, current_spin, ipert)
           enddo
           call cft3s (aux1, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, - 2)
           do ig = 1, npwq
              dvpsi (ig, ibnd) = dvpsi (ig, ibnd) + aux1 (nls (igkq (ig) ) )
           enddo
        end do
        call adddvscf (ipert, ik)

        !
        ! calculate elphmat(j,i)=<psi_{k+q,j}|dvscf_q*psi_{k,i}> for this pertur
        !
        do ibnd =1, nbnd
           do jbnd = 1, nbnd
              elphmat (jbnd, ibnd, ipert) = ZDOTC (npwq, evq (1, jbnd), 1, &
                   dvpsi (1, ibnd), 1)
           enddo
           !
        enddo
     enddo
#ifdef __PARA
     call reduce (2 * nbnd * nbnd * npe, elphmat)
#endif
     !
     !  save all e-ph matrix elements into el_ph_mat
     !
     do ipert = 1, npe
        do jbnd = 1, nbnd
           do ibnd = 1, nbnd
              el_ph_mat (ibnd, jbnd, ik, ipert + imode0) = elphmat (ibnd, jbnd, ipert)
           enddo
        enddo
     enddo
  enddo
  !
  deallocate (elphmat)
  deallocate (aux1)
  !
  return
end subroutine elphel
!