CHARMM Element source/machdep/paral2.src 1.1
C Utility routines for parallel CHARMM on MIMD machines.
C
##IFN PARALLEL
      SUBROUTINE NULL_PL3
      RETURN
##ELSE
      INTEGER FUNCTION MNOD()
C-----------------------------------------------------------------------
C
C Reports node number for this node
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
##IF IBMSP
      INTEGER NUMTASK, TASKID
      CALL MP_ENVIRON(NUMTASK,TASKID)
      MNOD = TASKID
C
##ELIF CM5
      INCLUDE '/usr/include/cm/cmmd_fort.h'
      MNOD = CMMD_SELF_ADDRESS()
C
##ELIF INTEL
      INTEGER MYNODE
      EXTERNAL MYNODE
      MNOD = MYNODE()
C
##ELIF SOCKET PVM PVMC
      INTEGER MNODS
      EXTERNAL MNODS
      MNOD = MNODS()
C
##ELIF T3D T3E
      integer get_my_pe
      mnod = get_my_pe()
C
##ELIF MPI
      INCLUDE 'mpi/mpif.h'
##IF INTEGER8
      INTEGER*4 STATUS,MMNOD,MCOMM
      MCOMM = MPI_COMM_WORLD
      CALL MPI_COMM_RANK(MCOMM,MMNOD,STATUS)
      MNOD=MMNOD
##ELSE
      INTEGER STATUS         
      CALL MPI_COMM_RANK(MPI_COMM_WORLD,MNOD,STATUS)
##ENDIF
C
##ENDIF
C
      RETURN
      END
      INTEGER FUNCTION NNOD()
C-----------------------------------------------------------------------
C
C How many nodes are in the parallel system.
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
##IF IBMSP
      INTEGER NUMTASK, TASKID
      CALL MP_ENVIRON(NUMTASK,TASKID)
      NNOD = NUMTASK
C
##ELIF T3D T3E
      integer GET_NUM_PES
      nnod = GET_NUM_PES()
##ELIF CM5
      INCLUDE '/usr/include/cm/cmmd_fort.h'
      NNOD = CMMD_PARTITION_SIZE()
C
##ELIF INTEL
##IF PARAFULL
      NNOD = NUMNODES()
##ELSE
      INTEGER NUMXPART
      EXTERNAL NUMXPART
      NNOD = NUMXPART()
##ENDIF
C
##ELIF SOCKET PVM PVMC
      INTEGER NNODS
      EXTERNAL NNODS
      NNOD = NNODS()
##ELIF MPI
      INCLUDE 'mpi/mpif.h'
##IF INTEGER8
      INTEGER*4 STATUS,MNNOD,MCOMM
      MCOMM = MPI_COMM_WORLD  
      CALL MPI_COMM_SIZE(MCOMM,MNNOD,STATUS)
      NNOD=MNNOD
##ELSE
      INTEGER STATUS
      CALL MPI_COMM_SIZE(MPI_COMM_WORLD,NNOD,STATUS)
##ENDIF
C
##ENDIF
C
      RETURN
      END
C
      INTEGER FUNCTION LMNOD()
C-----------------------------------------------------------------------
C
C Reports node number for this node. Local routine: IPPMAP ready
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
C
      LMNOD=MYNOD
C
      RETURN
      END
C
      INTEGER FUNCTION LNNOD()
C-----------------------------------------------------------------------
C
C Reports number of nodes. Local routine: IPPMAP ready
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
C
      LNNOD=NUMNOD
C
      RETURN
      END
C
##IF PARASCAL
      INTEGER FUNCTION NUMXPART()
C-----------------------------------------------------------------------
C
C Reports number of nodes on for PARASCAL on hypercube:
C load 0 1 2 3 4 5 6 7 8 9 <program-name> doesn't work!
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
C
      NUMXPART = (NPART+1)*NPART/2
      RETURN
      END
##ENDIF
C
      SUBROUTINE SETIOM()
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
C     SET I/O Mode for same of the machines or open fixed name input
C     file for those which cannot perform REWIND on redirected STDIN.
C     (Intel gamma for example)
C
##IF INTEL
##IFN DELTA PARAGON
C     opening for unit 5
C
C     getcube -tnn > output
C     load charmm23.exe
C     waitcube < charmm.inp # doesn't work (rewind not possible on Gamma)
C     or (delta /rewind works/):
C     mexec "-t(n,m)" -f charmm23.exe < input > output &
C
      OPEN(UNIT=5,FILE='charmm.inp',STATUS='OLD')
##ENDIF
##ENDIF
##IF CM5
      INCLUDE '/usr/include/cm/cmmd_fort.h'
      INTEGER I
C
      I=CMMD_set_open_mode(CMMD_sync_bc)
      I=CMMD_set_io_mode(6,CMMD_independent)
      I=CMMD_set_io_mode(5,CMMD_independent)
      I=CMMD_set_io_mode(0,CMMD_independent)
##ENDIF
C
      RETURN
      END
      SUBROUTINE APSYNC()
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
C     Actual Global Sync routine (sending one byte around)
C
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
C
      INTEGER DUMMY,IDIM,NEXT,ISYNC,IBIT,TYPE,ME
      INTEGER IEOR,IAND
C
C     Has to send at least one byte to really block the program
C
      TYPE=1
      ME=MYNOD
      DO IDIM = 1, NODDIM
         IBIT=2**(IDIM-1)
         NEXT=IEOR(ME,IBIT)
         ISYNC=IAND(ME,IBIT)
         IF(ISYNC.EQ.IBIT) THEN
            CALL GREC(NEXT,TYPE,DUMMY,1)
            CALL GSEN(NEXT,TYPE,DUMMY,1)
         ELSE
            CALL GSEN(NEXT,TYPE,DUMMY,1)
            CALL GREC(NEXT,TYPE,DUMMY,1)
         ENDIF
         TYPE=TYPE+1
      ENDDO
C
      RETURN
      END
C
      SUBROUTINE AGGSYNC()
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
C     Actual General Global SYNC routine (sending one byte around)
C
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
C
      INTEGER ITYPE,DUMMY,ISTEP,ME
      INTEGER ISYNC
C
      ITYPE=1
      ME=MYNOD
      ISYNC=MOD(ME,2)
      DO ISTEP = 1, NUMNOD-1
         IF(ISYNC.EQ.1) THEN
            CALL GRECL(IPPMAP,ITYPE,DUMMY,1)
            CALL GSENR(IPPMAP,ITYPE,DUMMY,1)
         ELSE          
            CALL GSENR(IPPMAP,ITYPE,DUMMY,1)
            CALL GRECL(IPPMAP,ITYPE,DUMMY,1)
         ENDIF
         ITYPE=ITYPE+1
      ENDDO
C
      RETURN
      END
C
      SUBROUTINE GREC(FROM,TYPE,BUF,LEN)
C-----------------------------------------------------------------------
C
C     General RECeive routine. Its purpose is to ease replacement
C     for various platforms (Intel, CM-5, SP1, TCP/IP, ...)
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
      INTEGER FROM, TYPE, LEN
##IFN T3D T3E GNU ULTRA
      BYTE BUF(*)                  !##.not.SGIF90
      integer (kind=1) buf(*)      !## SGIF90
##ELSE
      INTEGER BUF(*)
##ENDIF
C
##IF COMMEASURE
      REAL*8 NSEND,NRECV,LATENCY
      COMMON/NTSTAT/NSEND,NRECV,LATENCY
##ENDIF
C
CCC      integer i
CCC      write(*,'(3i3,30z2)')'from,type,len,buf',
CCC     $     from,type,len,(buf(i),i=1,min(len,30))
C
##IF MPI
      INCLUDE 'mpi/mpif.h'
##IF INTEGER8
      INTEGER*4 M_STAT(MPI_STATUS_SIZE),IERR,REQ,COUNT,
     &          size,MFROM,MTYPE,MLEN,MCOMM,MBYTE
      MFROM=FROM
      MTYPE=TYPE
      MLEN=LEN
      MCOMM = MPI_COMM_WORLD
      MBYTE = MPI_BYTE
      IF(MLEN.GT.0) THEN
         CALL MPI_IRECV(BUF,MLEN,MBYTE,MFROM,MTYPE,
     @                  MCOMM,REQ,IERR)
         IF(IERR.NE.MPI_SUCCESS) 
     @         CALL WRNDIE(-4,'<GREC>','MPI IRECV ERROR IN GREC')
         CALL MPI_WAIT(REQ,M_STAT,IERR)
         IF(IERR.NE.MPI_SUCCESS) 
     @         CALL WRNDIE(-4,'<GREC>','MPI WAIT ERROR IN GREC')
      ENDIF
##ELSE
      INTEGER M_STAT(MPI_STATUS_SIZE),IERR,REQ,COUNT,size
      IF(LEN.GT.0) THEN
         CALL MPI_IRECV(BUF,LEN,MPI_BYTE,FROM,TYPE,MPI_COMM_WORLD,REQ,
     @                  IERR)
         IF(IERR.NE.MPI_SUCCESS) 
     @         CALL WRNDIE(-4,'<GREC>','MPI IRECV ERROR IN GREC')
         CALL MPI_WAIT(REQ,M_STAT,IERR)
         IF(IERR.NE.MPI_SUCCESS) 
     @         CALL WRNDIE(-4,'<GREC>','MPI WAIT ERROR IN GREC')
      ENDIF
##ENDIF
##ELIF IBMSP
      INTEGER N
      CALL MP_BRECV(BUF,LEN,FROM,TYPE,N)
      IF(LEN.NE.N) CALL DIE
C
##ELIF CM5
      INCLUDE '/usr/include/cm/cmmd_fort.h'
      INTEGER ITMP
      ITMP=CMMD_RECEIVE_BLOCK(FROM,TYPE,BUF,LEN)
C
##ELIF T3D T3E
##IF SHMEM
Cc##INCLUDE '~/charmm_fcm/shmem.fcm'
      integer i
      call barrier()
cdir$ suppress
      do i = 1, len
         buf( i ) = sh_recv( i )
      end do
##ELSE
##IF T3D
      INCLUDE "/usr/include/mpp/fpvm3.h"
##ELSE 
      INCLUDE "fpvm3.h"
##ENDIF
      INTEGER IOK,STRIDE
      PARAMETER (STRIDE=1)

      CALL PVMFRECV(-1, TYPE, IOK)
      IF (IOK.LT.0) CALL PVMFPERROR(' GREC recv failed >',iok)
      CALL PVMFUNPACK(BYTE1, BUF, LEN, STRIDE, IOK)
      IF (IOK.LT.0) CALL PVMFPERROR(' GREC upck failed >',iok)
##ENDIF
##ELIF INTEL
      CALL CRECV(TYPE,BUF,LEN)
C
##ELIF SOCKET
c added length control since 0 length send/receives appear not to work
c w/ sockets on wither AMD or with Linux7.1. clbiii
      if(len.gt.0)CALL FRECV(FROM,BUF,LEN)
##ELIF PVMC
##INCLUDE '~/charmm_fcm/parallel.fcm'
      INCLUDE 'pvm/fpvm3.h'
      INTEGER BUFID,INFO
      CALL PVMFRECV(TIDS(FROM),TYPE,BUFID)
      CALL PVMFUNPACK(BYTE1,BUF,LEN,1,INFO)
##ELIF PVM
      CALL FRECV_PVM(FROM,TYPE,BUF,LEN)
C
##ELIF TERRA
      CALL WRNDIE(-5,'<GREC>','Should never be called on TERRA.')
C
##ENDIF
C
##IF COMMEASURE
      NRECV=NRECV+LEN
      LATENCY=LATENCY+1.0D0
##ENDIF
C
      RETURN
      END
C
      SUBROUTINE GRECMAP(FROM,TYPE,BUF,LEN)
C-----------------------------------------------------------------------
C
C     General RECeive routine. Its purpose is to ease replacement
C     for various platforms (Intel, CM-5, SP1, TCP/IP, ...)
C     This one maps the processor ID according to IPPMAP array
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
C
      INTEGER FROM, TYPE, LEN
##IFN T3D T3E GNU ULTRA
      BYTE BUF(*)
##ELSE
      INTEGER BUF(*)
##ENDIF
C
##IF COMMEASURE
      REAL*8 NSEND,NRECV,LATENCY
      COMMON/NTSTAT/NSEND,NRECV,LATENCY
##ENDIF
C
CCC      integer i
CCC      write(*,'(3i3,30z2)')'from,type,len,buf',
CCC     $     from,type,len,(buf(i),i=1,min(len,30))
C
##IF MPI
      INCLUDE 'mpi/mpif.h'
##IF GNUALPHA
      INTEGER*4 M_STAT(MPI_STATUS_SIZE),IERR,REQ,COUNT,
     &          size,MFROM,MTYPE,MLEN
      MFROM=IPPMAP(FROM)
      MTYPE=TYPE
      MLEN=LEN
      IF(MLEN.GT.0) THEN
         CALL MPI_IRECV(BUF,MLEN,MPI_BYTE,MFROM,MTYPE,
     @                  MPI_COMM_WORLD,REQ,IERR)
         IF(IERR.NE.MPI_SUCCESS) 
     @         CALL WRNDIE(-4,'<GREC>','MPI IRECV ERROR IN GREC')
         CALL MPI_WAIT(REQ,M_STAT,IERR)
         IF(IERR.NE.MPI_SUCCESS) 
     @         CALL WRNDIE(-4,'<GREC>','MPI WAIT ERROR IN GREC')
      ENDIF
##ELSE
      INTEGER M_STAT(MPI_STATUS_SIZE),IERR,REQ,COUNT,size
      IF(LEN.GT.0) THEN
         CALL MPI_IRECV(BUF,LEN,MPI_BYTE,IPPMAP(FROM),TYPE,
     $        MPI_COMM_WORLD,REQ,IERR)
         IF(IERR.NE.MPI_SUCCESS) 
     @         CALL WRNDIE(-4,'<GREC>','MPI IRECV ERROR IN GREC')
         CALL MPI_WAIT(REQ,M_STAT,IERR)
         IF(IERR.NE.MPI_SUCCESS) 
     @         CALL WRNDIE(-4,'<GREC>','MPI WAIT ERROR IN GREC')
      ENDIF
##ENDIF
##ELIF IBMSP
      INTEGER N
      CALL MP_BRECV(BUF,LEN,IPPMAP(FROM),TYPE,N)
      IF(LEN.NE.N) CALL DIE
C
##ELIF CM5
      INCLUDE '/usr/include/cm/cmmd_fort.h'
      INTEGER ITMP
      ITMP=CMMD_RECEIVE_BLOCK(IPPMAP(FROM),TYPE,BUF,LEN)
C
##ELIF T3D T3E
##IF SHMEM
Cc##INCLUDE '~/charmm_fcm/shmem.fcm'
      integer i
      call barrier()
cdir$ suppress
      do i = 1, len
         buf( i ) = sh_recv( i )
      end do
##ELSE
##IF T3D
      INCLUDE "/usr/include/mpp/fpvm3.h"
##ELSE 
      INCLUDE "fpvm3.h"
##ENDIF
      INTEGER IOK,STRIDE
      PARAMETER (STRIDE=1)

      CALL PVMFRECV(-1, TYPE, IOK)
      IF (IOK.LT.0) CALL PVMFPERROR(' GREC recv failed >',iok)
      CALL PVMFUNPACK(BYTE1, BUF, LEN, STRIDE, IOK)
      IF (IOK.LT.0) CALL PVMFPERROR(' GREC upck failed >',iok)
##ENDIF
##ELIF INTEL
      CALL CRECV(TYPE,BUF,LEN)
C
##ELIF SOCKET
      CALL FRECV(IPPMAP(FROM),BUF,LEN)
##ELIF PVMC
##INCLUDE '~/charmm_fcm/parallel.fcm'
      INCLUDE 'pvm/fpvm3.h'
      INTEGER BUFID,INFO
      CALL PVMFRECV(TIDS(IPPMAP(FROM)),TYPE,BUFID)
      CALL PVMFUNPACK(BYTE1,BUF,LEN,1,INFO)
##ELIF PVM
      CALL FRECV_PVM(IPPMAP(FROM),TYPE,BUF,LEN)
C
##ELIF TERRA
      CALL WRNDIE(-5,'<GREC>','Should never be called on TERRA.')
C
##ENDIF
C
##IF COMMEASURE
      NRECV=NRECV+LEN
      LATENCY=LATENCY+1.0D0
##ENDIF
C
      RETURN
      END
C
      SUBROUTINE GRECL(IPPMAP,TYPE,BUF,LEN)
C-----------------------------------------------------------------------
C
C     RECeive from Left routine. Used in general communication routines
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
      INTEGER TYPE, LEN, IPPMAP(0:*)
##IFN T3D T3E GNU ULTRA
      BYTE BUF(*)                  !##.not.SGIF90
      integer (kind=1) buf(*)      !## SGIF90
##ELSE
      INTEGER BUF(*)
##ENDIF
##IF TERRA
      INTEGER RC
      CALL TCOMMLIBRECVLEFT(BUF,LEN,RC)
##ELSE
      INTEGER ME,NUMNOD,LMNOD,LNNOD
      ME=LMNOD()
      NUMNOD=LNNOD()
      CALL GREC(IPPMAP(MOD(ME+NUMNOD-1,NUMNOD)),TYPE,BUF,LEN)
##ENDIF
      RETURN
      END
C
      SUBROUTINE GRECR(IPPMAP,TYPE,BUF,LEN)
C-----------------------------------------------------------------------
C
C     RECeive from Right routine. Used in general communication routines
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
      INTEGER TYPE, LEN, IPPMAP(0:*)
##IFN T3D T3E GNU ULTRA
      BYTE BUF(*)                  !##.not.SGIF90
      integer (kind=1) buf(*)      !## SGIF90
##ELSE
      INTEGER BUF(*)
##ENDIF
##IF TERRA
      INTEGER RC
      CALL TCOMMLIBRECVRIGHT(BUF,LEN,RC)
##ELSE
      INTEGER NUMNOD,LNNOD,ME,LMNOD
      ME=LMNOD()
      NUMNOD=LNNOD()
      CALL GREC(IPPMAP(MOD(ME+1,NUMNOD)),TYPE,BUF,LEN)
##ENDIF
      RETURN
      END
C
##IF T3D T3E
      SUBROUTINE GRECC(FROM,TYPE,BUF,LEN)
C-----------------------------------------------------------------------
C
C     T3D RECeive routine for character arrays.
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
      INTEGER FROM, TYPE, LEN
      CHARACTER*(*) BUF(1)
##IF COMMEASURE
      REAL*8 NSEND,NRECV,LATENCY
      COMMON/NSTAT/NSEND,NRECV,LATENCY
##ENDIF
C
##IF SHMEM
Cc##INCLUDE '~/charmm_fcm/shmem.fcm'
      integer i
      call barrier()
cdir$ suppress
      do i = 1, len
         buf( i ) = sh_recv( i )
      end do
##ELSE
##IF T3D
      INCLUDE "/usr/include/mpp/fpvm3.h"
##ELSE 
      INCLUDE "fpvm3.h"
##ENDIF
      INTEGER IOK,STRIDE,I
      PARAMETER (STRIDE=1)

      CALL PVMFRECV(-1, TYPE, IOK)
      IF (IOK.LT.0) CALL PVMFPERROR(' GRECC recv failed >',iok)
      CALL PVMFUNPACK(BYTE1, BUF( 1 ), LEN, STRIDE, IOK)
      IF (IOK.LT.0) CALL PVMFPERROR(' GRECC upck failed >',iok)
##ENDIF
C
##IF COMMEASURE
      NRECV=NRECV+LEN
      LATENCY=LATENCY+1.0D0
##ENDIF
C
      RETURN
      END
##ENDIF
      SUBROUTINE GSEN(TO,TYPE,BUF,LEN)
C-----------------------------------------------------------------------
C
C     General SENd routine. Its purpose is to ease replacement
C     for various platforms (Intel, CM-5, SP1, TCP/IP, ...)
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
      INTEGER TO, TYPE, LEN
##IFN T3D T3E GNU ULTRA
      BYTE BUF(*)                  !##.not.SGIF90
      integer (kind=1) buf(*)      !## SGIF90
##ELSE
      INTEGER BUF(*)
##ENDIF
C
##IF COMMEASURE
      REAL*8 NSEND,NRECV,LATENCY
      COMMON/NTSTAT/NSEND,NRECV,LATENCY
##ENDIF
C
##IF MPI
      INCLUDE 'mpi/mpif.h'
##IF INTEGER8
      INTEGER*4 IERR,REQ,M_STAT(MPI_STATUS_SIZE),
     &          MTO,MTYPE,MLEN,MCOMM,MBYTE
      MTO=TO
      MTYPE=TYPE
      MLEN=LEN
      MCOMM = MPI_COMM_WORLD
      MBYTE = MPI_BYTE
      IF(MLEN.GT.0) THEN
         CALL MPI_ISEND(BUF,MLEN,MBYTE,MTO,MTYPE,
     @                  MCOMM,REQ,IERR)
         IF(IERR.NE.MPI_SUCCESS) 
     @         CALL WRNDIE(-4,'<GSEN>','MPI ISEND ERROR IN GSEN')
         CALL MPI_WAIT(REQ,M_STAT,IERR)
         IF(IERR.NE.MPI_SUCCESS) 
     @         CALL WRNDIE(-4,'<GSEN>','MPI WAIT ERROR IN GSEN')
      ENDIF
##ELSE
      INTEGER IERR,REQ,M_STAT(MPI_STATUS_SIZE)
      IF(LEN.GT.0) THEN
         CALL MPI_ISEND(BUF,LEN,MPI_BYTE,TO,TYPE,MPI_COMM_WORLD,REQ,
     @                  IERR)
         IF(IERR.NE.MPI_SUCCESS) 
     @         CALL WRNDIE(-4,'<GSEN>','MPI ISEND ERROR IN GSEN')
         CALL MPI_WAIT(REQ,M_STAT,IERR)
         IF(IERR.NE.MPI_SUCCESS) 
     @         CALL WRNDIE(-4,'<GSEN>','MPI WAIT ERROR IN GSEN')
      ENDIF
##ENDIF
C
##ELIF IBMSP
      CALL MP_BSEND(BUF,LEN,TO,TYPE)
##ELIF CM5
      INCLUDE '/usr/include/cm/cmmd_fort.h'
      INTEGER ITMP
      ITMP=CMMD_SEND_BLOCK(TO,TYPE,BUF,LEN)
C
##ELIF T3D T3E
##IF SHMEM
Cc##INCLUDE '~/charmm_fcm/shmem.fcm'
      integer i
      do i = 1, len
         sh_send( i ) = buf( i )
      end do
      call shmem_put( sh_recv, sh_send, len, to )
      call barrier()
##ELSE
##IF T3D
      INCLUDE "/usr/include/mpp/fpvm3.h"
##ELSE 
      INCLUDE "fpvm3.h"
##ENDIF
      INTEGER IOK

      CALL PVMFINITSEND( PVMINPLACE, IOK)
      IF (IOK.LT.0) CALL PVMFPERROR(' GSEN init failed >',iok)
      CALL PVMFPACK( BYTE1, BUF, LEN, 1, IOK)
      IF (IOK.LT.0) CALL PVMFPERROR(' GSEN pack failed >',iok)
      CALL PVMFSEND( TO, TYPE, IOK)
      IF (IOK.LT.0) CALL PVMFPERROR(' GSEN send failed >',iok)
##ENDIF
##ELIF TERRA
      CALL WRNDIE(-5,'<GSEN>','Should never be called on TERRA.')
C
##ELIF PVMC
##INCLUDE '~/charmm_fcm/parallel.fcm'
      INCLUDE 'pvm/fpvm3.h'
      INTEGER BUFID,INFO
##IF SGIMP CSPP ALPHAMP
      CALL PVMFINITSEND(PVMRAW,BUFID)
##ELSE
      CALL PVMFINITSEND(PVMDATAINPLACE,BUFID)
##ENDIF
      CALL PVMFPACK(BYTE1,BUF,LEN,1,INFO)
      CALL PVMFSEND(TIDS(TO),TYPE,BUFID)
C
##ELIF PVM
      CALL FSEND_PVM(TO,TYPE,BUF,LEN)
C
##ELIF INTEL
      CALL CSEND(TYPE,BUF,LEN,TO,0)
C
##ELIF SOCKET
c added length control since 0 length send/receives appear not to work
c w/ sockets on wither AMD or with Linux7.1. clbiii
      if(len .gt.0)CALL FSEND(TO,BUF,LEN)
C
##ENDIF
C
##IF COMMEASURE
      NSEND=NSEND+LEN
      LATENCY=LATENCY+1.0D0
##ENDIF
C
      RETURN
      END
C
      SUBROUTINE GSENMAP(TO,TYPE,BUF,LEN)
C-----------------------------------------------------------------------
C
C     General SENd routine. Its purpose is to ease replacement
C     for various platforms (Intel, CM-5, SP1, TCP/IP, ...)
C     This one maps the processor ID according to IPPMAP array
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
C
      INTEGER TO, TYPE, LEN
##IFN T3D T3E GNU ULTRA
      BYTE BUF(*)
##ELSE
      INTEGER BUF(*)
##ENDIF
C
##IF COMMEASURE
      REAL*8 NSEND,NRECV,LATENCY
      COMMON/NTSTAT/NSEND,NRECV,LATENCY
##ENDIF
C
##IF MPI
      INCLUDE 'mpi/mpif.h'
##IF GNUALPHA
      INTEGER*4 IERR,REQ,M_STAT(MPI_STATUS_SIZE),
     &          MTO,MTYPE,MLEN
      MTO=IPPMAP(TO)
      MTYPE=TYPE
      MLEN=LEN
      IF(MLEN.GT.0) THEN
         CALL MPI_ISEND(BUF,MLEN,MPI_BYTE,MTO,MTYPE,
     @                  MPI_COMM_WORLD,REQ,IERR)
         IF(IERR.NE.MPI_SUCCESS) 
     @         CALL WRNDIE(-4,'<GSEN>','MPI ISEND ERROR IN GSEN')
         CALL MPI_WAIT(REQ,M_STAT,IERR)
         IF(IERR.NE.MPI_SUCCESS) 
     @         CALL WRNDIE(-4,'<GSEN>','MPI WAIT ERROR IN GSEN')
      ENDIF
##ELSE
      INTEGER IERR,REQ,M_STAT(MPI_STATUS_SIZE)
      IF(LEN.GT.0) THEN
         CALL MPI_ISEND(BUF,LEN,MPI_BYTE,IPPMAP(TO),TYPE,
     $        MPI_COMM_WORLD,REQ,IERR)
         IF(IERR.NE.MPI_SUCCESS) 
     @         CALL WRNDIE(-4,'<GSEN>','MPI ISEND ERROR IN GSEN')
         CALL MPI_WAIT(REQ,M_STAT,IERR)
         IF(IERR.NE.MPI_SUCCESS) 
     @         CALL WRNDIE(-4,'<GSEN>','MPI WAIT ERROR IN GSEN')
      ENDIF
##ENDIF
C
##ELIF IBMSP
      CALL MP_BSEND(BUF,LEN,IPPMAP(TO),TYPE)
##ELIF CM5
      INCLUDE '/usr/include/cm/cmmd_fort.h'
      INTEGER ITMP
      ITMP=CMMD_SEND_BLOCK(IPPMAP(TO),TYPE,BUF,LEN)
C
##ELIF T3D T3E
##IF SHMEM
Cc##INCLUDE '~/charmm_fcm/shmem.fcm'
      integer i
      do i = 1, len
         sh_send( i ) = buf( i )
      end do
      call shmem_put( sh_recv, sh_send, len, ippmap(to) )
      call barrier()
##ELSE
##IF T3D
      INCLUDE "/usr/include/mpp/fpvm3.h"
##ELSE 
      INCLUDE "fpvm3.h"
##ENDIF
      INTEGER IOK

      CALL PVMFINITSEND( PVMINPLACE, IOK)
      IF (IOK.LT.0) CALL PVMFPERROR(' GSEN init failed >',iok)
      CALL PVMFPACK( BYTE1, BUF, LEN, 1, IOK)
      IF (IOK.LT.0) CALL PVMFPERROR(' GSEN pack failed >',iok)
      CALL PVMFSEND( IPPMAP(TO), TYPE, IOK)
      IF (IOK.LT.0) CALL PVMFPERROR(' GSEN send failed >',iok)
##ENDIF
##ELIF TERRA
      CALL WRNDIE(-5,'<GSEN>','Should never be called on TERRA.')
C
##ELIF PVMC
##INCLUDE '~/charmm_fcm/parallel.fcm'
      INCLUDE 'pvm/fpvm3.h'
      INTEGER BUFID,INFO
##IF SGIMP CSPP ALPHAMP
      CALL PVMFINITSEND(PVMRAW,BUFID)
##ELSE
      CALL PVMFINITSEND(PVMDATAINPLACE,BUFID)
##ENDIF
      CALL PVMFPACK(BYTE1,BUF,LEN,1,INFO)
      CALL PVMFSEND(TIDS(IPPMAP(TO)),TYPE,BUFID)
C
##ELIF PVM
      CALL FSEND_PVM(TO,TYPE,BUF,LEN)
C
##ELIF INTEL
      CALL CSEND(TYPE,BUF,LEN,IPPMAP(TO),0)
C
##ELIF SOCKET
      if(len.gt.0)CALL FSEND(IPPMAP(TO),BUF,LEN)
C
##ENDIF
C
##IF COMMEASURE
      NSEND=NSEND+LEN
      LATENCY=LATENCY+1.0D0
##ENDIF
C
      RETURN
      END
C
      SUBROUTINE GSENL(IPPMAP,TYPE,BUF,LEN)
C-----------------------------------------------------------------------
C
C     SENd to Left routine. Used in general communication routines
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
      INTEGER TYPE, LEN, IPPMAP(0:*)
##IFN T3D T3E GNU ULTRA
      BYTE BUF(*)                  !##.not.SGIF90
      integer (kind=1) buf(*)      !## SGIF90
##ELSE
      INTEGER BUF(*)
##ENDIF
##IF TERRA
      INTEGER RC
      CALL TCOMMLIBSENDLEFT(BUF,LEN,RC)
##ELSE
      INTEGER ME,NUMNOD,LMNOD,LNNOD
      ME=LMNOD()
      NUMNOD=LNNOD()
      CALL GSEN(IPPMAP(MOD(ME+NUMNOD-1,NUMNOD)),TYPE,BUF,LEN)
##ENDIF
      RETURN
      END
C
      SUBROUTINE GSENR(IPPMAP,TYPE,BUF,LEN)
C-----------------------------------------------------------------------
C
C     SENd to Right routine. Used in general communication routines
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
      INTEGER TYPE, LEN, IPPMAP(0:*)
##IFN T3D T3E GNU ULTRA
      BYTE BUF(*)                  !##.not.SGIF90
      integer (kind=1) buf(*)      !## SGIF90
##ELSE
      INTEGER BUF(*)
##ENDIF
##IF TERRA
      INTEGER RC
      CALL TCOMMLIBSENDRIGHT(BUF,LEN,RC)
##ELSE
      INTEGER NUMNOD,LNNOD,ME,LMNOD
      ME=LMNOD()
      NUMNOD=LNNOD()
      CALL GSEN(IPPMAP(MOD(ME+1,NUMNOD)),TYPE,BUF,LEN)
##ENDIF
      RETURN
      END
C
##IF T3D T3E
      SUBROUTINE GSENC(TO,TYPE,BUF,LEN)
C-----------------------------------------------------------------------
C
C     T3D SENd routine. Its purpose is to ease replacement
C     for various platforms (Intel, CM-5, SP1, TCP/IP, ...)
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
      INTEGER TO, TYPE, LEN
      CHARACTER*(*) BUF(*)
##IF COMMEASURE
      REAL*8 NSEND,NRECV
      COMMON/STAT/NSEND,NRECV
##ENDIF
##IF SHMEM
Cc##INCLUDE '~/charmm_fcm/shmem.fcm'
      integer i
      do i = 1, len
         sh_send( i ) = buf( i )
      end do
      call shmem_put( sh_recv, sh_send, len, to )
      call barrier()
##ELSE
##IF T3D
      INCLUDE "/usr/include/mpp/fpvm3.h"
##ELSE 
      INCLUDE "fpvm3.h"
##ENDIF
      INTEGER IOK,I
      CALL PVMFINITSEND( PVMINPLACE, IOK)
      IF (IOK.LT.0) CALL PVMFPERROR(' GSENC init failed >',iok)
      CALL PVMFPACK( BYTE1, BUF( 1 ), LEN, 1, IOK)
      IF (IOK.LT.0) CALL PVMFPERROR(' GSENC pack failed >',iok)
      CALL PVMFSEND( TO, TYPE, IOK)
      IF (IOK.LT.0) CALL PVMFPERROR(' GSENC send failed >',iok)
##ENDIF
##IF COMMEASURE
      NSEND=NSEND+LEN
##ENDIF
C
      RETURN
      END
##ENDIF
C
      SUBROUTINE GRECSEN(TO,TYPE,RBUF,RLEN,SBUF,SLEN)
C-----------------------------------------------------------------------
C
C     General RECeive and SENd routine. Its purpose is to ease replacement
C     for various platforms (Intel, CM-5, SP1, TCP/IP, ...)
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
      INTEGER TO, TYPE, RLEN, SLEN
      REAL*8 RBUF(*),SBUF(*)
C
##IF COMMEASURE
      REAL*8 NSEND,NRECV,LATENCY
      COMMON/NTSTAT/NSEND,NRECV,LATENCY
##ENDIF
C
##IF MPI
      INCLUDE 'mpi/mpif.h'
##IF INTEGER8
      INTEGER*4 REQ(2),STATUS_ARRAY(MPI_STATUS_SIZE,2),IERR,
     &          MTO,MTYPE,MRLEN,MSLEN,WTALL,MDP,MCOMM
##ELSE
      INTEGER REQ(2),STATUS_ARRAY(MPI_STATUS_SIZE,2),IERR,
     &        MTO,MTYPE,MRLEN,MSLEN,WTALL,MDP,MCOMM
##ENDIF
##IF BUFFERED (buffers)
##INCLUDE '~/charmm_fcm/lobos.fcm'
      INTEGER ILENB,JLENB,IPT,JPT,IL,JL,IX    !##.not.INTEGER8
      INTEGER*4 ILENB,JLENB,IPT,JPT,IL,JL,IX  !##INTEGER8
CCC      integer mnod
C
C
C     Some networks perform better with the fixed small buffer size
C     To activate the code put BUFFERED in the pref.dat
C     Works only for CMPI MPI keywords!
C
      MDP = MPI_DOUBLE_PRECISION
      MCOMM = MPI_COMM_WORLD
      MTO=TO
      MTYPE=TYPE
      WTALL=2
      IPT=1
      JPT=1
      ILENB=SLEN
      JLENB=RLEN
 10   IX=0
      IF(JLENB.GT.NICBUF) THEN
         JL=NICBUF
      ELSE
         JL=JLENB
      ENDIF
      IF(ILENB.GT.NICBUF) THEN
         IL=NICBUF
      ELSE
         IL=ILENB
      ENDIF
CCC      write(*,*)'me,to,il,jl=',mnod(),to,il,jl
      IF(IL.GT.0) THEN
         IX=IX+1
         CALL MPI_ISEND(SBUF(IPT),IL,MDP,MTO,MTYPE,
     $        MCOMM,REQ(IX),IERR)
         IF(IERR.NE.MPI_SUCCESS) 
     @        CALL WRNDIE(-4,'<GRECSEN>','MPI ISEND ERROR')
      ENDIF
      IF(JL.GT.0) THEN
         IX=IX+1
         CALL MPI_IRECV(RBUF(JPT),JL,MDP,MTO,MTYPE,
     $        MCOMM,REQ(IX),IERR)
         IF(IERR.NE.MPI_SUCCESS) 
     @        CALL WRNDIE(-4,'<GRECSEN>','MPI IRECV ERROR')
      ENDIF
      IF(IX.GT.0) THEN
         CALL MPI_WAITALL(WTALL,REQ,STATUS_ARRAY,IERR)
         IF(IERR.NE.MPI_SUCCESS) 
     @        CALL WRNDIE(-4,'<GRECSEN>','MPI WAITALL ERROR')
      ENDIF
      IPT=IPT+IL
      JPT=JPT+JL
      ILENB=ILENB-IL
      JLENB=JLENB-JL
      IF((ILENB.GT.0).OR.(JLENB.GT.0)) GOTO 10
##ELSE (buffers)
C
      MTO=TO
      MTYPE=TYPE
      MRLEN=RLEN
      MSLEN=SLEN
      WTALL=2
      MDP = MPI_DOUBLE_PRECISION
      MCOMM = MPI_COMM_WORLD

      CALL MPI_ISEND(SBUF,MSLEN,MDP,MTO,MTYPE,
     $     MCOMM,REQ(1),IERR)
      IF(IERR.NE.MPI_SUCCESS) 
     @     CALL WRNDIE(-4,'<GRECSEN>','MPI ISEND ERROR')
      CALL MPI_IRECV(RBUF,MRLEN,MDP,MTO,MTYPE,
     $     MCOMM,REQ(2),IERR)
      IF(IERR.NE.MPI_SUCCESS) 
     @     CALL WRNDIE(-4,'<GRECSEN>','MPI IRECV ERROR')
      CALL MPI_WAITALL(WTALL,REQ,STATUS_ARRAY,IERR)
      IF(IERR.NE.MPI_SUCCESS) 
     @     CALL WRNDIE(-4,'<GRECSEN>','MPI WAITALL ERROR')
C
##ENDIF (buffers)
C
##ELIF IBMSP
      INTEGER ITMP,JTMP,N
      CALL MP_RECV(RBUF,RLEN*8,TO,TYPE,ITMP)
      CALL MP_SEND(SBUF,SLEN*8,TO,TYPE,JTMP)
      CALL MP_WAIT(ITMP,N)
      CALL MP_WAIT(JTMP,N)
C
##ELIF T3D T3E
##IF SHMEM
Cc##INCLUDE '~/charmm_fcm/shmem.fcm'
      integer i
Cc slen is measured in R*8???
      do i = 1, slen
         sh_send( i ) = sbuf( i )
      end do
Cc len now undefined!?
      call shmem_put( sh_recv, sh_send, len, to )
      call barrier()
cdir$ suppress
      do i = 1, rlen
         rbuf( i ) = sh_recv( i )
      end do
##ELSE
##IF T3D
      INCLUDE "/usr/include/mpp/fpvm3.h"
##ELSE 
      INCLUDE "fpvm3.h"
##ENDIF
      INTEGER IOK,STRIDE
      integer slent,rlent
      PARAMETER (STRIDE=1)

c      rlent = rlen/8
c      slent = slen/8
c      call shmem_get_f(rbuf,rbuf,rlent,to)
c      call shmem_put_f(sbuf,sbuf,slent,to)

      CALL PVMFINITSEND( PVMINPLACE, IOK)
      CALL PVMFPACK( BYTE1, SBUF, SLEN*8, 1, IOK)
      CALL PVMFSEND( TO, TYPE, IOK)

      CALL PVMFRECV(-1, TYPE, IOK)
      CALL PVMFUNPACK(BYTE1, RBUF, RLEN*8, STRIDE, IOK)

      CALL PSYNC()
##ENDIF
##ELIF CM5
      INCLUDE '/usr/include/cm/cmmd_fort.h'
      INTEGER ITMP,JTMP
      ITMP=CMMD_SEND_AND_RECEIVE(TO,TYPE,RBUF,RLEN*8,
     $                           TO,TYPE,SBUF,SLEN*8)
C
##ELIF INTEL
##IF DELTA
      INCLUDE 'fmesh.h'
##ELIF PARAGON
      INCLUDE 'fnx.h'
##ELSE
      INCLUDE 'fcube.h'
##ENDIF
      BYTE DUM
      INTEGER ITMP
      ITMP=IRECV(TYPE+FORCE_TYPE,RBUF,RLEN*8)
      CALL CSEND(TYPE,DUM,0,TO,0)
      CALL CRECV(TYPE,DUM,0)
      CALL CSEND(TYPE+FORCE_TYPE,SBUF,SLEN*8,TO,0)
      CALL MSGWAIT(ITMP)
##ELIF SOCKET
      CALL WRNDIE(-5,'<GRECSEN>','Asynch. sockets not supported yet.')
##ELIF TERRA
      CALL WRNDIE(-5,'<GRECSEN>','Asynch. communication not supported.')
##ELIF PVMC
##INCLUDE '~/charmm_fcm/parallel.fcm'
      INCLUDE 'pvm/fpvm3.h'
      INTEGER BUFID,INFO,RBUFID
C
##IFN ALPHAMP
      CALL PVMFPROBE(TIDS(TO),TYPE,RBUFID)
      IF(RBUFID.NE.0)THEN
         CALL PVMFNRECV(TIDS(TO),TYPE,RBUFID)
         CALL PVMFUNPACK(BYTE1,RBUF,RLEN*8,1,INFO)
      ENDIF
##ENDIF
C
##IF SGIMP CSPP ALPHAMP
      CALL PVMFINITSEND(PVMRAW,BUFID)
##ELSE
      CALL PVMFINITSEND(PVMDATAINPLACE,BUFID)
##ENDIF
      CALL PVMFPACK(BYTE1,SBUF,SLEN*8,1,INFO)
      CALL PVMFSEND(TIDS(TO),TYPE,BUFID)
C
 1    CONTINUE
      CALL PVMFPROBE(TIDS(TO),TYPE,RBUFID)
      IF(RBUFID.NE.0)THEN
         CALL PVMFNRECV(TIDS(TO),TYPE,RBUFID)
         CALL PVMFUNPACK(BYTE1,RBUF,RLEN*8,1,INFO)
      ELSE
         GOTO 1
      ENDIF
##ENDIF
##IF COMMEASURE
      NSEND=NSEND+SLEN*8
      NRECV=NRECV+RLEN*8
      LATENCY=LATENCY+2.0D0
##ENDIF
C
      RETURN
      END
C
      SUBROUTINE GRECSENR(ME,IPPMAP,NUMNOD,ITYPE,W,JLEN,X,ILEN)
C-----------------------------------------------------------------------
C
C     General RECeive and SENd routine. Its purpose is to ease replacement
C     for various platforms (Intel, CM-5, SP1, TCP/IP, ...)
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
      INTEGER ME,IPPMAP(0:*),NUMNOD,ITYPE,ILEN,JLEN
      REAL*8 W(*),X(*)
C
##IF TIMESTAMP
##INCLUDE '~/charmm_fcm/timer.fcm'
##ENDIF
C      
##IF COMMEASURE
      REAL*8 NSEND,NRECV,LATENCY
      COMMON/NTSTAT/NSEND,NRECV,LATENCY
##ENDIF
C
      INTEGER LEFT,RIGHT
C
##IF MPI
C
      INCLUDE 'mpi/mpif.h'
##IF INTEGER8
      INTEGER*4 ITMP,REQ(2),STATUS_ARRAY(MPI_STATUS_SIZE,2),IERR,
     &          MILEN,MITYPE,MJLEN,WTALL,MDP,MCOMM
##ELSE
      INTEGER ITMP,REQ(2),STATUS_ARRAY(MPI_STATUS_SIZE,2),IERR,
     &        MILEN,MITYPE,MJLEN,WTALL,MDP,MCOMM
##ENDIF
##IF BUFFERED (buffersr)
##INCLUDE '~/charmm_fcm/lobos.fcm'
      INTEGER NBUFF,ILENB,JLENB,ILEFT,JLEFT,IPT,JPT,IL,JL,I   !##.not.INTEGER8
      INTEGER*4 NBUFF,ILENB,JLENB,ILEFT,JLEFT,IPT,JPT,IL,JL,I !##INTEGER8
C

      MDP = MPI_DOUBLE_PRECISION
      MCOMM = MPI_COMM_WORLD
      MITYPE=ITYPE

      LEFT=IPPMAP(MOD(ME+NUMNOD-1,NUMNOD))
      RIGHT=IPPMAP(MOD(ME+1,NUMNOD))
C
C     Some networks perform better with the fixed small buffer size
C
      NBUFF=(MAX(ILEN,JLEN)-1)/NICBUF+1
      ILENB=ILEN/NBUFF
      JLENB=JLEN/NBUFF
      ILEFT=ILEN-ILENB*NBUFF
      JLEFT=JLEN-JLENB*NBUFF
      IPT=1
      JPT=1
      DO I=1,NBUFF
         IL=ILENB
         IF(I.LE.ILEFT)IL=IL+1
         JL=JLENB
         IF(I.LE.JLEFT)JL=JL+1
         IF((IL.NE.0).AND.(JL.NE.0)) THEN
##IF TIMESTAMP
            IF(IAND(ISTMPL,4).NE.0) CALL TIME_STAMP(ME,
     $           'GRECSENR>before irecv,R='//CHAR(0),JL*8)
##ENDIF
            CALL MPI_IRECV(W(JPT),JL,MDP,LEFT,MITYPE,
     $           MCOMM,REQ(1),IERR)
            IF(IERR.NE.MPI_SUCCESS) 
     @           CALL WRNDIE(-4,'<GRECSENR>','MPI IRECV ERROR')
            CALL MPI_ISEND(X(IPT),IL,MDP,RIGHT,MITYPE,
     $           MCOMM,REQ(2),IERR)
            IF(IERR.NE.MPI_SUCCESS) 
     @           CALL WRNDIE(-4,'<GRECSENR>','MPI ISEND ERROR')
            WTALL=2
            CALL MPI_WAITALL(WTALL,REQ,STATUS_ARRAY,IERR)
            IF(IERR.NE.MPI_SUCCESS) 
     @           CALL WRNDIE(-4,'<GRECSENR>','MPI WAITALL ERROR')
##IF TIMESTAMP
            IF(IAND(ISTMPL,4).NE.0) CALL TIME_STAMP(ME,
     $           'GRECSENR>after waitall,S='//CHAR(0),IL*8)
##ENDIF
         ELSE IF((IL.NE.0).AND.(JL.EQ.0)) THEN
##IF TIMESTAMP
            IF(IAND(ISTMPL,4).NE.0) CALL TIME_STAMP(ME,
     $           'GRECSENR>before isend,S='//CHAR(0),IL*8)
##ENDIF
            CALL MPI_ISEND(X(IPT),IL,MDP,RIGHT,MITYPE,
     $           MCOMM,REQ(1),IERR)
            IF(IERR.NE.MPI_SUCCESS) 
     @           CALL WRNDIE(-4,'<GRECSENR>','MPI ISEND ERROR')
            WTALL=1
            CALL MPI_WAITALL(WTALL,REQ,STATUS_ARRAY,IERR)
##IF TIMESTAMP
            IF(IAND(ISTMPL,4).NE.0) CALL TIME_STAMP(ME,
     $           'GRECSENR>after waitall,S='//CHAR(0),IL*8)
##ENDIF
            IF(IERR.NE.MPI_SUCCESS) 
     @           CALL WRNDIE(-4,'<GRECSENR>','MPI WAITALL ERROR')
         ELSE IF((IL.EQ.0).AND.(JL.NE.0)) THEN
##IF TIMESTAMP
            IF(IAND(ISTMPL,4).NE.0) CALL TIME_STAMP(ME,
     $           'GRECSENR>before irecv,R='//CHAR(0),JL*8)
##ENDIF
            CALL MPI_IRECV(W(JPT),JL,MDP,LEFT,MITYPE,
     $           MCOMM,REQ(1),IERR)
            IF(IERR.NE.MPI_SUCCESS) 
     @           CALL WRNDIE(-4,'<GRECSENR>','MPI IRECV ERROR')
            WTALL=1
            CALL MPI_WAITALL(WTALL,REQ,STATUS_ARRAY,IERR)
##IF TIMESTAMP
            IF(IAND(ISTMPL,4).NE.0) CALL TIME_STAMP(ME,
     $           'GRECSENR>after waitall,R='//CHAR(0),JL*8)
##ENDIF
            IF(IERR.NE.MPI_SUCCESS) 
     @           CALL WRNDIE(-4,'<GRECSENR>','MPI WAITALL ERROR')
         ENDIF
         IPT=IPT+IL
         JPT=JPT+JL
      ENDDO
##ELSE (buffersr)
C

      MDP = MPI_DOUBLE_PRECISION
      MCOMM = MPI_COMM_WORLD
      MILEN=ILEN
      MJLEN=JLEN
      MITYPE=ITYPE

      LEFT=IPPMAP(MOD(ME+NUMNOD-1,NUMNOD))
      RIGHT=IPPMAP(MOD(ME+1,NUMNOD))
C
##IF TIMESTAMP
      IF(IAND(ISTMPL,4).NE.0) CALL TIME_STAMP(ME,
     $     'GRECSENR>before irecv,R='//CHAR(0),JLEN*8)
##ENDIF
      CALL MPI_IRECV(W,MJLEN,MDP,LEFT,MITYPE,
     $     MCOMM,REQ(1),IERR)
      IF(IERR.NE.MPI_SUCCESS) 
     @     CALL WRNDIE(-4,'<GRECSENR>','MPI IRECV ERROR')
      CALL MPI_ISEND(X,MILEN,MDP,RIGHT,MITYPE,
     $     MCOMM,REQ(2),IERR)
      IF(IERR.NE.MPI_SUCCESS) 
     @     CALL WRNDIE(-4,'<GRECSENR>','MPI ISEND ERROR')
      WTALL=2
      CALL MPI_WAITALL(WTALL,REQ,STATUS_ARRAY,IERR)
##IF TIMESTAMP
      IF(IAND(ISTMPL,4).NE.0) CALL TIME_STAMP(ME,
     $     'GRECSENR>after waitall,S='//CHAR(0),ILEN*8)
##ENDIF
      IF(IERR.NE.MPI_SUCCESS) 
     @     CALL WRNDIE(-4,'<GRECSENR>','MPI WAITALL ERROR')
##ENDIF (buffersr)
C
##ELIF IBMSP
      INTEGER ITMP,JTMP,N
      LEFT=IPPMAP(MOD(ME+NUMNOD-1,NUMNOD))
      RIGHT=IPPMAP(MOD(ME+1,NUMNOD))
      CALL MP_RECV(W,8*JLEN,LEFT,ITYPE,ITMP)
      CALL MP_SEND(X,8*ILEN,RIGHT,ITYPE,JTMP)
      CALL MP_WAIT(ITMP,N)
      CALL MP_WAIT(JTMP,N)
C
##ELSE
      CALL WRNDIE(-5,'<GRECSENR>','Async. ring not supported yet.')
##ENDIF
##IF COMMEASURE
      NSEND=NSEND+ILEN*8
      NRECV=NRECV+JLEN*8
      LATENCY=LATENCY+1.0D0
##ENDIF

      RETURN
      END
C
      SUBROUTINE SWAPD3(N,X,XL,XR,PL,PR)
C-----------------------------------------------------------------------
C
C     Double receive & send routine. Works in 2 direction concurrently
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
##INCLUDE '~/charmm_fcm/stream.fcm'
C
      INTEGER N,PL,PR
      REAL*8 X(*),XL(*),XR(*)
      INTEGER TYPE    !##.not.INTEGER8
      INTEGER*4 TYPE  !##INTEGER8
C
##IF COMMEASURE
      REAL*8 NSEND,NRECV,LATENCY
      COMMON/NTSTAT/NSEND,NRECV,LATENCY
##ENDIF
C
##IF CMPI
##IF MPI
      INCLUDE 'mpi/mpif.h'
##IF INTEGER8
      INTEGER*4 REQ(4),STATUS_ARRAY(MPI_STATUS_SIZE,4),IERR,WTALL,
     &          MN,MPL,MPR,MDP,MCOMM
##ELSE
      INTEGER REQ(4),STATUS_ARRAY(MPI_STATUS_SIZE,4),IERR,WTALL,
     &          MN,MPL,MPR,MDP,MCOMM
##ENDIF
C
      MN=N
      MPL=PL
      MPR=PR
      TYPE=1
      CALL MPI_ISEND(X,MN,MDP,MPL,TYPE,
     $     MCOMM,REQ(1),IERR)
      IF(IERR.NE.MPI_SUCCESS) 
     @     CALL WRNDIE(-4,'<GRECSEN>','MPI ISEND ERROR')
      CALL MPI_IRECV(XL,MN,MDP,MPL,TYPE,
     $     MCOMM,REQ(2),IERR)
      IF(IERR.NE.MPI_SUCCESS) 
     @     CALL WRNDIE(-4,'<GRECSEN>','MPI IRECV ERROR')
      CALL MPI_ISEND(X,MN,MDP,MPR,TYPE,
     $     MCOMM,REQ(1),IERR)
      IF(IERR.NE.MPI_SUCCESS) 
     @     CALL WRNDIE(-4,'<GRECSEN>','MPI ISEND ERROR')
      CALL MPI_IRECV(XR,MN,MDP,MPR,TYPE,
     $     MCOMM,REQ(2),IERR)
      IF(IERR.NE.MPI_SUCCESS) 
     @     CALL WRNDIE(-4,'<GRECSEN>','MPI IRECV ERROR')
      WTALL=4
      CALL MPI_WAITALL(WTALL,REQ,STATUS_ARRAY,IERR)
      IF(IERR.NE.MPI_SUCCESS) 
     @     CALL WRNDIE(-4,'<GRECSEN>','MPI WAITALL ERROR')
C
##ELSE
      CALL WRNDIE(-5,'<SWAPD3>','Not compiled with MPI.')
##ENDIF
##ELSE
      CALL WRNDIE(-5,'<SWAPD3>','Not compiled with CMPI.')
##ENDIF
##IF COMMEASURE
      NSEND=NSEND+SLEN*8
      NRECV=NRECV+RLEN*8
      LATENCY=LATENCY+2.0D0
##ENDIF
C
      RETURN
      END
C
      SUBROUTINE NXYMSH(NX,NY)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
C     Returns number of rows and columns in mesh. Valid for 2**N
C
      INTEGER NX, NY, IDIM, NPTWO
      EXTERNAL NPTWO
C
      IDIM=NPTWO()
      NX=2**(IDIM/2)
      NY=2**(IDIM-IDIM/2)
      RETURN
      END
C
      SUBROUTINE MESH(MYNOD,IPPMAP)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
C     Remaps cube network node numbers to mesh network
C
      INTEGER  MYNOD,IPPMAP(0:*)
      INTEGER I,NODDIM,NUMNOD
      INTEGER  MNOD, NNOD
      EXTERNAL MNOD, NNOD
C
      INTEGER NX,NY
      INTEGER DIMX,DIMY,IBASE,STOD,STEV
      INTEGER INOD,J,K,IBIT,BIT(0:15),BWT(0:15)
C
      CALL NXYMSH(NX,NY)
      NUMNOD=NX*NY
      IF(NUMNOD.NE.NNOD()) CALL DIE
C
      DIMX=-1
      DIMY=-1
      DO I = 0,16
         IF(2**I.EQ.NX) DIMX=I
         IF(2**I.EQ.NY) DIMY=I
      ENDDO
      IF(DIMX.LT.0) CALL DIE
      IF(DIMY.LT.0) CALL DIE
      NODDIM=DIMX+DIMY
C
      STEV=0
      STOD=-DIMX
C
      DO I = DIMX+DIMY,15
         BWT(I)=0
      ENDDO
C
      DO I = DIMX+DIMY-1,0,-1
         IF(STOD.LT.0) THEN
            BWT(I)=STEV
            STEV=STEV+1
            IF(STOD.GT.-DIMX-DIMY) STOD=-STOD
         ELSE
            BWT(I)=STOD
            STOD=STOD+1
            IF(STEV.LT.DIMX) STOD=-STOD
         ENDIF
      ENDDO
C
CC      WRITE(6,47) (BWT(I),I=13,0,-1)
CC  47  FORMAT(' BTW:  ',16I5)
C
      DO I = 0,DIMX+DIMY-1
         BWT(I)=2**BWT(I)
      ENDDO
C
CC      WRITE(6,47) (BWT(I),I=13,0,-1)
C
      MYNOD=-1
      DO INOD = 0,NUMNOD-1
         DO IBIT = 0,15
            IBASE=2**IBIT
            BIT(IBIT)=(MOD(INOD,IBASE*2)-MOD(INOD,IBASE))/IBASE
         ENDDO
CC         WRITE(6,48) INOD,(BIT(J),J=15,0,-1)
CC  48     FORMAT(' INOD,BIT: ',I5,5X,16I2)
         J=0
         DO K = 0,15
            J=J+BWT(K)*BIT(K)
         ENDDO
CC         WRITE(6,49) INOD,J
CC  49     FORMAT(' INOD,J:  ',2I8)
         IPPMAP(INOD)=J
         IF(J.EQ.MNOD()) MYNOD=INOD
      ENDDO
      IF(MYNOD.LT.0) CALL DIE
C
      RETURN
      END
C
      SUBROUTINE AGVDGB(X,IPARPT,NUMNOD,MYNOD,IPPMAP)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
C     This routine performs Actual General Vector Distributed Global Broadcast
C     using communication between neighbours only
C
C     This is synchronized version i.e. it can use synchronized communication.
C     Useful for TCP/IP and like.
C
C     X      - vector to be broadcasted [input&output]
C     IPARPT - vector of length NUMNODES+1 for offsets on each node [input]
C     NUMNOD - number of nodes [input]
C     MYNOD  - who am I [input]
C     IPPMAP - mapping numbers of nodes (cube -> mesh) [1 to 1 on cube] [input]
C
C
C     It is based on ideas from Robert van de Geijn's
C     papers and is optimal in the amount of communication to be performed
C     but needs NUMNOD-1 startups (latency!)
C
      REAL*8 X(*)
      INTEGER IPARPT(0:*), IPPMAP(0:*), MYNOD, NUMNOD
      INTEGER ITYPE,IBIT,ISTEP,ME,IOFF,JOFF
      INTEGER ILEN,JLEN,ISYNC
C
      ITYPE=1
      ME=MYNOD
      ISYNC=MOD(ME,2)
      DO ISTEP = 1, NUMNOD-1
         IBIT=NUMNOD-ISTEP+ME
         IOFF=MOD(IBIT+1,NUMNOD)
         JOFF=MOD(IBIT,NUMNOD)
         ILEN=IPARPT(IOFF+1)-IPARPT(IOFF)
         JLEN=IPARPT(JOFF+1)-IPARPT(JOFF)
         IF(ISYNC.EQ.1) THEN
            CALL GRECL(IPPMAP,ITYPE,X(IPARPT(JOFF)+1),8*JLEN)
            CALL GSENR(IPPMAP,ITYPE,X(IPARPT(IOFF)+1),8*ILEN)
         ELSE          
            CALL GSENR(IPPMAP,ITYPE,X(IPARPT(IOFF)+1),8*ILEN)
            CALL GRECL(IPPMAP,ITYPE,X(IPARPT(JOFF)+1),8*JLEN)
         ENDIF
         ITYPE=ITYPE+1
      ENDDO
C
      RETURN
      END
C
      SUBROUTINE AGVDGBA(X,IPARPT,NUMNOD,MYNOD,IPPMAP)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
C     This routine performs:
C     Actual General Vector Distributed Global Broadcast Asynchronous
C     using communication between neighbours only
C
C     This is asynchronized version
C
C     X      - vector to be broadcasted [input&output]
C     IPARPT - vector of length NUMNODES+1 for offsets on each node [input]
C     NUMNOD - number of nodes [input]
C     MYNOD  - who am I [input]
C     IPPMAP - mapping numbers of nodes (cube -> mesh) [1 to 1 on cube] [input]
C
C
C     It is based on ideas from Robert van de Geijn's
C     papers and is optimal in the amount of communication to be performed
C     but needs NUMNOD-1 startups (latency!)
C
      REAL*8 X(*)
      INTEGER IPARPT(0:*), IPPMAP(0:*), MYNOD, NUMNOD
      INTEGER ITYPE,IBIT,ISTEP,ME,IOFF,JOFF
      INTEGER ILEN,JLEN
C
      ITYPE=1
      ME=MYNOD
      DO ISTEP = 1, NUMNOD-1
         IBIT=NUMNOD-ISTEP+ME
         IOFF=MOD(IBIT+1,NUMNOD)
         JOFF=MOD(IBIT,NUMNOD)
         ILEN=IPARPT(IOFF+1)-IPARPT(IOFF)
         JLEN=IPARPT(JOFF+1)-IPARPT(JOFF)
         CALL GRECSENR(ME,IPPMAP,NUMNOD,ITYPE,X(IPARPT(JOFF)+1),JLEN,
     $        X(IPARPT(IOFF)+1),ILEN)
         ITYPE=ITYPE+1
      ENDDO
C
      RETURN
      END
C
      SUBROUTINE IGVDGB(X,IPARPT,NUMNOD,MYNOD,IPPMAP)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
C     This routine performs Integer General Vector Distributed Global Broadcast
C     using communication between neighbours only. 
C
C     This is synchronized version i.e. it can use synchronized communication.
C     Useful for TCP/IP and like.
C
C     X      - vector to be broadcasted [input&output]
C     IPARPT - vector of length NUMNODES+1 for offsets on each node [input]
C     NUMNOD - number of nodes [input]
C     MYNOD  - who am I [input]
C     IPPMAP - mapping numbers of nodes (cube -> mesh) [1 to 1 on cube] [input]
C
C
C     It is based on ideas from Robert van de Geijn's
C     papers and is optimal in the amount of communication to be performed
C     but needs NUMNOD-1 startups (latency!)
C
      INTEGER X(*)
      INTEGER IPARPT(0:*), IPPMAP(0:*), MYNOD, NUMNOD
      INTEGER ITYPE,IBIT,ISTEP,ME,IOFF,JOFF
      INTEGER ILEN,JLEN,ISYNC
C
      ITYPE=1
      ME=MYNOD
      ISYNC=MOD(ME,2)
      DO ISTEP = 1, NUMNOD-1
         IBIT=NUMNOD-ISTEP+ME
         IOFF=MOD(IBIT+1,NUMNOD)
         JOFF=MOD(IBIT,NUMNOD)
         ILEN=IPARPT(IOFF+1)-IPARPT(IOFF)
         JLEN=IPARPT(JOFF+1)-IPARPT(JOFF)
##IF INTEGER8
         IF(ISYNC.EQ.1) THEN
            CALL GRECL(IPPMAP,ITYPE,X(IPARPT(JOFF)+1),8*JLEN)
            CALL GSENR(IPPMAP,ITYPE,X(IPARPT(IOFF)+1),8*ILEN)
         ELSE          
            CALL GSENR(IPPMAP,ITYPE,X(IPARPT(IOFF)+1),8*ILEN)
            CALL GRECL(IPPMAP,ITYPE,X(IPARPT(JOFF)+1),8*JLEN)
         ENDIF
##ELSE
         IF(ISYNC.EQ.1) THEN
            CALL GRECL(IPPMAP,ITYPE,X(IPARPT(JOFF)+1),4*JLEN)
            CALL GSENR(IPPMAP,ITYPE,X(IPARPT(IOFF)+1),4*ILEN)
         ELSE          
            CALL GSENR(IPPMAP,ITYPE,X(IPARPT(IOFF)+1),4*ILEN)
            CALL GRECL(IPPMAP,ITYPE,X(IPARPT(JOFF)+1),4*JLEN)
         ENDIF
##ENDIF
         ITYPE=ITYPE+1
      ENDDO
C
      RETURN
      END
C
      SUBROUTINE AGVDGS(X,IPARPT,W,NUMNOD,MYNOD,IPPMAP)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
C     This routine performs General Actual Vector Distributed Global Sum
C     using communication between neighbours only
C
C     This is synchronized version i.e. it can use synchronized communication.
C     Useful for TCP/IP and others.
C
C     X      - vector to be broadcasted [input&output]
C     IPARPT - vector of length NUMNODES+1 for offsets on each node [input]
C     NUMNOD - number of nodes [input]
C     MYNOD  - who am I [input]
C     IPPMAP - mapping numbers of nodes (cube -> mesh) [1 to 1 on cube] [input]
C
C     It is based on ideas from Robert van de Geijn's
C     papers and is optimal in the amount of communication to be performed
C     but needs NUMNOD-1 startups (latency!)
C
      REAL*8 X(*), W(*), ONE
      INTEGER IPARPT(0:*), IPPMAP(0:*), MYNOD, NUMNOD
      INTEGER ITYPE,IBIT,ISTEP,ME,IOFF,JOFF
      INTEGER ILEN,JLEN,ISYNC
C
      DATA ONE/1.0D0/
C
      ITYPE=1
      ME=MYNOD
      ISYNC=MOD(ME,2)
      DO ISTEP = 1, NUMNOD-1
         IBIT=NUMNOD-ISTEP+ME
         IOFF=MOD(IBIT,NUMNOD)
         JOFF=MOD(IBIT-1,NUMNOD)
         ILEN=IPARPT(IOFF+1)-IPARPT(IOFF)
         JLEN=IPARPT(JOFF+1)-IPARPT(JOFF)
         IF(ISYNC.EQ.1) THEN
            CALL GRECL(IPPMAP,ITYPE,W(IPARPT(JOFF)+1),8*JLEN)
            CALL GSENR(IPPMAP,ITYPE,X(IPARPT(IOFF)+1),8*ILEN)
         ELSE
            CALL GSENR(IPPMAP,ITYPE,X(IPARPT(IOFF)+1),8*ILEN)
            CALL GRECL(IPPMAP,ITYPE,W(IPARPT(JOFF)+1),8*JLEN)
         ENDIF
         ITYPE=ITYPE+1
         CALL DAXPY_s(JLEN,ONE,W(IPARPT(JOFF)+1),1,X(IPARPT(JOFF)+1),1)
      ENDDO
      RETURN
      END
C
      SUBROUTINE AGVDGSA(X,IPARPT,W,NUMNOD,MYNOD,IPPMAP)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
C     This routine performs:
C     Actual General Vector Distributed Global Sum Asynchronous
C     using communication between neighbours only
C
C     This is asynchronous version i.e. it can communicate in both
C     directions at once.
C
C     X      - vector to be broadcasted [input&output]
C     IPARPT - vector of length NUMNODES+1 for offsets on each node [input]
C     NUMNOD - number of nodes [input]
C     MYNOD  - who am I [input]
C     IPPMAP - mapping numbers of nodes (cube -> mesh) [1 to 1 on cube] [input]
C
C     It is based on ideas from Robert van de Geijn's
C     papers and is optimal in the amount of communication to be performed
C     but needs NUMNOD-1 startups (latency!)
C
      REAL*8 X(*), W(*), ONE
      INTEGER IPARPT(0:*), IPPMAP(0:*), MYNOD, NUMNOD
      INTEGER ITYPE,IBIT,ISTEP,ME,IOFF,JOFF
      INTEGER ILEN,JLEN
C
      DATA ONE/1.0D0/
C
      ITYPE=1
      ME=MYNOD
      DO ISTEP = 1, NUMNOD-1
         IBIT=NUMNOD-ISTEP+ME
         IOFF=MOD(IBIT,NUMNOD)
         JOFF=MOD(IBIT-1,NUMNOD)
         ILEN=IPARPT(IOFF+1)-IPARPT(IOFF)
         JLEN=IPARPT(JOFF+1)-IPARPT(JOFF)
         CALL GRECSENR(ME,IPPMAP,NUMNOD,ITYPE,W(IPARPT(JOFF)+1),JLEN,
     $        X(IPARPT(IOFF)+1),ILEN)
         ITYPE=ITYPE+1
         CALL DAXPY_s(JLEN,ONE,W(IPARPT(JOFF)+1),1,X(IPARPT(JOFF)+1),1)
      ENDDO
      RETURN
      END
C
      SUBROUTINE IGVDGO(X,IPARPT,W,NUMNOD,MYNOD,IPPMAP)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
C     This routine performs Integer General Actual Vector Distributed Global Or
C     using communication between neighbours only
C
C     This is synchronized version i.e. it can use synchronized communication.
C     Useful for TCP/IP and others.
C
C     X      - vector to be broadcasted [input&output]
C     IPARPT - vector of length NUMNODES+1 for offsets on each node [input]
C     NUMNOD - number of nodes [input]
C     MYNOD  - who am I [input]
C     IPPMAP - mapping numbers of nodes (cube -> mesh) [1 to 1 on cube] [input]
C
C     It is based on ideas from Robert van de Geijn's
C     papers and is optimal in the amount of communication to be performed
C     but needs NUMNOD-1 startups (latency!)
C
      INTEGER X(*), W(*)
      INTEGER IPARPT(0:*), IPPMAP(0:*), MYNOD, NUMNOD
      INTEGER ITYPE,IBIT,ISTEP,ME,IOFF,JOFF
      INTEGER ILEN,JLEN,ISYNC
C
      ITYPE=1
      ME=MYNOD
      ISYNC=MOD(ME,2)
      DO ISTEP = 1, NUMNOD-1
         IBIT=NUMNOD-ISTEP+ME
         IOFF=MOD(IBIT,NUMNOD)
         JOFF=MOD(IBIT-1,NUMNOD)
         ILEN=IPARPT(IOFF+1)-IPARPT(IOFF)
         JLEN=IPARPT(JOFF+1)-IPARPT(JOFF)
##IF INTEGER8
         IF(ISYNC.EQ.1) THEN
            CALL GRECL(IPPMAP,ITYPE,W(IPARPT(JOFF)+1),8*JLEN)
            CALL GSENR(IPPMAP,ITYPE,X(IPARPT(IOFF)+1),8*ILEN)
         ELSE          
            CALL GSENR(IPPMAP,ITYPE,X(IPARPT(IOFF)+1),8*ILEN)
            CALL GRECL(IPPMAP,ITYPE,W(IPARPT(JOFF)+1),8*JLEN)
         ENDIF
##ELSE
         IF(ISYNC.EQ.1) THEN
            CALL GRECL(IPPMAP,ITYPE,W(IPARPT(JOFF)+1),4*JLEN)
            CALL GSENR(IPPMAP,ITYPE,X(IPARPT(IOFF)+1),4*ILEN)
         ELSE          
            CALL GSENR(IPPMAP,ITYPE,X(IPARPT(IOFF)+1),4*ILEN)
            CALL GRECL(IPPMAP,ITYPE,W(IPARPT(JOFF)+1),4*JLEN)
         ENDIF
##ENDIF
         ITYPE=ITYPE+1
         CALL IARROR(JLEN,W(IPARPT(JOFF)+1),X(IPARPT(JOFF)+1))
      ENDDO
      RETURN
      END
C
      SUBROUTINE IARROR(N,W,X)
C-----------------------------------------------------------------------
C
C     Integer ARRay OR routine
C
      INTEGER I, N
      INTEGER W(*), X(*)
      DO I = 1, N
         X(I) = IOR(X(I), W(I))
      ENDDO
      RETURN
      END
C
      SUBROUTINE IGVDGS(X,IPARPT,W,NUMNOD,MYNOD,IPPMAP)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
C     This routine performs Integer General Actual Vector Distributed Global Sum
C     using communication between neighbours only
C
C     This is synchronized version i.e. it can use synchronized communication.
C     Useful for TCP/IP and others.
C
C     X      - vector to be broadcasted [input&output]
C     IPARPT - vector of length NUMNODES+1 for offsets on each node [input]
C     NUMNOD - number of nodes [input]
C     MYNOD  - who am I [input]
C     IPPMAP - mapping numbers of nodes (cube -> mesh) [1 to 1 on cube] [input]
C
C     It is based on ideas from Robert van de Geijn's
C     papers and is optimal in the amount of communication to be performed
C     but needs NUMNOD-1 startups (latency!)
C
      INTEGER X(*), W(*)
      INTEGER IPARPT(0:*), IPPMAP(0:*), MYNOD, NUMNOD
      INTEGER ITYPE,IBIT,ISTEP,ME,IOFF,JOFF
      INTEGER ILEN,JLEN,ISYNC
C
      ITYPE=1
      ME=MYNOD
      ISYNC=MOD(ME,2)
      DO ISTEP = 1, NUMNOD-1
         IBIT=NUMNOD-ISTEP+ME
         IOFF=MOD(IBIT,NUMNOD)
         JOFF=MOD(IBIT-1,NUMNOD)
         ILEN=IPARPT(IOFF+1)-IPARPT(IOFF)
         JLEN=IPARPT(JOFF+1)-IPARPT(JOFF)
##IF GNUALPHA
         IF(ISYNC.EQ.1) THEN
            CALL GRECL(IPPMAP,ITYPE,W(IPARPT(JOFF)+1),8*JLEN)
            CALL GSENR(IPPMAP,ITYPE,X(IPARPT(IOFF)+1),8*ILEN)
         ELSE          
            CALL GSENR(IPPMAP,ITYPE,X(IPARPT(IOFF)+1),8*ILEN)
            CALL GRECL(IPPMAP,ITYPE,W(IPARPT(JOFF)+1),8*JLEN)
         ENDIF
##ELSE
         IF(ISYNC.EQ.1) THEN
            CALL GRECL(IPPMAP,ITYPE,W(IPARPT(JOFF)+1),4*JLEN)
            CALL GSENR(IPPMAP,ITYPE,X(IPARPT(IOFF)+1),4*ILEN)
         ELSE          
            CALL GSENR(IPPMAP,ITYPE,X(IPARPT(IOFF)+1),4*ILEN)
            CALL GRECL(IPPMAP,ITYPE,W(IPARPT(JOFF)+1),4*JLEN)
         ENDIF
##ENDIF
         ITYPE=ITYPE+1
         CALL IARRSUM(JLEN,W(IPARPT(JOFF)+1),X(IPARPT(JOFF)+1))
      ENDDO
      RETURN
      END
C
      SUBROUTINE IARRSUM(N,W,X)
C-----------------------------------------------------------------------
C
C     Integer ARRay SUM routine
C
      INTEGER I, N
      INTEGER W(*), X(*)
      DO I = 1, N
         X(I) = X(I)+W(I)
      ENDDO
      RETURN
      END
C
      SUBROUTINE AVDGBR(X,IPARPT,NDIM,MYNOD,IPPMAP)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
C     This routine performs Actual Vector Distributed Global BRoadcast
C     using recursive doubling scheme.
C
C     This is synchronized version i.e. it can use synchronized communication.
C     Useful for TCP/IP and like.
C
C     X - vector to be broadcasted [input&output]
C     IPARPT - vector of length NUMNODES+1 for offsets on each node [input]
C     NDIM - dimension of hypercube [input]
C     MYNOD - who am I [input]
C     IPPMAP - mapping numbers of nodes (cube -> mesh) [1 to 1 on cube] [input]
C
C     It is based on ideas from book Fox, et al, and Robert van de Geijn's
C     papers & private communications and gives theoretically 
C     best possible performance. 
C
      REAL*8 X(*)
      INTEGER IPARPT(0:*), IPPMAP(0:*), MYNOD, NDIM
      INTEGER ITYPE,IBIT,IDIM,ME,IOFF,JOFF
      INTEGER ILEN,JLEN,NEXT,ISYNC
      INTEGER IEOR,IAND
C
##IF CM5
C     This CM-5 implementation uses active messages block transfers
C     
      INCLUDE '/usr/include/cm/cmmd_fort.h'
C     
      INTEGER RPORT2, itmp
      EXTERNAL rhfunc_sr 
      INTEGER CMAML_number_of_rports
      INTEGER CMAML_set_rport_addr
      INTEGER CMAML_set_rport_byte_counter
      INTEGER CMAML_set_rport_total_nbytes
      INTEGER CMAML_set_rport_handler
      INTEGER CMAML_set_rport_handler_arg2
C     
##INCLUDE '~/charmm_fcm/cm5.fcm'
C
      rport2 = CMAML_number_of_rports() - 2
      itmp = CMAML_disable_interrupts()
      itmp = CMAML_allocate_this_rport(RPORT2)
      if (itmp.lt.1)
     $     CALL WRNDIE(-5,'<AVDGBR>','PORT not available for CMAML.')
C     
      itype = 1
      ibit = 1
      do idim = 1, ndim
         call CMMD_sync_with_nodes()
         next = ieor(mynod, ibit)
         ioff = iand(next, -ibit)
         joff = iand(mynod, -ibit)
         ilen = iparpt(ioff + ibit) - iparpt(ioff)
         jlen = iparpt(joff + ibit) - iparpt(joff)
         if (ilen.ne.0) then
            itmp = CMAML_set_rport_addr(RPORT2, x(iparpt(ioff) + 1))
            itmp = CMAML_set_rport_byte_counter(RPORT2, 8*ilen)
            itmp = CMAML_set_rport_total_nbytes(RPORT2, 8*ilen)
            itmp = CMAML_set_rport_handler(RPORT2, rhfunc_sr)
            itmp = CMAML_set_rport_handler_arg2(RPORT2, 0)
            recv_completed = 0
         else
            recv_completed = 1
         endif
         call CMMD_sync_with_nodes()
         if (jlen.ne.0) then
            itmp = CMAML_scopy(ippmap(next), RPORT2, 0, 0,
     &           x(iparpt(joff) + 1), 8*jlen, 0, 0)
         endif
C     
 10      continue
         call CMAML_poll()
         if (recv_completed.eq.0) goto 10 
C     
         ibit = 2*ibit 
         itype = itype + 1
      enddo
C     
      itmp = CMAML_free_rport(RPORT2)
      itmp = CMAML_enable_interrupts()
C     
##ELSE
      ITYPE=1
      ME=MYNOD
      DO IDIM = 1, NDIM
         IBIT=2**(IDIM-1)
         NEXT=IEOR(ME,IBIT)
         ISYNC=IAND(ME,IBIT)
         IOFF=IAND(NEXT,-IBIT)
         JOFF=IAND(ME,-IBIT)
         ILEN=IPARPT(IOFF+IBIT)-IPARPT(IOFF)
         JLEN=IPARPT(JOFF+IBIT)-IPARPT(JOFF)
         IF(ISYNC.EQ.IBIT) THEN
            CALL GREC(IPPMAP(NEXT),ITYPE,X(IPARPT(IOFF)+1),8*ILEN)
            CALL GSEN(IPPMAP(NEXT),ITYPE,X(IPARPT(JOFF)+1),8*JLEN)
         ELSE
            CALL GSEN(IPPMAP(NEXT),ITYPE,X(IPARPT(JOFF)+1),8*JLEN)
            CALL GREC(IPPMAP(NEXT),ITYPE,X(IPARPT(IOFF)+1),8*ILEN)
         ENDIF
         ITYPE=ITYPE+1
      ENDDO
##ENDIF
C
C Problems with Intel communication
##IF INTEL T3E
      CALL PSYNC()
##ENDIF
C
      RETURN
      END
C
      SUBROUTINE IVDGBR(X,IPARPT,NDIM,MYNOD,IPPMAP)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
C     This routine performs Integer Vector Distributed Global BRoadcast
C     using recursive doubling scheme.
C
C     This is synchronized version i.e. it can use synchronized communication.
C     Useful for TCP/IP and like.
C
C     X - vector to be broadcasted [input&output]
C     IPARPT - vector of length NUMNODES+1 for offsets on each node [input]
C     NDIM - dimension of hypercube [input]
C     MYNOD - who am I [input]
C     IPPMAP - mapping numbers of nodes (cube -> mesh) [1 to 1 on cube] [input]
C
C     It is based on ideas from book Fox, et al, and Robert van de Geijn's
C     papers & private communications and gives theoretically 
C     best possible performance. 
C
      INTEGER X(*)
      INTEGER IPARPT(0:*), IPPMAP(0:*), MYNOD, NDIM
      INTEGER ITYPE,IBIT,IDIM,ME,IOFF,JOFF
      INTEGER ILEN,JLEN,NEXT,ISYNC
      INTEGER IEOR,IAND
C
##IF CM5
C     This CM-5 implementation uses active messages block transfers
C     
      INCLUDE '/usr/include/cm/cmmd_fort.h'
C     
      INTEGER RPORT2, itmp
      EXTERNAL rhfunc_sr 
      INTEGER CMAML_number_of_rports
      INTEGER CMAML_set_rport_addr
      INTEGER CMAML_set_rport_byte_counter
      INTEGER CMAML_set_rport_total_nbytes
      INTEGER CMAML_set_rport_handler
      INTEGER CMAML_set_rport_handler_arg2
C     
##INCLUDE '~/charmm_fcm/cm5.fcm'
C
      rport2 = CMAML_number_of_rports() - 2
      itmp = CMAML_disable_interrupts()
      itmp = CMAML_allocate_this_rport(RPORT2)
      if (itmp.lt.1)
     $     CALL WRNDIE(-5,'<AVDGBR>','PORT not available for CMAML.')
C     
      itype = 1
      ibit = 1
      do idim = 1, ndim
         call CMMD_sync_with_nodes()
         next = ieor(mynod, ibit)
         ioff = iand(next, -ibit)
         joff = iand(mynod, -ibit)
         ilen = iparpt(ioff + ibit) - iparpt(ioff)
         jlen = iparpt(joff + ibit) - iparpt(joff)
         if (ilen.ne.0) then
            itmp = CMAML_set_rport_addr(RPORT2, x(iparpt(ioff) + 1))
            itmp = CMAML_set_rport_byte_counter(RPORT2, 4*ilen)
            itmp = CMAML_set_rport_total_nbytes(RPORT2, 4*ilen)
            itmp = CMAML_set_rport_handler(RPORT2, rhfunc_sr)
            itmp = CMAML_set_rport_handler_arg2(RPORT2, 0)
            recv_completed = 0
         else
            recv_completed = 1
         endif
         call CMMD_sync_with_nodes()
         if (jlen.ne.0) then
            itmp = CMAML_scopy(ippmap(next), RPORT2, 0, 0,
     &           x(iparpt(joff) + 1), 4*jlen, 0, 0)
         endif
C     
 10      continue
         call CMAML_poll()
         if (recv_completed.eq.0) goto 10 
C     
         ibit = 2*ibit 
         itype = itype + 1
      enddo
C     
      itmp = CMAML_free_rport(RPORT2)
      itmp = CMAML_enable_interrupts()
C     
##ELSE
      ITYPE=1
      ME=MYNOD
      DO IDIM = 1, NDIM
         IBIT=2**(IDIM-1)
         NEXT=IEOR(ME,IBIT)
         ISYNC=IAND(ME,IBIT)
         IOFF=IAND(NEXT,-IBIT)
         JOFF=IAND(ME,-IBIT)
         ILEN=IPARPT(IOFF+IBIT)-IPARPT(IOFF)
         JLEN=IPARPT(JOFF+IBIT)-IPARPT(JOFF)
##IF T3D T3E INTEGER8
         IF(ISYNC.EQ.IBIT) THEN
            CALL GREC(IPPMAP(NEXT),ITYPE,X(IPARPT(IOFF)+1),8*ILEN)
            CALL GSEN(IPPMAP(NEXT),ITYPE,X(IPARPT(JOFF)+1),8*JLEN)
         ELSE
            CALL GSEN(IPPMAP(NEXT),ITYPE,X(IPARPT(JOFF)+1),8*JLEN)
            CALL GREC(IPPMAP(NEXT),ITYPE,X(IPARPT(IOFF)+1),8*ILEN)
         ENDIF
##ELSE
         IF(ISYNC.EQ.IBIT) THEN
            CALL GREC(IPPMAP(NEXT),ITYPE,X(IPARPT(IOFF)+1),4*ILEN)
            CALL GSEN(IPPMAP(NEXT),ITYPE,X(IPARPT(JOFF)+1),4*JLEN)
         ELSE
            CALL GSEN(IPPMAP(NEXT),ITYPE,X(IPARPT(JOFF)+1),4*JLEN)
            CALL GREC(IPPMAP(NEXT),ITYPE,X(IPARPT(IOFF)+1),4*ILEN)
         ENDIF
##ENDIF
         ITYPE=ITYPE+1
      ENDDO
##ENDIF
C
C Problems with Intel communication
##IF INTEL T3E
      CALL PSYNC()
##ENDIF
C
      RETURN
      END
C
C
      SUBROUTINE AVDGBRA(X,IPARPT,NDIM,MYNOD,IPPMAP)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
C     This routine performs Actual Vector Distributed Global BRoadcast
C     using recursive doubling scheme.
C
C     This is asynchronous version i.e. it can overlap send and receive
C     operations for maximum speed.
C
C     X - vector to be broadcasted [input&output]
C     IPARPT - vector of length NUMNODES+1 for offsets on each node [input]
C     NDIM - dimension of hypercube [input]
C     MYNOD - who am I [input]
C     IPPMAP - mapping numbers of nodes (cube -> mesh) [1 to 1 on cube] [input]
C
C     It is based on ideas from book Fox, et al, and Robert van de Geijn's
C     papers & private communications and gives theoretically 
C     best possible performance. 
C
      REAL*8 X(*)
      INTEGER IPARPT(0:*), IPPMAP(0:*), MYNOD, NDIM
      INTEGER ITYPE,IBIT,IDIM,ME,IOFF,JOFF
      INTEGER ILEN,JLEN,NEXT
      INTEGER IEOR,IAND
C
      ITYPE=1
      ME=MYNOD
      DO IDIM = 1, NDIM
         IBIT=2**(IDIM-1)
         NEXT=IEOR(ME,IBIT)
         IOFF=IAND(NEXT,-IBIT)
         JOFF=IAND(ME,-IBIT)
         ILEN=IPARPT(IOFF+IBIT)-IPARPT(IOFF)
         JLEN=IPARPT(JOFF+IBIT)-IPARPT(JOFF)
         CALL GRECSEN(IPPMAP(NEXT),ITYPE,X(IPARPT(IOFF)+1),ILEN,
     $                                   X(IPARPT(JOFF)+1),JLEN)
         ITYPE=ITYPE+1
      ENDDO
C
C Problems with Intel communication
##IF INTEL T3E
      CALL PSYNC()
##ENDIF
C
      RETURN
      END
C
      SUBROUTINE AVDGSUM(X,IPARPT,W,NDIM,MYNOD,IPPMAP)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
C     This routine performs Actual Vector Distributed Global SUM
C     using recursive halving scheme.
C
C     This is synchronized version i.e. it can use synchronized communication.
C     Useful for TCP/IP and others.
C
C     X - vector to be broadcasted [input&output]
C     IPARPT - vector of length NUMNODES+1 for offsets on each node [input]
C     NDIM - dimension of hypercube [input]
C     MYNOD - who am I [input]
C     IPPMAP - mapping numbers of nodes (cube -> mesh) [1 to 1 on cube] [input]
C
C     It is based on ideas from book Fox, et al, and Robert van de Geijn's
C     papers & private communications and gives theoretically 
C     best possible performance. 
C
      REAL*8 X(*), W(*), ONE
      INTEGER IPARPT(0:*), IPPMAP(0:*), MYNOD, NDIM
      INTEGER ITYPE,IBIT,IDIM,ME,IOFF,JOFF
      INTEGER ILEN,JLEN,NEXT,ISYNC
      INTEGER IEOR,IAND
C
##IF CM5
C     This CM-5 version uses active messages block transfers (rports)
      
      INCLUDE '/usr/include/cm/cmmd_fort.h'
      
      INTEGER i, RPORT2, itmp
      EXTERNAL rhfunc_sr 
      INTEGER CMAML_number_of_rports
      INTEGER CMAML_set_rport_addr
      INTEGER CMAML_set_rport_byte_counter
      INTEGER CMAML_set_rport_total_nbytes
      INTEGER CMAML_set_rport_handler
      INTEGER CMAML_set_rport_handler_arg2
C     
##INCLUDE '~/charmm_fcm/cm5.fcm'
C     
      rport2 = CMAML_number_of_rports() - 1
      itmp = CMAML_disable_interrupts()
      itmp = CMAML_allocate_this_rport(RPORT2)
      if (itmp.lt.1)
     $     CALL WRNDIE(-5,'<AVDGSUM>','PORT not available for CMAML.')
C     
      itype = 1
      do idim = ndim, 1, -1
         ibit=2**(idim-1)
         call CMMD_sync_with_nodes()
         next = ieor(mynod, ibit)
         ioff = iand(next,-ibit)
         joff = iand(mynod,-ibit)
         ilen = iparpt(ioff+ibit)-iparpt(ioff)
         jlen = iparpt(joff+ibit)-iparpt(joff)
         if (jlen.ne.0) then
            itmp = CMAML_set_rport_addr(RPORT2, w(iparpt(joff) + 1))
            itmp = CMAML_set_rport_byte_counter(RPORT2, 8*jlen)
            itmp = CMAML_set_rport_total_nbytes(RPORT2, 8*jlen)
            itmp = CMAML_set_rport_handler(RPORT2, rhfunc_sr)
            itmp = CMAML_set_rport_handler_arg2(RPORT2, 0)
            recv_completed = 0
         else
            recv_completed = 1
         end if
         call CMMD_sync_with_nodes()
         if (ilen.ne.0) then
            itmp = CMAML_scopy(ippmap(next), RPORT2, 0, 0, 
     &           x(iparpt(ioff) + 1), 8*ilen, 0, 0)
         end if
         
 10      continue
         call CMAML_poll()
         if (recv_completed.eq.0) goto 10 
         
         do i = 1, jlen
            x(iparpt(joff)+i) = x(iparpt(joff)+i) + w(iparpt(joff)+i)
         enddo
         itype = itype + 1
      end do
      
      itmp = CMAML_free_rport(RPORT2)
      itmp = CMAML_enable_interrupts()
C     
##ELSE
      DATA ONE/1.0D0/
C
      ITYPE=1
      ME=MYNOD
      DO IDIM = NDIM, 1, -1
         IBIT=2**(IDIM-1)
         NEXT=IEOR(ME,IBIT)
         ISYNC=IAND(ME,IBIT)
         IOFF=IAND(NEXT,-IBIT)
         JOFF=IAND(ME,-IBIT)
         ILEN=IPARPT(IOFF+IBIT)-IPARPT(IOFF)
         JLEN=IPARPT(JOFF+IBIT)-IPARPT(JOFF)
         IF(ISYNC.EQ.IBIT) THEN
            CALL GREC(IPPMAP(NEXT),ITYPE,W(IPARPT(JOFF)+1),8*JLEN)
            CALL GSEN(IPPMAP(NEXT),ITYPE,X(IPARPT(IOFF)+1),8*ILEN)
         ELSE
            CALL GSEN(IPPMAP(NEXT),ITYPE,X(IPARPT(IOFF)+1),8*ILEN)
            CALL GREC(IPPMAP(NEXT),ITYPE,W(IPARPT(JOFF)+1),8*JLEN)
         ENDIF
         ITYPE=ITYPE+1
         CALL DAXPY_s(JLEN,ONE,W(IPARPT(JOFF)+1),1,X(IPARPT(JOFF)+1),1)
      ENDDO
C
##ENDIF
C Problems with Intel communication
##IF INTEL T3E
      CALL PSYNC()
##ENDIF
C
      RETURN
      END
C
      SUBROUTINE IVDGOR(X,IPARPT,W,NDIM,MYNOD,IPPMAP)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
C     This routine performs Integer Vector Distributed Global OR
C     using recursive halving scheme.
C
C     This is synchronized version i.e. it can use synchronized communication.
C     Useful for TCP/IP and others.
C
C     X - vector to be broadcasted [input&output]
C     IPARPT - vector of length NUMNODES+1 for offsets on each node [input]
C     NDIM - dimension of hypercube [input]
C     MYNOD - who am I [input]
C     IPPMAP - mapping numbers of nodes (cube -> mesh) [1 to 1 on cube] [input]
C
C     It is based on ideas from book Fox, et al, and Robert van de Geijn's
C     papers & private communications and gives theoretically 
C     best possible performance. 
C
      INTEGER X(*), W(*)
      INTEGER IPARPT(0:*), IPPMAP(0:*), MYNOD, NDIM
      INTEGER ITYPE,IBIT,IDIM,ME,IOFF,JOFF
      INTEGER ILEN,JLEN,NEXT,ISYNC
      INTEGER IEOR,IAND
C
##IF CM5
C     This CM-5 version uses active messages block transfers (rports)
      
      INCLUDE '/usr/include/cm/cmmd_fort.h'
      
      INTEGER i, RPORT2, itmp
      EXTERNAL rhfunc_sr 
      INTEGER CMAML_number_of_rports
      INTEGER CMAML_set_rport_addr
      INTEGER CMAML_set_rport_byte_counter
      INTEGER CMAML_set_rport_total_nbytes
      INTEGER CMAML_set_rport_handler
      INTEGER CMAML_set_rport_handler_arg2
C     
##INCLUDE '~/charmm_fcm/cm5.fcm'
C     
      rport2 = CMAML_number_of_rports() - 1
      itmp = CMAML_disable_interrupts()
      itmp = CMAML_allocate_this_rport(RPORT2)
      if (itmp.lt.1)
     $     CALL WRNDIE(-5,'<AVDGSUM>','PORT not available for CMAML.')
C     
      itype = 1
      do idim = ndim, 1, -1
         ibit=2**(idim-1)
         call CMMD_sync_with_nodes()
         next = ieor(mynod, ibit)
         ioff = iand(next,-ibit)
         joff = iand(mynod,-ibit)
         ilen = iparpt(ioff+ibit)-iparpt(ioff)
         jlen = iparpt(joff+ibit)-iparpt(joff)
         if (jlen.ne.0) then
            itmp = CMAML_set_rport_addr(RPORT2, w(iparpt(joff) + 1))
            itmp = CMAML_set_rport_byte_counter(RPORT2, 4*jlen)
            itmp = CMAML_set_rport_total_nbytes(RPORT2, 4*jlen)
            itmp = CMAML_set_rport_handler(RPORT2, rhfunc_sr)
            itmp = CMAML_set_rport_handler_arg2(RPORT2, 0)
            recv_completed = 0
         else
            recv_completed = 1
         end if
         call CMMD_sync_with_nodes()
         if (ilen.ne.0) then
            itmp = CMAML_scopy(ippmap(next), RPORT2, 0, 0, 
     &           x(iparpt(ioff) + 1), 4*ilen, 0, 0)
         end if
         
 10      continue
         call CMAML_poll()
         if (recv_completed.eq.0) goto 10 
         
         do i = 1, jlen
            x(iparpt(joff)+i) = ior(x(iparpt(joff)+i),w(iparpt(joff)+i))
         enddo
         itype = itype + 1
      end do
      
      itmp = CMAML_free_rport(RPORT2)
      itmp = CMAML_enable_interrupts()
C     
##ELSE
C
      ITYPE=1
      ME=MYNOD
      DO IDIM = NDIM, 1, -1
         IBIT=2**(IDIM-1)
         NEXT=IEOR(ME,IBIT)
         ISYNC=IAND(ME,IBIT)
         IOFF=IAND(NEXT,-IBIT)
         JOFF=IAND(ME,-IBIT)
         ILEN=IPARPT(IOFF+IBIT)-IPARPT(IOFF)
         JLEN=IPARPT(JOFF+IBIT)-IPARPT(JOFF)
##IF T3D T3E INTEGER8
         IF(ISYNC.EQ.IBIT) THEN
            CALL GREC(IPPMAP(NEXT),ITYPE,W(IPARPT(JOFF)+1),8*JLEN)
            CALL GSEN(IPPMAP(NEXT),ITYPE,X(IPARPT(IOFF)+1),8*ILEN)
         ELSE
            CALL GSEN(IPPMAP(NEXT),ITYPE,X(IPARPT(IOFF)+1),8*ILEN)
            CALL GREC(IPPMAP(NEXT),ITYPE,W(IPARPT(JOFF)+1),8*JLEN)
         ENDIF
##ELSE
         IF(ISYNC.EQ.IBIT) THEN
            CALL GREC(IPPMAP(NEXT),ITYPE,W(IPARPT(JOFF)+1),4*JLEN)
            CALL GSEN(IPPMAP(NEXT),ITYPE,X(IPARPT(IOFF)+1),4*ILEN)
         ELSE
            CALL GSEN(IPPMAP(NEXT),ITYPE,X(IPARPT(IOFF)+1),4*ILEN)
            CALL GREC(IPPMAP(NEXT),ITYPE,W(IPARPT(JOFF)+1),4*JLEN)
         ENDIF
##ENDIF
         ITYPE=ITYPE+1
         CALL IARROR(JLEN,W(IPARPT(JOFF)+1),X(IPARPT(JOFF)+1))
      ENDDO
C
##ENDIF
C Problems with Intel communication
##IF INTEL T3E
      CALL PSYNC()
##ENDIF
C
      RETURN
      END
C
      SUBROUTINE IVDGSUM(X,IPARPT,W,NDIM,MYNOD,IPPMAP)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
C     This routine performs Integer Vector Distributed Global SUM
C     using recursive halving scheme.
C
C     This is synchronized version i.e. it can use synchronized communication.
C     Useful for TCP/IP and others.
C
C     X - vector to be broadcasted [input&output]
C     IPARPT - vector of length NUMNODES+1 for offsets on each node [input]
C     NDIM - dimension of hypercube [input]
C     MYNOD - who am I [input]
C     IPPMAP - mapping numbers of nodes (cube -> mesh) [1 to 1 on cube] [input]
C
C     It is based on ideas from book Fox, et al, and Robert van de Geijn's
C     papers & private communications and gives theoretically 
C     best possible performance. 
C
      INTEGER X(*), W(*)
      INTEGER IPARPT(0:*), IPPMAP(0:*), MYNOD, NDIM
      INTEGER ITYPE,IBIT,IDIM,ME,IOFF,JOFF
      INTEGER ILEN,JLEN,NEXT,ISYNC
      INTEGER IEOR,IAND
C
C
      ITYPE=1
      ME=MYNOD
      DO IDIM = NDIM, 1, -1
         IBIT=2**(IDIM-1)
         NEXT=IEOR(ME,IBIT)
         ISYNC=IAND(ME,IBIT)
         IOFF=IAND(NEXT,-IBIT)
         JOFF=IAND(ME,-IBIT)
         ILEN=IPARPT(IOFF+IBIT)-IPARPT(IOFF)
         JLEN=IPARPT(JOFF+IBIT)-IPARPT(JOFF)
##IF T3D T3E GNUALPHA
         IF(ISYNC.EQ.IBIT) THEN
            CALL GREC(IPPMAP(NEXT),ITYPE,W(IPARPT(JOFF)+1),8*JLEN)
            CALL GSEN(IPPMAP(NEXT),ITYPE,X(IPARPT(IOFF)+1),8*ILEN)
         ELSE
            CALL GSEN(IPPMAP(NEXT),ITYPE,X(IPARPT(IOFF)+1),8*ILEN)
            CALL GREC(IPPMAP(NEXT),ITYPE,W(IPARPT(JOFF)+1),8*JLEN)
         ENDIF
##ELSE
         IF(ISYNC.EQ.IBIT) THEN
            CALL GREC(IPPMAP(NEXT),ITYPE,W(IPARPT(JOFF)+1),4*JLEN)
            CALL GSEN(IPPMAP(NEXT),ITYPE,X(IPARPT(IOFF)+1),4*ILEN)
         ELSE
            CALL GSEN(IPPMAP(NEXT),ITYPE,X(IPARPT(IOFF)+1),4*ILEN)
            CALL GREC(IPPMAP(NEXT),ITYPE,W(IPARPT(JOFF)+1),4*JLEN)
         ENDIF
##ENDIF
         ITYPE=ITYPE+1
         CALL IARRSUM(JLEN,W(IPARPT(JOFF)+1),X(IPARPT(JOFF)+1))
      ENDDO
C
      RETURN
      END
C
      SUBROUTINE AVDGSA(X,IPARPT,W,NDIM,MYNOD,IPPMAP)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
C     This routine performs Actual Vector Distributed Global SUM
C     using recursive halving scheme.
C
C     This is asynchronous version for concurent send/recieve
C
C     X - vector to be broadcasted [input&output]
C     IPARPT - vector of length NUMNODES+1 for offsets on each node [input]
C     NDIM - dimension of hypercube [input]
C     MYNOD - who am I [input]
C     IPPMAP - mapping numbers of nodes (cube -> mesh) [1 to 1 on cube] [input]
C
C     It is based on ideas from book Fox, et al, and Robert van de Geijn's
C     papers & private communications and gives theoretically 
C     best possible performance. 
C
      REAL*8 X(*), W(*), ONE
      INTEGER IPARPT(0:*), IPPMAP(0:*), MYNOD, NDIM
      INTEGER ITYPE,IBIT,IDIM,ME,IOFF,JOFF
      INTEGER ILEN,JLEN,NEXT
      INTEGER IEOR,IAND
C
      DATA ONE/1.0D0/
C
      ITYPE=1
      ME=MYNOD
      DO IDIM = NDIM, 1, -1
         IBIT=2**(IDIM-1)
         NEXT=IEOR(ME,IBIT)
         IOFF=IAND(NEXT,-IBIT)
         JOFF=IAND(ME,-IBIT)
         ILEN=IPARPT(IOFF+IBIT)-IPARPT(IOFF)
         JLEN=IPARPT(JOFF+IBIT)-IPARPT(JOFF)
         CALL GRECSEN(IPPMAP(NEXT),ITYPE,W(IPARPT(JOFF)+1),JLEN,
     $                                   X(IPARPT(IOFF)+1),ILEN)
         ITYPE=ITYPE+1
         CALL DAXPY_s(JLEN,ONE,W(IPARPT(JOFF)+1),1,X(IPARPT(JOFF)+1),1)
      ENDDO
C
C Problems with Intel communication
##IF INTEL T3E
      CALL PSYNC()
##ENDIF
C
      RETURN
      END
C
      SUBROUTINE DAXPY_s(N,ONE,W,IONE,X,JONE)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
C     Simplified BLAS routine with the same name
C
      INTEGER I, N, IONE, JONE
      REAL*8 ONE, W(*), X(*)
      DO I = 1, N
         X(I) = X(I) + W(I)
      ENDDO
      RETURN
      END
C
C=======================
C     Block communication for force decomposition....
C     Maybe we don't put PARASCAL keyword here because this
C     routines are used also in small communication only test
C     programs without any changes
C
C=======================
##IF PARASCAL (force_decomp)
      SUBROUTINE BLOCBRA(IB,X,IPARPT,NDIM,MYNOD,IPPMAP,IPMAT,JPMAT)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
C     This routine performs Actual Vector Distributed Global BRoadcast
C     using recursive doubling scheme.
C
C     This is asynchronous version i.e. it can overlap send and receive
C     operations for maximum speed.
C
C     X - vector to be broadcasted [input&output]
C     IPARPT - vector of length NUMNODES+1 for offsets on each node [input]
C     NDIM - dimension of hypercube [input]
C     MYNOD - who am I [input]
C     IPPMAP - mapping numbers of nodes (cube -> mesh) [1 to 1 on cube] [input]
C     IPMAT - block organization
C     JPMAT - mapping of MYNOD to ME which is valid within one block.
C
C     It is based on ideas from book Fox, et al, and Robert van de Geijn's
C     papers & private communications and gives theoretically 
C     best possible performance. 
C
      REAL*8 X(*)
      INTEGER IPARPT(0:*), IPPMAP(0:*), MYNOD, NDIM
      INTEGER IB, IPMAT(4,*),JPMAT(4,0:*)
      INTEGER ITYPE,IBIT,IDIM,ME,IOFF,JOFF
      INTEGER ILEN,JLEN,NEXT
      INTEGER IEOR,IAND
C
      ITYPE=IB*NDIM
      ME=JPMAT(IB,MYNOD)
C     Am I in this block?
      IF(ME.EQ.-1) RETURN
      DO IDIM = 1, NDIM
         IBIT=2**(IDIM-1)
         NEXT=IEOR(ME,IBIT)
         IOFF=IAND(NEXT,-IBIT)
         JOFF=IAND(ME,-IBIT)
         ILEN=IPARPT(IOFF+IBIT)-IPARPT(IOFF)
         JLEN=IPARPT(JOFF+IBIT)-IPARPT(JOFF)
         CALL GRECSEN(IPPMAP(IPMAT(IB,NEXT+1)),ITYPE,X(IPARPT(IOFF)+1),
     $                ILEN,X(IPARPT(JOFF)+1),JLEN)
         ITYPE=ITYPE+1
      ENDDO
C
      RETURN
      END
C
      SUBROUTINE BLOCKBRO(IB,X,IPARPT,NDIM,MYNOD,IPPMAP,IPMAT,JPMAT)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
C     This routine performs Actual Vector Distributed Global BRoadcast
C     using recursive doubling scheme.
C
C     This is version for force decomposition, where vectors are organized in
C     blocks and each block is further divided on parts. This is very similar
C     to AVDGBR with some additional mapping.
C
C
C     This is synchronized version i.e. it can use synchronized communication.
C     Useful for TCP/IP and like.
C
C     X - vector to be broadcasted [input&output]
C     IPARPT - vector of length NUMNODES+1 for offsets on each node [input]
C     NDIM - dimension of hypercube [input]
C     MYNOD - who am I [input]
C     IPPMAP - mapping numbers of nodes (cube -> mesh) [1 to 1 on cube] [input]
C     IPMAT - block organization
C     JPMAT - mapping of MYNOD to ME which is valid within one block.
C
C     It is based on ideas from book Fox, et al, and Robert van de Geijn's
C     papers & private communications and gives theoretically 
C     best possible performance. 
C
      REAL*8 X(*)
      INTEGER IPARPT(0:*), IPPMAP(0:*), MYNOD, NDIM
      INTEGER IB, IPMAT(4,*),JPMAT(4,0:*)
      INTEGER ITYPE,IBIT,IDIM,ME,IOFF,JOFF
      INTEGER ILEN,JLEN,NEXT,ISYNC
      INTEGER IEOR,IAND
C
      ITYPE=IB*NDIM
      ME=JPMAT(IB,MYNOD)
C     Am I in this block?
      IF(ME.EQ.-1) RETURN
      DO IDIM = 1, NDIM
         IBIT=2**(IDIM-1)
         NEXT=IEOR(ME,IBIT)
         ISYNC=IAND(ME,IBIT)
         IOFF=IAND(NEXT,-IBIT)
         JOFF=IAND(ME,-IBIT)
         ILEN=IPARPT(IOFF+IBIT)-IPARPT(IOFF)
         JLEN=IPARPT(JOFF+IBIT)-IPARPT(JOFF)
         IF(ISYNC.EQ.IBIT) THEN
            CALL GREC(IPPMAP(IPMAT(IB,NEXT+1)),ITYPE,
     $           X(IPARPT(IOFF)+1),8*ILEN)
            CALL GSEN(IPPMAP(IPMAT(IB,NEXT+1)),ITYPE,
     $           X(IPARPT(JOFF)+1),8*JLEN)
         ELSE
            CALL GSEN(IPPMAP(IPMAT(IB,NEXT+1)),ITYPE,
     $           X(IPARPT(JOFF)+1),8*JLEN)
            CALL GREC(IPPMAP(IPMAT(IB,NEXT+1)),ITYPE,
     $           X(IPARPT(IOFF)+1),8*ILEN)
         ENDIF
         ITYPE=ITYPE+1
      ENDDO
C
      RETURN
      END
C
      SUBROUTINE BLKBRA(X,IPARPT,NPART,MNSTEP,ISTEP,IDIM,MYNOD,
     $                  ISCHED,JSCHED,IPMAT,IPPMAP)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
C     This routine performs Actual Vector Distributed Global SUM
C     using recursive halving scheme.
C
C     This is asynchronous version for concurent send/receive
C
C     X - vector to be broadcasted [input&output]
C     IPARPT - vector of length NUMNODES+1 for offsets on each node [input]
C     NDIM - dimension of hypercube [input]
C     MYNOD - who am I [input]
C     IPPMAP - mapping numbers of nodes (cube -> mesh) [1 to 1 on cube] [input]
C     NPART - How many blocks
C     ISTEP - communication step
C     IPMAT - block organization
C     ISCHED - mapping of MYNOD to ME which is valid within one block
C              on a given step.
C     JSCHED - To which block MYNOD processor belong on a given step.
C
C
      REAL*8 X(*)
      INTEGER IPARPT(0:*),IPPMAP(0:*),MYNOD,IDIM,ISTEP,NPART,MNSTEP
      INTEGER IB,IPMAT(8,*),ISCHED(MNSTEP,0:*),JSCHED(MNSTEP,0:*)
      INTEGER ITYPE,IBIT,ME,IOFF,JOFF
      INTEGER ILEN,JLEN,NEXT
      INTEGER IEOR,IAND
C
      ME=ISCHED(ISTEP,MYNOD)
C     Am I in this step?
      IF(ME.EQ.-1) RETURN
      IB=JSCHED(ISTEP,MYNOD)
C
      ITYPE=IB
C
      IBIT=2**(IDIM-1)
      NEXT=IEOR(ME,IBIT)
      IOFF=IAND(NEXT,-IBIT)+(IB-1)*NPART
      JOFF=IAND(ME,-IBIT)+(IB-1)*NPART
      ILEN=IPARPT(IOFF+IBIT)-IPARPT(IOFF)
      JLEN=IPARPT(JOFF+IBIT)-IPARPT(JOFF)
      CALL GRECSEN(IPPMAP(IPMAT(IB,NEXT+1)),ITYPE,X(IPARPT(IOFF)+1),
     $             ILEN,X(IPARPT(JOFF)+1),JLEN)
C
      RETURN
      END
C
      SUBROUTINE BLKSA(X,IPARPT,W,NPART,MNSTEP,ISTEP,IDIM,MYNOD,
     $                 ISCHED,JSCHED,IPMAT,IPPMAP)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
C     This routine performs Actual Vector Distributed Global SUM
C     using recursive halving scheme.
C
C     This is asynchronous version for concurent send/receive
C
C     X - vector to be broadcasted [input&output]
C     IPARPT - vector of length NUMNODES+1 for offsets on each node [input]
C     NDIM - dimension of hypercube [input]
C     MYNOD - who am I [input]
C     IPPMAP - mapping numbers of nodes (cube -> mesh) [1 to 1 on cube] [input]
C     NPART - How many blocks
C     ISTEP - communication step
C     IPMAT - block organization
C     ISCHED - mapping of MYNOD to ME which is valid within one block
C              on a given step.
C     JSCHED - To which block MYNOD processor belong on a given step.
C
C
      REAL*8 X(*), W(*), ONE
      INTEGER IPARPT(0:*),IPPMAP(0:*),MYNOD,IDIM,ISTEP,NPART,MNSTEP
      INTEGER IB,IPMAT(8,*),ISCHED(MNSTEP,0:*),JSCHED(MNSTEP,0:*)
      INTEGER ITYPE,IBIT,ME,IOFF,JOFF
      INTEGER ILEN,JLEN,NEXT
      INTEGER IEOR,IAND
C
      DATA ONE/1.0D0/
C
      ME=ISCHED(ISTEP,MYNOD)
C     Am I in this step?
      IF(ME.EQ.-1) RETURN
      IB=JSCHED(ISTEP,MYNOD)
C
      ITYPE=IB
C
      IBIT=2**(IDIM-1)
      NEXT=IEOR(ME,IBIT)
      IOFF=IAND(NEXT,-IBIT)+(IB-1)*NPART
      JOFF=IAND(ME,-IBIT)+(IB-1)*NPART
      ILEN=IPARPT(IOFF+IBIT)-IPARPT(IOFF)
      JLEN=IPARPT(JOFF+IBIT)-IPARPT(JOFF)
      CALL GRECSEN(IPPMAP(IPMAT(IB,NEXT+1)),ITYPE,W(IPARPT(JOFF)+1),
     $             JLEN,X(IPARPT(IOFF)+1),ILEN)
      CALL DAXPY_s(JLEN,ONE,W(IPARPT(JOFF)+1),1,X(IPARPT(JOFF)+1),1)
C
      RETURN
      END
C
C
      SUBROUTINE BLOCKSA(IB,X,IPARPT,W,NDIM,MYNOD,IPPMAP,IPMAT,JPMAT)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
C     This routine performs Actual Vector Distributed Global SUM
C     using recursive halving scheme.
C
C     This is asynchronous version for concurent send/recieve
C
C     X - vector to be broadcasted [input&output]
C     IPARPT - vector of length NUMNODES+1 for offsets on each node [input]
C     NDIM - dimension of hypercube [input]
C     MYNOD - who am I [input]
C     IPPMAP - mapping numbers of nodes (cube -> mesh) [1 to 1 on cube] [input]
C     IPMAT - block organization
C     JPMAT - mapping of MYNOD to ME which is valid within one block.
C
C     It is based on ideas from book Fox, et al, and Robert van de Geijn's
C     papers & private communications and gives theoretically 
C     best possible performance. 
C
      REAL*8 X(*), W(*), ONE
      INTEGER IPARPT(0:*), IPPMAP(0:*), MYNOD, NDIM
      INTEGER IB, IPMAT(4,4),JPMAT(4,0:9)
      INTEGER ITYPE,IBIT,IDIM,ME,IOFF,JOFF
      INTEGER ILEN,JLEN,NEXT
      INTEGER IEOR,IAND
C
      DATA ONE/1.0D0/
C
      ITYPE=IB*NDIM
      ME=JPMAT(IB,MYNOD)
C     Am I in this block?
      IF(ME.EQ.-1) RETURN
      DO IDIM = NDIM, 1, -1
         IBIT=2**(IDIM-1)
         NEXT=IEOR(ME,IBIT)
         IOFF=IAND(NEXT,-IBIT)
         JOFF=IAND(ME,-IBIT)
         ILEN=IPARPT(IOFF+IBIT)-IPARPT(IOFF)
         JLEN=IPARPT(JOFF+IBIT)-IPARPT(JOFF)
         CALL GRECSEN(IPPMAP(IPMAT(IB,NEXT+1)),ITYPE,W(IPARPT(JOFF)+1),
     $                JLEN,X(IPARPT(IOFF)+1),ILEN)
         ITYPE=ITYPE+1
         CALL DAXPY_s(JLEN,ONE,W(IPARPT(JOFF)+1),1,X(IPARPT(JOFF)+1),1)
      ENDDO
      RETURN
      END
C
      SUBROUTINE BLOCKSUM(IB,X,IPARPT,W,NDIM,MYNOD,IPPMAP,IPMAT,JPMAT)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
C
C     This routine performs Actual Vector Distributed Global SUM
C     using recursive halving scheme. 
C
C     This is version for force decomposition, where vectors are organized in
C     blocks and each block is further divided on parts. This is very similar
C     to AVDGSUM with some additional mapping.
C
C     This is synchronized version i.e. it can use synchronized communication.
C     Useful for TCP/IP and others.
C
C     X - vector to be broadcasted [input&output]
C     IPARPT - vector of length NUMNODES+1 for offsets on each node [input]
C     NDIM - dimension of hypercube [input]
C     MYNOD - who am I [input]
C     IPPMAP - mapping numbers of nodes (cube -> mesh) [1 to 1 on cube] [input]
C     IPMAT - block organization
C     JPMAT - mapping of MYNOD to ME which is valid within one block.
C
C     It is based on ideas from book Fox, et al, and Robert van de Geijn's
C     papers & private communications and gives theoretically 
C     best possible performance on hypercube architecture.
C
      REAL*8 X(*), W(*), ONE
      INTEGER IPARPT(0:*), IPPMAP(0:*), MYNOD, NDIM
      INTEGER IB, IPMAT(4,4),JPMAT(4,0:9)
      INTEGER ITYPE,IBIT,IDIM,ME,IOFF,JOFF
      INTEGER ILEN,JLEN,NEXT,ISYNC
      INTEGER IEOR,IAND
C
      DATA ONE/1.0D0/
C
      ITYPE=IB*NDIM
      ME=JPMAT(IB,MYNOD)
C     Am I in this block?
      IF(ME.EQ.-1) RETURN
CCCC      write(*,*)'BLOCKSUM>me,mynod,ib',me,mynod,ib
      
      DO IDIM = NDIM, 1, -1
         IBIT=2**(IDIM-1)
         NEXT=IEOR(ME,IBIT)
         ISYNC=IAND(ME,IBIT)
         IOFF=IAND(NEXT,-IBIT)
         JOFF=IAND(ME,-IBIT)
         ILEN=IPARPT(IOFF+IBIT)-IPARPT(IOFF)
         JLEN=IPARPT(JOFF+IBIT)-IPARPT(JOFF)
         IF(ISYNC.EQ.IBIT) THEN
            CALL GREC(IPPMAP(IPMAT(IB,NEXT+1)),ITYPE,
     $           W(IPARPT(JOFF)+1),8*JLEN)
            CALL GSEN(IPPMAP(IPMAT(IB,NEXT+1)),ITYPE,
     $           X(IPARPT(IOFF)+1),8*ILEN)
         ELSE
            CALL GSEN(IPPMAP(IPMAT(IB,NEXT+1)),ITYPE,
     $           X(IPARPT(IOFF)+1),8*ILEN)
            CALL GREC(IPPMAP(IPMAT(IB,NEXT+1)),ITYPE,
     $           W(IPARPT(JOFF)+1),8*JLEN)
         ENDIF
         ITYPE=ITYPE+1
         CALL DAXPY_s(JLEN,ONE,W(IPARPT(JOFF)+1),1,X(IPARPT(JOFF)+1),1)
      ENDDO
      RETURN
      END
C
##ENDIF (force_decomp)
C
      SUBROUTINE APSYNG()
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
##INCLUDE '~/charmm_fcm/number.fcm'
C
C     Actual Global Sync routine General 
C     Just do Global sum of one number
C
      REAL*8 DUMMY,WDUMMY
C
      DUMMY=ONE
      CALL GGSUM0(1,DUMMY,WDUMMY)
C
      RETURN
      END
C
      SUBROUTINE APSYNH()
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
C     Actual Global Sync routine General (sending one byte around)
C
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
C
      INTEGER DUMMY,NEXT,IBIT,TYPE,ME
C
C     Has to send at least something to really block the program
C
      TYPE=100
      ME=MYNOD
      IBIT=1
 10   CONTINUE
         IF (ME.LT.IBIT) THEN
            NEXT=ME+IBIT
            IF (NEXT.LT.NUMNOD) THEN 
               CALL GSEN(IPPMAP(NEXT),TYPE,DUMMY,1)
            ENDIF
         ELSE IF (ME.LT.2*IBIT) THEN
            NEXT=ME-IBIT
            CALL GREC(IPPMAP(NEXT),TYPE,DUMMY,1)
         ENDIF
         IBIT=IBIT*2
         TYPE=TYPE+1
      IF (IBIT.LT.NUMNOD) GOTO 10
C
      RETURN
      END
C
      SUBROUTINE GGSUM0(N,X,W)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
C     General Global SUM (valid for any number of nodes)
C
C     Inefficient and simple, but it is true global sum
C
      INTEGER N
      REAL*8 X(*),W(*)
C
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
C
      INTEGER R1, R2, S, TYPE, IONE, JONE
      REAL*8 ONE
      DATA IONE/1/, JONE/1/, ONE/1.0D0/
C
      TYPE=1
      R1=2*MYNOD+1
      R2=R1+1
      S=(MYNOD-1)/2
C
      IF(R1.LT.NUMNOD) THEN
         CALL GREC(IPPMAP(R1),TYPE,W,8*N)
         CALL DAXPY_S(N,ONE,W,IONE,X,JONE)
      ENDIF
C
      IF(R2.LT.NUMNOD) THEN
         CALL GREC(IPPMAP(R2),TYPE,W,8*N)
         CALL DAXPY_S(N,ONE,W,IONE,X,JONE)
      ENDIF
C
      IF(MYNOD.NE.0) THEN
         CALL GSEN(IPPMAP(S),TYPE,X,8*N)
      ENDIF
C
      CALL GBRDCS(N,X)
C
      RETURN
      END
C
      SUBROUTINE GBRDCS(N,X)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
C     General BRoaDCaSt routine from node 0 
C      (valid for any number of nodes)
C
      INTEGER N
      REAL*8 X(*)
C
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
C
      INTEGER NEXT,IBIT,TYPE,ME,I,IMAX
C
      TYPE=1
      ME=MYNOD
      IBIT=1
      IMAX=0
      IF(NUMNOD.GT.2**NODDIM)IMAX=1
CCC      DO I = 1, NUMNOD/2
      DO I = 1,NODDIM+IMAX
         IF (ME.LT.IBIT) THEN
            NEXT=ME+IBIT
CCC            IF (NEXT.LT.NUMNOD) CALL GSEN(IPPMAP(NEXT),TYPE,X,8*N)
            IF (NEXT.LT.NUMNOD) CALL GSEN(IPPMAP(NEXT),TYPE,X,N)
         ELSE IF (ME.LT.2*IBIT) THEN
            NEXT=ME-IBIT
CCC            CALL GREC(IPPMAP(NEXT),TYPE,X,8*N)
            CALL GREC(IPPMAP(NEXT),TYPE,X,N)
         ENDIF
         IBIT=IBIT*2
      ENDDO
C
      RETURN
C
C=========================================
C End of BLOCK stuff
C
C===================================
##ENDIF ! PARALLEL
      END

