c
c  ** pbegin : initialise parallel processing
c
      subroutine pbeginf(nWorkgroups)
      implicit none
      include 'mpif.h'

      Integer MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Common /tcg2mpi_data/ MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Save   /tcg2mpi_data/

      integer nWorkgroups, nProcs, workgroupSize
      integer myWorldRank, myWorkgroup, myWorkgroupRank
c
c  this version assumes a process has NOT been set aside for
c  global index serving
c
      integer ierr1, ierr2, ierr3, ierr4, ierr, n, nw, i, k, iw, ifrom
      integer ic
      Integer me

      logical debug
      data debug/.false./

      logical tcg2mpi_test_verb

      call tcg2mpi_push_verb(1)

      call MPI_INIT(ierr)
      if(ierr.ne.0)call tcg2mpi_errmsg('init',ierr)

c      call mpi_comm_size( MPI_COMM_WORLD, n , ierr1)
c      call mpi_comm_rank( MPI_COMM_WORLD, me, ierr2)
c
c      If( ierr1 .NE. 0 .OR .ierr2 .NE. 0 ) Then
c         call tcg2mpi_errmsg( 'bad initialisation', -1 )
c      End If

c     Divide processors up according to the number of workgroups requested
      call MPI_COMM_SIZE(MPI_COMM_WORLD, nProcs, ierr)

      if (mod(nProcs, nWorkgroups) .ne. 0) call tcg2mpi_errmsg
     &     ('no. procs not divisible by no. workgroups', ierr)
      workgroupSize = nProcs / nWorkgroups

c     Use integer division to assign the first group of procs to
c     workgroup 0, the second to workgroup 1, and so on up to nWorkgroups-1.
      call MPI_COMM_RANK(MPI_COMM_WORLD, myWorldRank, ierr)
      myWorkgroup = myWorldRank / workgroupSize

c     Use mod function to calculate the rank of each processor with respect
c     to its workgroup, starting from 0.
      myWorkgroupRank = mod(myWorldRank, workgroupSize)

c     Create a new communicator for each workgroup
      call MPI_COMM_SPLIT(MPI_COMM_WORLD, myWorkgroup, 
     &     myWorkgroupRank, MPI_COMM_WORKGROUP, ierr)
      if(ierr.ne.0)call tcg2mpi_errmsg('comm_split workgroup',ierr)

c     Create a new communicator for the processors 
c     that share the same rank within their respective workgroups.
      call MPI_COMM_SPLIT(MPI_COMM_WORLD, myWorkgroupRank, 
     &     myWorkgroup, MPI_COMM_COUNTERPARTS, ierr)
      if(ierr.ne.0)call tcg2mpi_errmsg('comm_split ',ierr)

      if (debug) then
         if (myWorldRank .eq. 0) write(6,'(A,I5,3X, A,I5,3X, A,I5)')
     &        "nProcs:", nProcs,
     &        "nWorkgroups:", nWorkgroups,
     &        "workgroupSize:", workgroupSize
         write(6,'(A,I5,3X, A,I5,3X, A,I5)') 
     &        "myWorldRank:", myWorldRank,
     &        "myWorkgroup:", myWorkgroup,
     &        "myWorkgroupRank:", myWorkgroupRank
      endif

      end

c========================================================================
c
c  ** pend
c
      subroutine pend
      implicit none
      integer ierr
      include 'mpif.h'

      Integer MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Common /tcg2mpi_data/ MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Save   /tcg2mpi_data/

      integer icode

c      integer i,n
c      call MPI_COMM_RANK(MPI_COMM_WORLD,i,ierr)
c      call MPI_COMM_SIZE(MPI_COMM_WORLD,n,ierr)

      call MPI_FINALIZE(icode)
      if(icode.ne.0)call tcg2mpi_errmsg('finalize',ierr)

c      write(6,*)'pg_end called  mpi_finalize',i,icode

      call exit(0)
      end
c
c========================================================================
c
c ** worldnodeid() : return index (0 - (nworldnodes-1)) for the
c                   current process in MPI_COMM_WORLD
c
c    Consider whether you should use nodeid() instead before calling this.
c
      integer function worldnodeid()
      implicit none
      include 'mpif.h'

      integer ierr

      call MPI_COMM_RANK(MPI_COMM_WORLD,worldnodeid,ierr)
      if(ierr .ne. 0)call tcg2mpi_errmsg('worldnodeid comm_rank',ierr)

      return
      end
c
c========================================================================
c
c ** nworldnodes() : return number of nodes in MPI_COMM_WORLD
c
c    Consider whether you should use nnodes() instead before calling this.
c
      integer function nworldnodes()
      implicit none
      include 'mpif.h'

      integer ierr

      call MPI_COMM_SIZE(MPI_COMM_WORLD,nworldnodes,ierr)
      if(ierr .ne. 0)call tcg2mpi_errmsg('nworldnodes comm_size',ierr)
      return
      end
c
c========================================================================
c
c ** nodeid() : return index (0 - (nnodes-1)) for the 
c                   current process in its WORKGROUP
c
      integer function nodeid()
      implicit none
      include 'mpif.h'

      Integer MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Common /tcg2mpi_data/ MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Save   /tcg2mpi_data/
      integer ierr

      call MPI_COMM_RANK(MPI_COMM_WORKGROUP,nodeid,ierr)
      if(ierr .ne. 0)call tcg2mpi_errmsg('nodeid comm_rank',ierr)

      return
      end
c
c========================================================================
c
c ** nnodes() : return number of nodes in WORKGROUP
c
      integer function nnodes()
      implicit none
      include 'mpif.h'

      Integer MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Common /tcg2mpi_data/ MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Save   /tcg2mpi_data/
      integer ierr
      call MPI_COMM_SIZE(MPI_COMM_WORKGROUP,nnodes,ierr)
      if(ierr .ne. 0)call tcg2mpi_errmsg('nnodes comm_size',ierr)
      return
      end
c
c========================================================================
c
c ** workgroupid() : return index (0 - (nworkgroups-1)) of the  
c                    workgroup the current process belongs to
c
      integer function workgroupid()
      implicit none
      include 'mpif.h'

      Integer MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Common /tcg2mpi_data/ MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Save   /tcg2mpi_data/
      integer ierr

      call MPI_COMM_RANK(MPI_COMM_COUNTERPARTS,workgroupid,ierr)
      if(ierr .ne. 0)call tcg2mpi_errmsg('workgroupid comm_rank',ierr)

      return
      end
c
c
c========================================================================
c
c ** nworkgroups() : return number of workgroups in the task farm
c
      integer function nworkgroups()
      implicit none
      include 'mpif.h'

      Integer MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Common /tcg2mpi_data/ MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Save   /tcg2mpi_data/
      integer ierr

      call MPI_COMM_SIZE(MPI_COMM_COUNTERPARTS,nworkgroups,ierr)
      if(ierr .ne. 0)call tcg2mpi_errmsg('nworkgroups comm_rank',ierr)

      return
      end
c
c========================================================================
c
c  ** subroutine dgop : double precision workgroup sum
c
c     Double Group OPeration.
c     x(1:n) is a vector present on each process in WORKGROUP. 
c     dgop 'sums' x accross all nodes in WORKGROUP using the commutative 
c     operator op.
c     The result is broadcast to all nodes in WORKGROUP. 
c     Operations to include '+', '*', 'max', 'min', 'absmax', 'absmin'.
c
      subroutine dgop(TYPE, X, N, OP)
      implicit none
      integer TYPE, N
      real*8 X(N)
      character*(*) OP     

      include 'mpif.h'

      Integer MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Common /tcg2mpi_data/ MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Save   /tcg2mpi_data/
      logical tcg2mpi_test_verb

      integer maxbuf, nbatch, ix, ilen, iop, i, ierr
      parameter(maxbuf=10000)
      real*8 buff(maxbuf)

      if(tcg2mpi_test_verb(2))write(6,*)'dgop type=',type,' length=',n

      if(op .eq. '+')then
         iop = MPI_SUM
      else
         write(6,*)op
         call tcg2mpi_errmsg('unsuported op',-1)
      endif
c
      nbatch = ( (n-1) / maxbuf ) + 1
      ix = 1
      do i = 1,nbatch
         ilen = min(maxbuf,1+n-ix)
         call MPI_ALLREDUCE (x(ix),buff,ilen,
     &        MPI_DOUBLE_PRECISION,
     &        iop, MPI_COMM_WORKGROUP, ierr)

         if(ierr .ne. 0)
     &        call tcg2mpi_errmsg('dgop mpi_allreduce',ierr)

         call dcopy( ilen, buff, 1, x(ix), 1 )
         ix = ix + maxbuf
      enddo

      return
      end
c

c========================================================================
c
c  ** subroutine igop : integer workgroup sum
c
c     Integer Group OPeration.
c     x(1:n) is a vector present on each process in WORKGROUP. 
c     igop 'sums' x accross all nodes in WOKRGROUP using the commutative 
c     operator op.
c     The result is broadcast to all nodes in WORKGROUP.
c     Operations to include '+', '*', 'max', 'min', 'absmax', 'absmin'.
c
c
      subroutine igop(TYPE, X, N, OP)
      implicit none
      integer TYPE, N
      integer X(N)
      character*(*) OP     
c
      include 'mpif.h'

      Integer MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Common /tcg2mpi_data/ MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Save   /tcg2mpi_data/

      logical tcg2mpi_test_verb

      integer i, ii, maxbuf, nbatch, ix, iop, ilen, ierr

      parameter(maxbuf=10000)
      integer ibuff(maxbuf)

      if(op .eq. '+')then
         iop = MPI_SUM
      else
         write(6,*)op
         call tcg2mpi_errmsg('unsuported op',-1)
      endif
c
      if(tcg2mpi_test_verb(2))write(6,*)'igop type=',type,' length=',n

      nbatch = ( (n-1) / maxbuf ) + 1
      ix = 1
      do i = 1,nbatch
         ilen = min(maxbuf,1+n-ix)
         call MPI_ALLREDUCE (x(ix),ibuff,ilen,
     &        MPI_INTEGER,
     &        iop, MPI_COMM_WORKGROUP, ierr)

         if(ierr .ne. 0)call 
     &        tcg2mpi_errmsg('integer mpi_allreduce',ierr)

         do ii = 1,ilen
            x(ix+ii-1) = ibuff(ii)
         enddo

         ix = ix + maxbuf

      enddo
      return
      end
c
c========================================================================
c
c ** brdcst : byte-wise broadcast to WORKGROUP
c
      subroutine brdcst(TYPE, BUF, LENBUF, IFROM)
      implicit none
      INTEGER TYPE    
      INTEGER LENBUF  
      INTEGER IFROM   

      include 'mpif.h'

      Integer MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Common /tcg2mpi_data/ MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Save   /tcg2mpi_data/
      integer ierr
      integer buf(*)
      call MPI_BCAST (buf,lenbuf,MPI_BYTE,ifrom,MPI_COMM_WORKGROUP,ierr) 
      if(ierr.ne.0)call tcg2mpi_errmsg('brdcst',ierr)
      return 
      end
c
c ======================================================================
c
c ** snd : send bytewise message to WORKGROUP
c
      subroutine snd(TYPE, BUF, LENBUF, NODE, SYNC)
      implicit none
      INTEGER TYPE      
      INTEGER LENBUF    
      INTEGER NODE      
      INTEGER SYNC      
      integer buf(*)

      integer ierr
      include 'mpif.h'
      logical tcg2mpi_test_verb
      Integer MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Common /tcg2mpi_data/ MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Save   /tcg2mpi_data/

      if(tcg2mpi_test_verb(2))then
         write(6,*)'pg_snd to',node
      endif

      if(sync.ne.1)call tcg2mpi_errmsg('async',0)

      call MPI_SEND (buf,lenbuf,MPI_BYTE, node, type,
     &     MPI_COMM_WORKGROUP,ierr)            

      if (ierr.ne.0) call tcg2mpi_errmsg('send',ierr)
      return
      end
c
c ======================================================================
c
c ** rcv : receive bytewise message from WORKGROUP
c
      SUBROUTINE rcv(TYPE, BUF, LENBUF, LENMES, NODESEL, 
     &     NODEFROM, SYNC)
      implicit none
      INTEGER TYPE      
      INTEGER LENBUF     
c      BYTE BUF(LENBUF)   
      integer BUF(*)   
      INTEGER LENMES     
      INTEGER NODESEL    
      INTEGER NODEFROM   
      INTEGER SYNC       
      
ctfp
      integer ga_nnodes
      integer ga_id_to_msg_id, iii
      integer isel, ifrom
      integer type8, sync8, lenbuf8, lenmes8
      include 'mpif.h'
      integer stat(MPI_STATUS_SIZE)
      integer ierr
      logical tcg2mpi_test_verb

      Integer MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Common /tcg2mpi_data/ MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Save   /tcg2mpi_data/

      if(tcg2mpi_test_verb(2))then
         write(6,*)'pg_rcv from',ifrom
      endif

      if(sync.ne.1)call tcg2mpi_errmsg('async',0)

      ifrom = nodesel
      if (nodesel.eq.-1) ifrom = MPI_ANY_SOURCE

      call mpi_recv  (buf ,lenbuf, MPI_BYTE, ifrom,
     &     type,MPI_COMM_WORKGROUP,stat,ierr)

      if (ierr.ne.0) call tcg2mpi_errmsg('recv',ierr)

***      call MPI_GET_SOURCE(stat,nodefrom)
      nodefrom = stat( MPI_source )

c @@ don't get the true length yet
      lenmes = 0
      return
      end

c     subroutine synch : synchronisation within WORKGROUP
c
      subroutine synch(code)
      implicit none
      integer code
      include 'mpif.h'
      Integer MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Common /tcg2mpi_data/ MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Save   /tcg2mpi_data/
      integer ierr
      logical tcg2mpi_test_verb

      if(tcg2mpi_test_verb(2))write(6,*)'barrier'
      call MPI_BARRIER(MPI_COMM_WORKGROUP,ierr)
      if(ierr.ne.0)call tcg2mpi_errmsg('barrier',ierr)

      return
      end

c     synchworkgroups : synchronise all workgroups together
c                       using counterpart communicator
c
c     Effectively this will synchronise all nodes as the slave nodes
c     (those not in counterparts when called from Tcl) never 
c     operate independently of the masters.
c
      subroutine synchworkgroups(code)
      implicit none
      integer code
      include 'mpif.h'
      Integer MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Common /tcg2mpi_data/ MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Save   /tcg2mpi_data/
      integer ierr
      logical tcg2mpi_test_verb

      if(tcg2mpi_test_verb(2))write(6,*)'counterparts barrier'
      call MPI_BARRIER(MPI_COMM_COUNTERPARTS,ierr)
      if(ierr.ne.0)call tcg2mpi_errmsg('counterparts barrier',ierr)

      return
      end

c
c ======================================================================
c  control of error messages
c
      block data tcg2mpi_verbodat
      implicit none
      integer iverb, ilevel
      common/tcgverbo/iverb(100),ilevel
      data iverb/100*0/
      data ilevel/0/
      end

      subroutine tcg2mpi_push_verb(i)
      implicit none
      integer iverb, ilevel
      common/tcgverbo/iverb(100),ilevel
      integer i
      ilevel = ilevel + 1
      if(ilevel.eq.101)then
         call tcg2mpi_errmsg('recursion gone mad',-1)
      endif
      iverb(ilevel)=i
      return
      end

      subroutine tcg2mpi_pop_verb
      implicit none
      integer iverb, ilevel
      common/tcgverbo/iverb(100),ilevel
      ilevel = ilevel -1
      if(ilevel.le.0)then
         call tcg2mpi_errmsg('pop_verb gone mad',-1)
      endif
      return
      end
c
      logical function tcg2mpi_test_verb(i)
      implicit none
      integer i
      integer iverb, ilevel
      common/tcgverbo/iverb(100),ilevel
      if(ilevel.eq.0)then
         call tcg2mpi_errmsg('bad initialisation or pop_verb error',-1)
      endif
      tcg2mpi_test_verb = i.le.iverb(ilevel)
      return
      end

c
c ======================================================================
c
c ** pg_err : parallel error handling
c
c   code : numerical error code
c
      subroutine parerr(code)
      implicit none
      integer code
      include 'mpif.h'
      integer ierr
      write(6,*) ' PARERR INVOKED, code = ', code
      call mpi_abort(MPI_COMM_WORLD,code,ierr)
      write(6,*) ' RETURN from mpi_abort',ierr
      end

      subroutine tcg2mpi_errmsg(s,i)
      implicit none
      integer nodeid, i, me
      character s*(*)
      write(6,*)'ERROR',s
      me=nodeid()
      write(6,*)'******* fatal error on node',me,' code =',i
      write(6,*)'       ',s
      call parerr(i)
      return
      end

c
c ======================================================================
c ** get world communicator (for, e.g., DL-FIND)
c
      integer function chemsh_comm_world()
      implicit none
      include 'mpif.h'

      chemsh_comm_world = MPI_COMM_WORLD

      return
      end

c
c ======================================================================
c ** get workgroup communicator (for external programs)
c
      integer function chemsh_comm_workgroup()
      implicit none
      include 'mpif.h'

      Integer MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Common /tcg2mpi_data/ MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Save   /tcg2mpi_data/

      chemsh_comm_workgroup = MPI_COMM_WORKGROUP

      return
      end
      
c
c ======================================================================
c ** get workgroup counterparts communicator (for, e.g., DL-FIND)
c

      integer function chemsh_comm_counterparts()
      implicit none
      include 'mpif.h'

      Integer MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Common /tcg2mpi_data/ MPI_COMM_WORKGROUP, MPI_COMM_COUNTERPARTS
      Save   /tcg2mpi_data/

      chemsh_comm_counterparts = MPI_COMM_COUNTERPARTS

      return
      end