!-----------------------------------------------------------------------
  ! Copyright (C) 2010-2015 Henry Lambert, Feliciano Giustino
  ! This file is distributed under the terms of the GNU General Public         
  ! License. See the file `LICENSE' in the root directory of the               
  ! present distribution, or http://www.gnu.org/copyleft.gpl.txt .
  !-----------------------------------------------------------------------
SUBROUTINE davcio_drho( drho, lrec, iunit, nrec, isw )
  !----------------------------------------------------------------------------
  !
  ! ... reads/writes variation of the charge with respect to a perturbation
  ! ... on a file.
  ! ... isw = +1 : gathers data from the processors, writes to a single file
  ! ... isw = -1 : reads data from a single file and distributes them
  !
  USE kinds,     ONLY : DP
  USE fft_base,  ONLY : dfftp, cgather_sym, cscatter_sym
  USE io_global, ONLY : ionode, ionode_id
  USE mp_global, ONLY : inter_pool_comm, me_pool, intra_image_comm
  USE mp,        ONLY : mp_bcast, mp_barrier
  USE gvect,     ONLY : nrx1, nrx2, nrx3, nrxx
  USE noncollin_module,  ONLY : nspin_mag
  !
  IMPLICIT NONE
  !
  INTEGER          :: iunit, lrec, nrec, isw
  COMPLEX(DP) :: drho(nrxx,nspin_mag)
  !
#ifdef __PARA
  !
  ! ... local variables
  !
  INTEGER                        :: is
  LOGICAL   :: exst
  COMPLEX(DP), ALLOCATABLE  :: ddrho(:,:)
  !
  !
  IF ( ionode ) INQUIRE (UNIT = iunit, OPENED = exst)
  CALL mp_bcast(exst,ionode_id, intra_image_comm)
  IF (.NOT.exst) RETURN

  ALLOCATE( ddrho( nrx1 * nrx2 * nrx3 , nspin_mag) )
  !
  IF ( isw == 1 ) THEN
     !
     ! ... First task is the only task allowed to write the file
     !
     DO is = 1, nspin_mag
        !   
        CALL cgather_sym( drho(:,is), ddrho(:,is) )
        !
     END DO
     !
     call mp_barrier(intra_image_comm)
     !
     IF ( ionode ) CALL davcio( ddrho, lrec, iunit, nrec, + 1 )
     !
  ELSE IF ( isw < 0 ) THEN
     !
     ! ... First task reads and broadcasts ddrho to all pools
     !
     IF ( ionode ) CALL davcio( ddrho, lrec, iunit, nrec, - 1 )
     !
     CALL mp_bcast( ddrho, ionode_id, inter_pool_comm )
     !
     ! ... distributes ddrho between between the tasks of the pool
     !
     DO is = 1, nspin_mag
        !   
        CALL cscatter_sym ( ddrho(:,is), drho(:,is) )
        !
     END DO
     !
  END IF
  !
  DEALLOCATE( ddrho )
  !
#else
  !
  CALL davcio( drho, lrec, iunit, nrec, isw )
! !
#endif
  !
  RETURN
  !
END SUBROUTINE davcio_drho