********************************************************************************
*                                                                              * 
*  2DHF version 1-2003                                                         *
*  Copyright (C) 1996  Jacek Kobus, Leif Laaksonen, Dage Sundholm              *
*                                                                              *
*  This software may be used and distributed according to the terms            *
*  of the GNU General Public License, see README and COPYING.                  *
*                                                                              *
********************************************************************************
c ### prepload ###
c
c     Prepares a share of exchange potentials for every available node 

      subroutine prepload(iorb)
      implicit integer*4 (i-n)
      implicit real*8 (a-h,o-z)
c_mpi
      include 'mpif.h'
      include 'commons_mpi.inc'
c_mpi
      include 'commons8.inc'
      dimension loadpernode(0:maxnumnode-1)
      parameter(iprint=0)

      call MPI_COMM_RANK( MPI_COMM_WORLD, myid, IERROR )
      call MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, IERROR )

      do n=0,maxnumnode-1
         do iorb1=1,norb
            load(iorb,n,iorb1)=0
         enddo
      enddo 

      if (numprocs.eq.1) return

c     There are Numprocs processors 0..Numprocs-1. Processor 0 deals
c     with the Coulomb potential for the orbital iorb. Exchange
c     potentials to be calculated for the specific iorb will be dealt
c     with by nodes 0..Numprocs-1. The work load should be devided
c     evenly between these processors provided that the pair of exchange
c     potentials corresponding to a pair of non-sigma orbitals are
c     calculated on the same node.

      npernode=nexchorb(iorb)/numprocs
      if (mod(nexchorb(iorb),numprocs).ne.0) npernode=npernode+1

      if (myid.eq.0.and.iprint.ne.0) then
         print *,'=========== npernode, numprocs ==========',
     &        npernode,numprocs
      endif

      node=0
      n=0
      kmax=0
      do k=1,nexchorb(iorb)
         if (iexchorb(iorb,k).gt.0) then
            n=n+1
            if (n.le.npernode) then 
               load(iorb,node,n)=iexchorb(iorb,k)
               kmax=k
            else
               n=1
               node=node+1
               if (node.le.numprocs-1) then
                  load(iorb,node,n)=iexchorb(iorb,k)
                  kmax=k
               else
                  goto 500
               endif
            endif
         elseif (iexchorb(iorb,k).lt.0) then
            n=n+1
c           npernode+1 to have a pair of non-sigma orbitals 
c           dealt by the same node
            if (n.le.npernode+1) then
               load(iorb,node,n)=iabs(iexchorb(iorb,k))
               kmax=k
            else
               n=1
               node=node+1
               if (node.le.numprocs-1) then
                  load(iorb,node,n)=iabs(iexchorb(iorb,k))
                  kmax=k
               else
                   goto 500
               endif
            endif
         endif

      if (myid.eq.0.and.iprint.ne.0) then
            print *,'k,iorb,iexchorb(iorb,k)',k,iorb,iexchorb(iorb,k)
            print *,'node,n,load(iorb,node,n)',node,n,load(iorb,node,n)
         endif
      enddo

 500  continue

c      n=npernode+1
      n0=n

      if (myid.eq.0.and.iprint.ne.0) then
         print *,'=========== npernode, kmax ==========',npernode,kmax
      endif

      do k=kmax+1,nexchorb(iorb)
         if (iexchorb(iorb,k).gt.0) then
            n=n+1
            if (n.le.npernode) then 
               load(iorb,node,n)=iexchorb(iorb,k)
               kmax=k
            else
               n=n0+1
               node=node+1
               if (node.le.numprocs-1) then
                  load(iorb,node,n)=iexchorb(iorb,k)
                  kmax=k
               endif
            endif
         elseif (iexchorb(iorb,k).lt.0) then
            n=n+1
            if (n.le.npernode+1) then
               load(iorb,node,n)=iabs(iexchorb(iorb,k))
               kmax=k
            else
               n=1
               node=node+1
               if (node.le.numprocs-1) then
                  load(iorb,node,n)=iabs(iexchorb(iorb,k))
                  kmax=k
               endif
            endif
         endif
         if (myid.eq.0.and.iprint.ne.0) then
            print *,'k,iorb,iexchorb(iorb,k)',k,iorb,iexchorb(iorb,k)
            print *,'node,n,load(iorb,node,n)',node,n,load(iorb,node,n)
         endif
      enddo

c     calculate the load per node, i.e. the number of exchange potentials 
c     being relaxed by a single node

      do i=0,numprocs-1
         loadpernode(i)=0
         do k=1,nexchorb(iorb)
            if (load(iorb,i,k).ne.0) 
     &           loadpernode(i)=loadpernode(i)+1
         enddo
      enddo
      
c     interchange rows 0 and numprocs-1 of load to have 0 node least loaded
      do k=1,nexchorb(iorb)
         li=load(iorb,0,k)
         load(iorb,0,k)=load(iorb,Numprocs-1,k)
         load(iorb,Numprocs-1,k)=li
      enddo
      li=loadpernode(0)
      loadpernode(0)=loadpernode(Numprocs-1)
      loadpernode(Numprocs-1)=li

      if (myid.eq.0.and.iprint.ne.0) then
         write(*,*)
         write(*,*)
         if (iorb.eq.1) write(*,1100) numprocs
         write(*,*)
         write(*,1110) iorb,nexchorb(iorb)
         write(*,1120)
         do i=0,Numprocs-1
            write(*,1122) loadpernode(i),i
         enddo
         write(*,*)
         
         do i=0,Numprocs-1
            write(*,1130) i
            do k=1,nexchorb(iorb)
               if (load(iorb,i,k).ne.0) write(*,1132) load(iorb,i,k)
            enddo
            write(*,*)
         enddo
         
         write(*,*) 'iexchorb(iorb,k)',
     &        (iexchorb(iorb,k),k=1,nexchorb(iorb))
         
         write(*,*) 'iexchadd ',(iexchadd(iorb,iexchorb(iorb,k)),
     &        k=1,nexchorb(iorb))
         write(*,*) 'iexchlen ',(iexchlen(iorb,iexchorb(iorb,k)),
     &        k=1,nexchorb(iorb))
      endif

c      stop 'prepload'

 1000 format(10i7)
 1100 format(5x,'number of processes:',i3)
 1110 format(5x,'load statistics for orbital',i3/
     &     15x,'nb of exchange potentials:  ',i4)
 1120 format(15x,'load per node:    ',$)
 1122 format(i6,'/',i2,$)
 1130 format(15x,'orbitals for node ',i2,': ',$)
 1132 format(i4,$)

      return
      end