!-----------------------------------------------------------------
      subroutine totrho(vr)
!-----------------------------------------------------------------
!     routine that  calculates the average of the charge density
!     rhor(ir) in the xy plane. The call from the main program
!     should appear just before the beginning of  verlet, after the
!     RHOOFR   call.  Exiting  RHOOFR,  RHOR  contains  charge
!     density in r space.
!     Works only in parallel mode. FG 27/01/03

!     I'VE JUST COPIED EVERYTHING FROM TOTPOT.F90
!     THEN ONE SHOULD READ DENSITY INSTEAD OF POTENTIAL
!     IN THE COMMENTS

      use para_mod
      use elct
      use parm
      implicit none
      include 'mpif.h'
      complex(kind=8) ci
      real(kind=8) vr(nnr,nspin),vave(nr3x,nspin)
      integer istart,istop,i1,i2,i3,is
      real(kind=8) vca(nr1x,nr2x,nr3x,nspin)
      integer root, comm, proc, ierr, displs(nproc), recvcount(nproc)
      real*8, allocatable:: vdist(:,:)
      allocate(vdist(nr1x*nr2x*nr3x,nspin))

      ci=(0.,1.)
      root=0
      comm=mpi_comm_world

      if (me.eq.1) then
         write(6,*) '--------------------------'
         write(6,*) 'entering subroutine TOTRHO'
         write(6,*) '--------------------------'
         write(6,*)
      end if

      do proc=1,nproc
        recvcount(proc) =  ncplane*npp(proc)
        if (proc.eq.1) then
           displs(proc)=0
        else
           displs(proc)=displs(proc-1) + recvcount(proc-1)
        end if
      end do

!     gather the total potential on the first node

      do is=1,nspin
        call mpi_barrier (comm,ierr)
        call mpi_gatherv (vr(1,is),recvcount(me),MPI_REAL8,        &
         &                vdist(:,is),recvcount,displs,MPI_REAL8,  &
         &                root,comm,ierr)
        if (ierr.ne.0) call error('mpi_gatherv','ierr<>0',ierr)
      end do

!------- RESHAPE VR INTO A CARTESIAN MATRIX VCA -----------------------------

         if (me.eq.1) then
           do is=1,nspin
             do i2=1,nr2x
             do i3=1,nr3x
               istart=1+(i2-1)*nr1x+(i3-1)*nr1x*nr2x
               istop =      i2*nr1x+(i3-1)*nr1x*nr2x
               vca(:,i2,i3,is)=vdist(istart:istop,is)
             end do
             end do
           end do
         endif

!------- AVERAGE CARTESIAN POTENTIAL ----------------------------------------------

!     do is=1,nspin
!       do i3=1,nr3x
!          vave(i3,is)=sum(vca(:,:,i3,is))/real(nr1x*nr2x)
!          write(103,'(e15.8)') vave(i3,is)
!       end do
!     end do

      if (me.eq.1) then

      do is=1,nspin
       do i1=1,nr1x
write(6,*) i1,vca(i1,1,1,1)
        do i2=1,nr2x
         do i3=1,nr3x
            write(103,'(e20.10)') vca(i1,i2,i3,is)
         end do
        end do
       end do
      end do

      endif


      if (me.eq.1) then
         write(6,*) '--------------------------'
         write(6,*) 'exiting subroutine  TOTRHO'
         write(6,*) '--------------------------'
      end if

      call mpi_barrier(comm,ierr)

      end