CHARMM Element source/machdep/paral1.src 1.1
C Utility routines for parallel CHARMM on MIMD machines.
C
C Machines currently supported:
C       Intel iPSC/860 (gamma-Touchstone)  (## INTEL)
C       Intel delta-Touchstone (both ## INTEL and DELTA)
C       Intel Paragon (both ## INTEL and PARAGON)
C       Ethernet cluster (diverse workstations) (## ETHER)
C       Token Ring (Apollo workstations) (## RING)
C       FDDI Ring or Ethernet (homogeneous workstations) (## SOCKET)
C       Thinking Machines CM-5 (## CM5)
C       IBM SPn (## IBMSP)
C       PVM with any platform which supports it.  (## PVM)
C           (NOTE: on majority of workstation clusters it is
C                  recommended to use ##SOCKET instead of ##PVM,
C                  but they are both fully functional)
C           PowerChallenge (PVM-3.3.7) doesn't support PVMDataInPlace
C                          - use (## SGIMP)
C       Convex Exemplar SPP-1000 (## CSPP) Also doesn't support
C                                          PVMDataInPlace
C       Cray - T3D (## T3D)
C       TERRA 2000 (## TERRA)
C
C Synchronous (one way) I/O controlled by (## SYNCHRON)
C
C Entire section controlled by ##IF PARALLEL
C
C Written by Bernie Brooks & Milan Hodoscek,
C          Original version From April, 1991
C          New Global sum routines from April 1992
C          (Based on methods by G. Fox and R. Van de Geijn)
C          (Intel iPSC/860 assistance from Stan Erwin)
C
C Completely rewritten for easy support of many platforms, June 1993
C CM5 support revamped, July 1994
C Some reorganization to better support different libraries on the same
C platform, July 1997
C
C This file is organized in 3 layers:
C 1. layer: routines which are called from the rest of CHARMM
C           This is implemented as direct call to MPI library
C           In case of bad implementation CMPI calls should be used.
C 2. layer: machine independent implementation of CMPI calls
C 3. layer: machine dependent implementation of layer 2
C For details see comments under LAYER labels below.
C
C NEW: This file was split in 3 - paral1.fcm, paral2.fcm, paral3.fcm
C      each LAYER is in its own file
C
C brbrooks@helix.nih.gov, milan@par10.mgsl.dcrt.nih.gov
C nagle@tammy.harvard.edu, tamayo@think.com (CM5)
C
C Please check with B. Brooks before publication for current
C appropriate acknowledgement.  Thanks.
C
##IFN PARALLEL
      SUBROUTINE PARCMD(COMLYN,COMLEN)
      CHARACTER COMLYN*(*)
      INTEGER   COMLEN
      CALL WRNDIE(-1,'<PARCMD>','Parallel code not compiled.')
      RETURN
      END
      SUBROUTINE NULL_PL
      RETURN
##ELSE
C
      SUBROUTINE PARCMD(COMLYN,COMLEN)
C-----------------------------------------------------------------------
C
C     PARAllel command parser for controlling parallel execution
C     Syntax:
C
C     PARAllel CONCurrent <int> ...
C
C         CONCurrent <int>   specify how many concurrent jobs
C                            to run in the system
C
C     PARAllel FIFO <int>    specify FIFO scheduler in LoBoS with
C                            static priority <int>
C
C     PARAllel BUFF <int>    specify buffer size for send/receive 
C                            calls. <int> is in REAL*8 units
C
C     PARAllel RPATH <nrep>  CYCL
C              RPATH <nrep>  specify the real number of replication
C                            for the whole system
C              CYCL          Flag for CYCLyc RPATH command (must match RPATH CYCL)
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
##INCLUDE '~/charmm_fcm/stream.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
##INCLUDE '~/charmm_fcm/lobos.fcm'
C
C . Passed variables.
      CHARACTER COMLYN*(*)
      INTEGER   COMLEN
      INTEGER I,MEOLD,NOLD
      INTEGER CLEN,IPGRP,ISTPRI
      CHARACTER*100 CCOM
      INTEGER LIREP,HIREP,NTREP
      LOGICAL QCYCL
C
##IF CONCURR (concurrent)
CCC      write(*,*)'PARCMD>ME-OLD,comlen,comlyn=',mynod,comlen,
CCC     $     comlyn(1:comlen)
      IPCURR=GTRMI(COMLYN,COMLEN,'CONC',1)
      IF((IPCURR.GT.1).AND.(PRNLEV.GE.2)) THEN
         WRITE(OUTU,'(A,I2,A)')'This job will run in ',IPCURR,
     $        ' separate groups.'
      ENDIF
      INTVCN=GTRMI(COMLYN,COMLEN,'COUN',-1)
      INTVMX=GTRMI(COMLYN,COMLEN,'MAXI',0)
      CALL PSND4(IPCURR,1)
      CALL PSND4(INTVCN,1)
      CALL PSND4(INTVMX,1)
      IF(IPCURR.EQ.1) GOTO 10
C
      MEOLD=MYNOD
      NOLD=NUMNOD
      CALL PARSEP(MYNOD,NUMNOD,IPPMAP,IPMGRP,IPCURR)
      CALL SETMSI('PARGROUP',IPMGRP)
C
      MYNODP = MYNOD+1
C
      PRNLEV=5
      WRNLEV=5
      IOLEV=1
      plnod0=prnlev
      IF(MYNOD.NE.0) plnod0=0
      call gbor(plnod0,1)
##IFN CM5
      IF(MYNOD.NE.0) THEN
         PRNLEV=-1
         WRNLEV=-5
         IOLEV=-1
      ENDIF
##ENDIF ! CM5
      DPCOMM=1
C
      DO I=1,MAXNODE
         INODE(I)=MOD(I,NUMNOD)
      ENDDO
      write(*,*)'PARCMD>ME-OLD,ME-NEW,NUMNOD-O,NUMNOD-N,IPMGRP=',
     $           meold,mynod,nold,numnod,ipmgrp
      write(*,*)'PARCMD>ME-OLD,prnlev,wrnlev,iolev,ipcurr,inode=',
     $     meold,prnlev,wrnlev,iolev,ipcurr,(inode(i),i=1,numnod)
C
      IF(IPMGRP.EQ.0) THEN
         CCOM='1.input'
         CLEN=7
      ELSE IF(IPMGRP.EQ.1) THEN
         CCOM='2.input'
         CLEN=7
      ELSE
         CALL WRNDIE(-5,'<PARCMD>','Only 2 groups allowed!')
      ENDIF
      CALL PUSTRM(CCOM,100,CLEN)
 10   CONTINUE
##ENDIF (concurrent)
C
##IF REPLICA
      NTREP=GTRMI(COMLYN,COMLEN,'RPAT',0)
      QCYCL=(INDXA(COMLYN,COMLEN,'CYCL').GT.0)
      IF(NTREP.GT.0)THEN
         CALL SETMSI('LIREP',LIREP)
         CALL SETMSI('HIREP',HIREP)
         CALL SETMSI('NTREP',NTREP)
         IF(PRNLEV.GE.2)THEN
            WRITE(*,*)'Replica path method will be used with',NTREP,
     $           ' replications'
            WRITE(*,*)'using',numnod,' procesors'
         ENDIF
      ENDIF
##ENDIF
C
##IF SCHED
C     This code is for Linux "real time" scheduling
C     with the static priority ISTPRI
      ISTPRI=GTRMI(COMLYN,COMLEN,'FIFO',0)
      CALL PSND4(ISTPRI,1)
      IF(ISTPRI.GT.0)CALL RTSCHED(ISTPRI)
##ENDIF
##IF BUFFERED
C     This code is to specify buffer size for send/receive
      NICBUF=GTRMI(COMLYN,COMLEN,'BUFF',2000)
      CALL PSND4(NICBUF,1)
##ENDIF
C
      RETURN
      END
C
C-----------------------------------------------------------------------
C Startup for parallel code.
C-----------------------------------------------------------------------
      SUBROUTINE PARSTRT
C  This code must be run when CHARMM is started to properly setup all
C  data communication.
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
      INTEGER I
C
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
##INCLUDE '~/charmm_fcm/stream.fcm'
##INCLUDE '~/charmm_fcm/number.fcm'
##IF TIMESTAMP
##INCLUDE '~/charmm_fcm/timer.fcm'
##ENDIF
C
##IF CMPI
CCC   This could go to parallel.fcm protected by CMPI ??
      INTEGER MPI_COMM_WORLD
      PARAMETER (MPI_COMM_WORLD=0)
##ELIF TERRA
      INCLUDE 'Trtlib_f.h'
##ELIF IBMSP
      INCLUDE  'mpif.h'
##ELSE
      INCLUDE 'mpi/mpif.h'
##ENDIF
##IF SOCKET
      INTEGER WDMAX
      PARAMETER (WDMAX=80)
      CHARACTER*(WDMAX) NAME(16), ENV0
      INTEGER LLEN(16), J, MYPID,LENH,LUNASS
      EXTERNAL LUNASS
      REAL*8 RPIDS(16)
##ENDIF
      INTEGER STATUS
C
C     Initialize
##IF CMPI
      CALL CMPI_INIT(STATUS)
##ELIF TERRA
C
      CALL TRTGETMYPENUMBER(MYNOD)
      CALL TRTGETNUMBERPES(NUMNOD)
      CALL CUBE(MYNOD,NUMNOD,IPPMAP)
C
##ELSE
      CALL MPI_INIT(STATUS)
      CALL MPI_COMM_RANK(MPI_COMM_WORLD,MYNOD,STATUS)
      CALL MPI_COMM_SIZE(MPI_COMM_WORLD,NUMNOD,STATUS)
      CALL MPI_BARRIER(MPI_COMM_WORLD,STATUS)
      CALL CUBE(MYNOD,NUMNOD,IPPMAP)
##ENDIF
C
      MYNODP = MYNOD+1
C
C     Normally there is only one parallel group
      IPCURR = 1
      IPMGRP = 0
      INTVCN=-1
C
##IFN CM5
      plnod0=prnlev
      IF(MYNOD.NE.0) plnod0=0
      call gbor(plnod0,1)
      IF(MYNOD.NE.0) THEN
         PRNLEV=-1
         WRNLEV=-5
         IOLEV=-1
      ENDIF
##ENDIF ! CM5
      DPCOMM=1
C
      CALL SETMSI('MYNODE',MYNOD)
      CALL SETMSI('NUMNODE',NUMNOD)
C
      DO I=1,MAXNODE
         INODE(I)=MOD(I,NUMNOD)
      ENDDO
C
C     Initialize timing array
C
      IMERI=0
      DO I = 1, TMGNUM
         TMERI(I)=ZERO
      ENDDO
##IF TIMESTAMP
      ICMRW=0
      ISTMPL=0
##ENDIF
C
C     Check priority of the main process and assign the same on the
C     rest of processes
##IF GNU
      I=0
      IF(MYNOD.EQ.0) CALL GETPPR(I)
      CALL PSND4(I,1)
      IF(MYNOD.NE.0) CALL SETPPR(I)
##ENDIF
C
##IF SOCKET
C LNI, March 1998. 
C Only needed when using SOCKET comm?
C Get pids of slave processes and save in a file on master node 
C so that runaway slaves can be killed in case of crash.
C File name has to be specified in env variable REMOVE_CH_SLAVES
C This file could be removed in PARFIN in case of a successful ending.
C
      IF (NUMNOD.LE.1) RETURN
      CALL FGETPID(MYPID)
C      write(*,*) '%%%: me&mypid',mynod,mypid
C GCOMB needs REAL*8 
      RPIDS(MYNODP)=MYPID
C      write(*,*)'%%%-to-gcomb:',mynod,NUMNOD,(RPIDS(I),I=1,NUMNOD)
      CALL GCOMB(RPIDS,NUMNOD)
C
      IF(MYNOD.GT.0) RETURN
C      write(*,*)'%%%-from-gcomb:',mynod,NUMNOD,(RPIDS(I),I=1,NUMNOD)
C
C Try to stay cool if there is a write problem, this is non-essential
      CALL FGETENV(ENV0,LENH,'REMOVE_CH_SLAVES'//CHAR(0))
      IF(LENH.LE.0) RETURN
      J=LUNASS(90)
      OPEN(J,NAME=ENV0(1:LENH),ERR=999,STATUS='UNKNOWN')
      CALL ENVNUM(I,NAME,LLEN)
      WRITE(J,'(A)') '#!/bin/sh'
      WRITE(J,'(A,I2,A)') '# script to kill possibly still running',
     &              NUMNOD-1,' slave(s) from'
      WRITE(J,'(A)') '# parallel CHARMM job on '//NAME(1)(1:LLEN(1))
      DO I=1,NUMNOD-1 
         MYPID=RPIDS(I+1)
         WRITE(J,'(A,I10)',ERR=998) 
     &   'rsh '//NAME(I+1)(1:LLEN(I+1))//' kill -9 ',MYPID
      ENDDO
998   CLOSE(UNIT=J,ERR=999)
999   RETURN
##ELSE
C
      RETURN
##ENDIF
      END
C
C-----------------------------------------------------------------------
C LAYER 1
C Routines that are called from the rest of CHARMM and for efficiency
C implemented as MPI or variuos other message passing libraries.
C-----------------------------------------------------------------------
C
      SUBROUTINE VDGSUM(X,Y,Z,MODE)
C-----------------------------------------------------------------------
C
C    Wraping routine for vector distributed global sum.
C
C    MODE=0 - Do just the inner sum for PARSCAL
C         (assume contributions are only on processors containing block)
C    MODE=1 - Do the complete global sum (outer followed by inner)
C             (assume values can be on any node)
C    MODE=2 - Do just the outer sum.
C             (assume values anywhere, but only needed by block nodes)
C
C The MODE values 0 and 1 are identical and MODE 2 is ignored for PARAFULL.
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
      REAL*8 X(*), Y(*), Z(*)
      INTEGER MODE
C
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
C
##INCLUDE '~/charmm_fcm/exfunc.fcm'
##INCLUDE '~/charmm_fcm/stack.fcm'
##IF TIMESTAMP
##INCLUDE '~/charmm_fcm/timer.fcm'
##ENDIF
C
##IF CMPI
      INTEGER MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD
      PARAMETER(MPI_DOUBLE_PRECISION=0,MPI_SUM=0,MPI_COMM_WORLD=0)
##ELIF TERRA
      INCLUDE 'tcommlib_f.h'
##ELIF IBMSP
      INCLUDE  'mpif.h'
##ELSE
      INCLUDE 'mpi/mpif.h'
##ENDIF
C
      INTEGER N, W, STATUS
      INTEGER TAG            !##TIMESTAMP
C
##IF PARAFULL
##IF PARASCAL
##ERROR 'Illegal parallel compile options'
##ENDIF
      IF(MODE.EQ.2) RETURN
##ELIF PARASCAL
##ELSE
##ERROR 'Illegal parallel compile options'
##ENDIF
C
##IF TIMESTAMP
      IF (IAND(ISTMPL,2).NE.0) CALL TIME_STAMP(MYNOD,
     $     'VDGSUM called from CHARMM'//CHAR(0),IPARPT(NUMNOD))
      IF (ICMRW.EQ.1) THEN
         N=IPARPT(NUMNOD)
         CALL GRDFMM(X,N,'VDGSUMX  ')
         CALL GRDFMM(Y,N,'VDGSUMY  ')
         CALL GRDFMM(Z,N,'VDGSUMZ  ')
         RETURN
      ELSE IF (ICMRW.EQ.2) THEN
         N=IPARPT(NUMNOD)
         CALL WRTFMR(X,'VDGSUMX  ',TAG,N)
         CALL WRTFMR(Y,'VDGSUMY  ',TAG,N)
         CALL WRTFMR(Z,'VDGSUMZ  ',TAG,N)
         RETURN
      ENDIF
##ENDIF
##IF CONCURR
      CALL PGRPIN(IPMGRP,MYNOD,NUMNOD,INTVCN,INTVMX)
##ENDIF
      IF(NUMNOD.EQ.1) RETURN
      N=IPARPT(NUMNOD)
      W=ALLSTK(IREAL8(N))
##IF CMPI (Charmm_MPI)
C     Uses CHARMM MPI implementation
      CALL CMPI_RED_SCAT(X,STACK(W),NPARPT,
     @                     MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,
     @                     STATUS)
      CALL COPYR8(STACK(W),X(IPARPT(MYNOD)+1),NPARPT(MYNOD))
      CALL CMPI_RED_SCAT(Y,STACK(W),NPARPT,
     @                     MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,
     @                     STATUS)
      CALL COPYR8(STACK(W),Y(IPARPT(MYNOD)+1),NPARPT(MYNOD))
      CALL CMPI_RED_SCAT(Z,STACK(W),NPARPT,
     @                     MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,
     @                     STATUS)
      CALL COPYR8(STACK(W),Z(IPARPT(MYNOD)+1),NPARPT(MYNOD))
##ELIF TERRA (Charmm_MPI)
      CALL TCOMMLIBREDUCESCATTER(X,X(IPARPT(MYNOD)+1),NPARPT,
     $     TCOMMLIB_OP_SUM, TCOMMLIB_D_DOUBLE,STATUS)
      CALL TCOMMLIBREDUCESCATTER(Y,Y(IPARPT(MYNOD)+1),NPARPT,
     $     TCOMMLIB_OP_SUM, TCOMMLIB_D_DOUBLE,STATUS)
      CALL TCOMMLIBREDUCESCATTER(Z,Z(IPARPT(MYNOD)+1),NPARPT,
     $     TCOMMLIB_OP_SUM, TCOMMLIB_D_DOUBLE,STATUS)
##ELIF CSPPMPI (Charmm_MPI)
C     Uses non-compliant in-place option
      CALL MPI_REDUCE_SCATTER(X,%VAL(0),NPARPT,
     @                     MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,
     @                     STATUS)
      CALL MPI_REDUCE_SCATTER(Y,%VAL(0),NPARPT,
     @                     MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,
     @                     STATUS)
      CALL MPI_REDUCE_SCATTER(Z,%VAL(0),NPARPT,
     @                     MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,
     @                     STATUS)
##ELSE (Charmm_MPI)
C     Uses compliant calls.
      CALL MPI_REDUCE_SCATTER(X,STACK(W),NPARPT,
     @                     MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,
     @                     STATUS)
      CALL COPYR8(STACK(W),X(IPARPT(MYNOD)+1),NPARPT(MYNOD))
      CALL MPI_REDUCE_SCATTER(Y,STACK(W),NPARPT,
     @                     MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,
     @                     STATUS)
      CALL COPYR8(STACK(W),Y(IPARPT(MYNOD)+1),NPARPT(MYNOD))
      CALL MPI_REDUCE_SCATTER(Z,STACK(W),NPARPT,
     @                     MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,
     @                     STATUS)
      CALL COPYR8(STACK(W),Z(IPARPT(MYNOD)+1),NPARPT(MYNOD))
##ENDIF (Charmm_MPI)
      CALL FRESTK(IREAL8(N))
C
      RETURN
      END
C
      SUBROUTINE VDGBR(X,Y,Z,MODE)
C-----------------------------------------------------------------------
C
C    Wraping routine for vector distributed global broadcast.
C
C
C    MODE=0 - Do just the inner broadcast.
C          (assume values are needed only on processors containing block).
C    MODE=1 - Do the complete global broadcast (inner followed by outer).
C          (assume values are needed everywhere)
C    MODE=2 - Do just the outer broadcast.
C          (assume values are correct for all nodes containing block).
C
C The MODE values 0 and 1 are identical and MODE 2 is ignored for PARAFULL.
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
      REAL*8 X(*), Y(*), Z(*)
      INTEGER MODE
C
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
##INCLUDE '~/charmm_fcm/stack.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
##IF TIMESTAMP
##INCLUDE '~/charmm_fcm/timer.fcm'
##ENDIF
C
##IF CMPI
      INTEGER MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD
      PARAMETER(MPI_DOUBLE_PRECISION=0,MPI_SUM=0,MPI_COMM_WORLD=0)
##ELIF TERRA
      INCLUDE 'tcommlib_f.h'
##ELIF IBMSP
      INCLUDE  'mpif.h'
##ELSE
      INCLUDE 'mpi/mpif.h'
##ENDIF
C
      INTEGER N,W,STATUS
      INTEGER TAG         !##TIMESTAMP
C
C     =====================================
C
##IF PARAFULL
##IF PARASCAL
##ERROR 'Illegal parallel compile options'
##ENDIF
      IF(MODE.EQ.2) RETURN
##ELIF PARASCAL
##ELSE
##ERROR 'Illegal parallel compile options'
##ENDIF
C
##IF TIMESTAMP
      IF (IAND(ISTMPL,2).NE.0) CALL TIME_STAMP(MYNOD,
     $     'VDGBR called from CHARMM'//CHAR(0),IPARPT(NUMNOD))
      IF (ICMRW.EQ.1) THEN
         N=IPARPT(NUMNOD)
         CALL GRDFMM(X,N,'VDGBRX   ')
         CALL GRDFMM(Y,N,'VDGBRY   ')
         CALL GRDFMM(Z,N,'VDGBRZ   ')
         RETURN
      ELSE IF (ICMRW.EQ.2) THEN
         N=IPARPT(NUMNOD)
         CALL WRTFMR(X,'VDGBRX   ',TAG,N)
         CALL WRTFMR(Y,'VDGBRY   ',TAG,N)
         CALL WRTFMR(Z,'VDGBRZ   ',TAG,N)
         RETURN
      ENDIF
##ENDIF
C
      IF(NUMNOD.EQ.1) RETURN
      N=IPARPT(NUMNOD)
      W=ALLSTK(IREAL8(N))
##IF CMPI
C     X is all vector not partial as in standard MPI.
C     In the current implementation NPARPT is ignored
      CALL CMPI_GATHERV(X,NPARPT(MYNOD),
     @                    MPI_DOUBLE_PRECISION,STACK(W),NPARPT,IPARPT,
     @                    MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,STATUS)
      CALL COPYR8(STACK(W),X,N)
      CALL CMPI_GATHERV(Y,NPARPT(MYNOD),
     @                    MPI_DOUBLE_PRECISION,STACK(W),NPARPT,IPARPT,
     @                    MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,STATUS)
      CALL COPYR8(STACK(W),Y,N)
      CALL CMPI_GATHERV(Z,NPARPT(MYNOD),
     @                    MPI_DOUBLE_PRECISION,STACK(W),NPARPT,IPARPT,
     @                    MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,STATUS)
      CALL COPYR8(STACK(W),Z,N)
##ELIF TERRA
      CALL TCOMMLIBALLGATHERV(X(IPARPT(MYNOD)+1),X,8,IPARPT,NPARPT,
     $     STATUS)
      CALL TCOMMLIBALLGATHERV(Y(IPARPT(MYNOD)+1),Y,8,IPARPT,NPARPT,
     $     STATUS)
      CALL TCOMMLIBALLGATHERV(Z(IPARPT(MYNOD)+1),Z,8,IPARPT,NPARPT,
     $     STATUS)
##ELIF CSPPMPI
C     Uses non-compliant inplace option
      CALL MPI_ALLGATHERV(%VAL(0),NPARPT(MYNOD),
     @                    MPI_DOUBLE_PRECISION,X,NPARPT,IPARPT,
     @                    MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,STATUS)
      CALL MPI_ALLGATHERV(%VAL(0),NPARPT(MYNOD),
     @                    MPI_DOUBLE_PRECISION,Y,NPARPT,IPARPT,
     @                    MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,STATUS)
      CALL MPI_ALLGATHERV(%VAL(0),NPARPT(MYNOD),
     @                    MPI_DOUBLE_PRECISION,Z,NPARPT,IPARPT,
     @                    MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,STATUS)
##ELSE
C     Uses compliant calls
      CALL MPI_ALLGATHERV(X(IPARPT(MYNOD)+1),NPARPT(MYNOD),
     @                    MPI_DOUBLE_PRECISION,STACK(W),NPARPT,IPARPT,
     @                    MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,STATUS)
      CALL COPYR8(STACK(W),X,N)
      CALL MPI_ALLGATHERV(Y(IPARPT(MYNOD)+1),NPARPT(MYNOD),
     @                    MPI_DOUBLE_PRECISION,STACK(W),NPARPT,IPARPT,
     @                    MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,STATUS)
      CALL COPYR8(STACK(W),Y,N)
      CALL MPI_ALLGATHERV(Z(IPARPT(MYNOD)+1),NPARPT(MYNOD),
     @                    MPI_DOUBLE_PRECISION,STACK(W),NPARPT,IPARPT,
     @                    MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,STATUS)
      CALL COPYR8(STACK(W),Z,N)
##ENDIF
C
      CALL FRESTK(IREAL8(N))
##IF CONCURR
      CALL PGRPIN(IPMGRP,MYNOD,NUMNOD,INTVCN,INTVMX)
##ENDIF
C
      RETURN
      END
C
      SUBROUTINE VDGBRE(X,KPARPT)
C-----------------------------------------------------------------------
C
C     Global vector distributed broadcast for single vector partitioned
C     in KPARPT array.
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
      REAL*8 X(*)
      INTEGER KPARPT(0:*)
C
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
##INCLUDE '~/charmm_fcm/heap.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
##IF TIMESTAMP
##INCLUDE '~/charmm_fcm/timer.fcm'
##ENDIF
C
##IF CMPI
      INTEGER MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD
      PARAMETER(MPI_DOUBLE_PRECISION=0,MPI_SUM=0,MPI_COMM_WORLD=0)
##ELIF TERRA
      INCLUDE 'tcommlib_f.h'
##ELIF IBMSP
      INCLUDE  'mpif.h'
##ELSE
      INCLUDE 'mpi/mpif.h'
##ENDIF
      INTEGER W,N,STATUS,I,COUNTS(0:MAXNODE)
      INTEGER TAG          !##TIMESTAMP
C
##IF TIMESTAMP
      IF (IAND(ISTMPL,2).NE.0) CALL TIME_STAMP(MYNOD,
     $     'VDGBRE called from CHARMM'//CHAR(0),KPARPT(NUMNOD))
      IF (ICMRW.EQ.1) THEN
         N=KPARPT(NUMNOD)
         CALL GRDFMM(X,N,'VDGBRE   ')
         RETURN
      ELSE IF (ICMRW.EQ.2) THEN
         N=KPARPT(NUMNOD)
         CALL WRTFMR(X,'VDGBRE   ',TAG,N)
         RETURN
      ENDIF
##ENDIF
C
      IF(NUMNOD.EQ.1) RETURN
##IF CMPI
C     For PME performance reason call inplace routine!
      CALL CIMPI_GATHERV(X,MPI_DOUBLE_PRECISION,KPARPT,
     $     MPI_COMM_WORLD,STATUS)
##ELIF TERRA
      DO I = 0, NUMNOD-1
         COUNTS(I)=KPARPT(I+1)-KPARPT(I)
      ENDDO
C
      CALL TCOMMLIBALLGATHERV(X(KPARPT(MYNOD)+1),
     $     X,8,KPARPT,COUNTS,STATUS)
##ELIF CSPPMPI
      DO I = 0, NUMNOD-1
         COUNTS(I)=KPARPT(I+1)-KPARPT(I)
      ENDDO
C
      CALL MPI_ALLGATHERV(%VAL(0),COUNTS(MYNOD),
     @                    MPI_DOUBLE_PRECISION,X,COUNTS,KPARPT,
     @                    MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,STATUS)
##ELSE
      N=KPARPT(NUMNOD)
      W=ALLHP(IREAL8(N))
C
      DO I = 0, NUMNOD-1
         COUNTS(I)=KPARPT(I+1)-KPARPT(I)
      ENDDO
C
      CALL MPI_ALLGATHERV(X(KPARPT(MYNOD)+1),COUNTS(MYNOD),
     @                    MPI_DOUBLE_PRECISION,HEAP(W),COUNTS,KPARPT,
     @                    MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,STATUS)
      CALL COPYR8(HEAP(W),X,N)
      CALL FREHP(W,IREAL8(N))
##ENDIF
C
      RETURN
      END
C
      SUBROUTINE GCOMB(X,N)
C-----------------------------------------------------------------------
C
C     Global COMBine routine.
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
      REAL*8 X(*)
C
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
C
##INCLUDE '~/charmm_fcm/exfunc.fcm'
##INCLUDE '~/charmm_fcm/stack.fcm'
##IF TIMESTAMP
##INCLUDE '~/charmm_fcm/timer.fcm'
##ENDIF
C
##IF CMPI
      INTEGER STATUS,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD
      PARAMETER(MPI_DOUBLE_PRECISION=0,MPI_SUM=0,MPI_COMM_WORLD=0)
##ELIF TERRA
      INCLUDE 'tcommlib_f.h'
      integer counts(0:MAXNODE),rc
##ELIF IBMSP
      INCLUDE  'mpif.h'
##ELSE
      INCLUDE 'mpi/mpif.h'
      INTEGER STATUS
##ENDIF
C
      INTEGER N, W, I
      INTEGER TAG      !##TIMESTAMP
C
##IF TIMESTAMP
      IF (IAND(ISTMPL,2).NE.0) CALL TIME_STAMP(MYNOD,
     $     'GCOMB called from CHARMM'//CHAR(0),N)
      IF (ICMRW.EQ.1) THEN
          CALL GRDFMM(X,N,'GCOMB    ')
          RETURN
      ELSE IF (ICMRW.EQ.2) THEN
          CALL WRTFMR(X,'GCOMB    ',TAG,N)
          RETURN
      ENDIF
##ENDIF
C
      IF(NUMNOD.EQ.1) RETURN
##IF CMPI
      W=ALLSTK(IREAL8(N))
      CALL CIMPI_ALLRED(X,STACK(W),N,MPI_DOUBLE_PRECISION,MPI_SUM,
     @                   MPI_COMM_WORLD,STATUS)
      CALL FRESTK(IREAL8(N))
##ELIF CMPIx
C     Just checking...
      W=ALLSTK(IREAL8(N))
      JPARPT(0)=0
      DO I = 1, NUMNOD
         JPARPT(I)=N*I/NUMNOD
      ENDDO
##IF GENCOMM
      CALL AGVDGS(X,JPARPT,STACK(W),NUMNOD,MYNOD,IPPMAP)
      CALL AGVDGB(X,JPARPT,NUMNOD,MYNOD,IPPMAP)
##ELSE
      CALL AVDGSUM(X,JPARPT,STACK(W),NODDIM,MYNOD,IPPMAP)
      CALL AVDGBR(X,JPARPT,NODDIM,MYNOD,IPPMAP)
##ENDIF
      CALL FRESTK(IREAL8(N))
##ELIF TERRA
      JPARPT(0)=0
      DO I = 1, NUMNOD
         JPARPT(I)=N*I/NUMNOD
      ENDDO
      DO I = 0, NUMNOD-1
         COUNTS(I)=JPARPT(I+1)-JPARPT(I)
      ENDDO
C
      CALL TCOMMLIBREDUCESCATTER(X,X(JPARPT(MYNOD)+1),COUNTS,
     $     TCOMMLIB_OP_SUM, TCOMMLIB_D_DOUBLE,RC)
      CALL TCOMMLIBALLGATHERV(X(JPARPT(MYNOD)+1),X,8,JPARPT,COUNTS,RC)
CC
CC      Check this one instead of above two:
CC      CALL tcommliballreduce(x,x,n,TCOMMLIB_OP_SUM,
CC     $     TCOMMLIB_D_DOUBLE,I)
##ELIF CSPPMPI
C     Uses non-compliant in-place option
      CALL MPI_ALLREDUCE(X,%VAL(0),N,MPI_DOUBLE_PRECISION,MPI_SUM,
     @                   MPI_COMM_WORLD,STATUS)
##ELSE
C     Uses compliant calls.
      W=ALLSTK(IREAL8(N))
      CALL MPI_ALLREDUCE(X,STACK(W),N,MPI_DOUBLE_PRECISION,MPI_SUM,
     @                   MPI_COMM_WORLD,STATUS)
      CALL COPYR8(STACK(W),X,N)
      CALL FRESTK(IREAL8(N))
C
##ENDIF
C
      RETURN
      END
C
      SUBROUTINE IGCOMB(X,N)
C-----------------------------------------------------------------------
C
C     Global COMBine routine. Integer version.
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
      REAL*8 X(*)
C
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
C
##INCLUDE '~/charmm_fcm/exfunc.fcm'
##INCLUDE '~/charmm_fcm/stack.fcm'
##IF TIMESTAMP
##INCLUDE '~/charmm_fcm/timer.fcm'
##ENDIF
C
##IF CMPI
      INTEGER STATUS,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD
      PARAMETER(MPI_INTEGER=1,MPI_SUM=0,MPI_COMM_WORLD=0)
##ELIF TERRA
      INCLUDE 'tcommlib_f.h'
      integer counts(0:MAXNODE),rc
##ELIF IBMSP
      INCLUDE  'mpif.h'
##ELSE
      INCLUDE 'mpi/mpif.h'
      INTEGER STATUS
##ENDIF
C
      INTEGER N, W, I
      INTEGER TAG      !##TIMESTAMP
C
##IF TIMESTAMP
      IF (IAND(ISTMPL,2).NE.0) CALL TIME_STAMP(MYNOD,
     $     'GCOMB called from CHARMM'//CHAR(0),N)
      IF (ICMRW.EQ.1) THEN
          CALL GRDFMM(X,N,'IGCOMB   ')
          RETURN
      ELSE IF (ICMRW.EQ.2) THEN
          CALL WRTFMR(X,'IGCOMB   ',TAG,N)
          RETURN
      ENDIF
##ENDIF
C
      IF(NUMNOD.EQ.1) RETURN
##IF CMPI
      W=ALLSTK(INTEG4(N))
      CALL CIMPI_ALLRED(X,STACK(W),N,MPI_INTEGER,MPI_SUM,
     @                   MPI_COMM_WORLD,STATUS)
      CALL FRESTK(INTEG4(N))
##ELIF TERRA
      JPARPT(0)=0
      DO I = 1, NUMNOD
         JPARPT(I)=N*I/NUMNOD
      ENDDO
      DO I = 0, NUMNOD-1
         COUNTS(I)=JPARPT(I+1)-JPARPT(I)
      ENDDO
C
      CALL TCOMMLIBREDUCESCATTER(X,X(JPARPT(MYNOD)+1),COUNTS,
     $     TCOMMLIB_OP_SUM, TCOMMLIB_D_UINT,RC)
      CALL TCOMMLIBALLGATHERV(X(JPARPT(MYNOD)+1),X,4,JPARPT,COUNTS,RC)
CC
CC      Check this one instead of above two:
CC      CALL tcommliballreduce(x,x,n,TCOMMLIB_OP_SUM,
CC     $     TCOMMLIB_D_DOUBLE,I)
##ELIF CSPPMPI
C     Uses non-compliant in-place option
      CALL MPI_ALLREDUCE(X,%VAL(0),N,MPI_INTEGER,MPI_SUM,
     @                   MPI_COMM_WORLD,STATUS)
##ELSE
C     Uses compliant calls.
      W=ALLSTK(INTEG4(N))
      CALL MPI_ALLREDUCE(X,STACK(W),N,MPI_INTEGER,MPI_SUM,
     @                   MPI_COMM_WORLD,STATUS)
      CALL COPYR8(STACK(W),X,N)
      CALL FRESTK(INTEG4(N))
C
##ENDIF
C
      RETURN
      END
C
      SUBROUTINE GBOR(X,N)
C-----------------------------------------------------------------------
C
C     Global Bitwise OR routine.
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
      INTEGER X(*), N
C
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
C
##INCLUDE '~/charmm_fcm/exfunc.fcm'
##INCLUDE '~/charmm_fcm/stack.fcm'
##IF TIMESTAMP
##INCLUDE '~/charmm_fcm/timer.fcm'
##ENDIF
C
##IF CMPI
      INTEGER MPI_INTEGER,MPI_BOR,MPI_COMM_WORLD
      PARAMETER(MPI_INTEGER=1,MPI_BOR=1,MPI_COMM_WORLD=0)
##ELIF TERRA
      INCLUDE 'tcommlib_f.h'
##ELIF IBMSP
      INCLUDE  'mpif.h'
##ELSE
      INCLUDE 'mpi/mpif.h'
##ENDIF
C
      INTEGER W, STATUS,I
      INTEGER TAG          !##TIMESTAMP
C
##IF TIMESTAMP
      IF (IAND(ISTMPL,2).NE.0) CALL TIME_STAMP(MYNOD,
     $     'GBOR called from CHARMM'//CHAR(0),N)
CCC      DO I=1,N,5
CCC         WRITE(*,1)MYNOD,I,I+4,X(I),X(I+1),X(I+2),X(I+3),X(I+4)
CCC      ENDDO
CCC 1    FORMAT('GBOR>ME=',I3,' X(',I5,'-',I5,')=',5Z10)
      IF (ICMRW.EQ.1) THEN
          CALL GIDFMM(X,N,'GBOR     ')
          RETURN
      ELSE IF (ICMRW.EQ.2) THEN
          CALL WRTFMI(X,'GBOR     ',TAG,N)
          RETURN
      ENDIF
C      
##ENDIF
C
      IF(NUMNOD.EQ.1) RETURN
      W=ALLSTK(INTEG4(N))
##IF CMPI
C     Call inplace routine
      CALL CIMPI_ALLRED(X,STACK(W),N,MPI_INTEGER,MPI_BOR,
     $     MPI_COMM_WORLD,STATUS)
##ELIF TERRA
      CALL TCOMMLIBALLREDUCE(X,X,N,TCOMMLIB_OP_BOR,
     $     TCOMMLIB_D_UINT,STATUS)
##ELIF CSPPMPI
C     Uses non-compliant in-place option
      CALL MPI_ALLREDUCE(X,%VAL(0),N,MPI_INTEGER,MPI_BOR,
     @                   MPI_COMM_WORLD,STATUS)
##ELSE
C     Uses compliant calls.
      CALL MPI_ALLREDUCE(X,STACK(W),N,MPI_INTEGER,MPI_BOR,
     @                   MPI_COMM_WORLD,STATUS)
      CALL COPYI4(STACK(W),X,N)
##ENDIF
C
      CALL FRESTK(INTEG4(N))
C
      RETURN
      END
C
      SUBROUTINE PSYNC()
C-----------------------------------------------------------------------
C
C     This is a wrapper routine for global sync.
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
##IF CMPI
      INTEGER MPI_COMM_WORLD
      PARAMETER(MPI_COMM_WORLD=0)
##ELIF IBMSP
      INCLUDE  'mpif.h'
##ELSE
      INCLUDE 'mpi/mpif.h'
##ENDIF
      INTEGER STATUS
C
      IF(NUMNOD.EQ.1) RETURN
      IF(NUMNOD.NE.1) RETURN !##LOBOS
##IF CMPI
      CALL CMPI_BARRIER(MPI_COMM_WORLD,STATUS)
##ELIF TERRA
      CALL TCOMMLIBBARRIER(STATUS)
##ELSE
      CALL MPI_BARRIER(MPI_COMM_WORLD,STATUS)
##ENDIF
C
      RETURN
      END
C
      SUBROUTINE SWAPD(N,X,Y,Z,XPL,YPL,ZPL,XPR,YPR,ZPR,LPEER,RPEER)
C-----------------------------------------------------------------------
C This routine performs a double swap of data between
C {X,Y,Z} and {XPL,YPL,ZPL}, and at the same time also {X,Y,Z} and {XPR,YPR,ZPR}.
C X,Y,Z is on MYNOD, XPL,YPL,ZPL is on LPEER,XPR,YPR,ZPR, is on RPEER.
C If RPEER=-1 or LPEER=-1, then skip the swap in this direction, ie call GRECSEN
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
      INTEGER N,LPEER,RPEER
      REAL*8 X(*),Y(*),Z(*),XPL(*),YPL(*),ZPL(*),XPR(*),YPR(*),ZPR(*)
C
      CALL SWAPD1(N,X,XPL,XPR,LPEER,RPEER)
      CALL SWAPD1(N,Y,YPL,YPR,LPEER,RPEER)
      CALL SWAPD1(N,Z,ZPL,ZPR,LPEER,RPEER)
C
      RETURN
      END
C
      SUBROUTINE PSNDC(ARRAY, LENGTH)
C-----------------------------------------------------------------------
C This routine performs on node 0 broadcast to all other nodes
C and receive from node 0 on all other nodes.
C Usually called after read on node 0. For character arrays.
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
      CHARACTER*(*) ARRAY(*)
      INTEGER LENGTH
C
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
##IF CMPI
      INTEGER MPI_COMM_WORLD,MPI_BYTE
      PARAMETER(MPI_COMM_WORLD=0,MPI_BYTE=2)
##ELIF IBMSP
      INCLUDE  'mpif.h'
##ELSE
      INCLUDE 'mpi/mpif.h'
##ENDIF
      INTEGER STATUS,LENST
C
      IF((NUMNOD.EQ.1).OR.(LENGTH.EQ.0)) RETURN
      LENST = LEN(ARRAY(1))
##IF T3D T3E
      CALL CT3C_BCAST(ARRAY,LENST,LENGTH)
##ELIF CMPI
      CALL CMPI_BCAST(ARRAY,LENST*LENGTH,MPI_BYTE,0,MPI_COMM_WORLD,
     $     STATUS)
##ELIF TERRA
      CALL TCOMMLIBBROADCAST(0,ARRAY,ARRAY,LENST*LENGTH,STATUS)
##ELSE
      CALL MPI_BCAST(ARRAY,LENST*LENGTH,MPI_BYTE,0,MPI_COMM_WORLD,
     $     STATUS)
##ENDIF
C
      RETURN
      END
C
      SUBROUTINE PSND4(ARRAY, LENGTH)
C-----------------------------------------------------------------------
C This routine performs on node 0 broadcast to all other nodes
C and receive from node 0 on all other nodes.
C Usually called after read on node 0. For single precision arrays.
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
      INTEGER ARRAY(*)
      INTEGER LENGTH
C
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
C
##IF T3D T3E
      integer offset, log2npes,type
##ENDIF
##IF CMPI
      INTEGER MPI_COMM_WORLD,MPI_BYTE
      PARAMETER(MPI_COMM_WORLD=0,MPI_BYTE=2)
##ELIF IBMSP
      INCLUDE  'mpif.h'
##ELSE
      INCLUDE 'mpi/mpif.h'          !##.not.TERRA
##ENDIF
      INTEGER STATUS
C
      IF(NUMNOD.EQ.1) RETURN
##IF T3D T3E
      call barrier()
      type=1
      offset = 1
      do i = 1, LOG2NPES()
        if( mynod .lt. offset ) then
            call gsen( mynod+offset, type, array, 8*length )
            call barrier()
        else if( mynod .lt. offset+offset ) then
            call grec( mynod-offset, type, array, 8*length )
            call barrier()
        else
            call barrier()
        end if
        offset = offset + offset
      end do
CC
CC       Later improvement: fixed the above in the below
CCC      CALL CT3I_BCAST(ARRAY,LENGTH)
##ELIF CMPI
##IF INTEGER8
C for case where integer is 8 bytes 
      CALL CMPI_BCAST(ARRAY,8*LENGTH,MPI_BYTE,0,MPI_COMM_WORLD,STATUS)
##ELSE
      CALL CMPI_BCAST(ARRAY,4*LENGTH,MPI_BYTE,0,MPI_COMM_WORLD,STATUS)
##ENDIF
##ELIF TERRA
      CALL TCOMMLIBBROADCAST(0,ARRAY,ARRAY,4*LENGTH,STATUS)
##ELSE
      CALL MPI_BCAST(ARRAY,4*LENGTH,MPI_BYTE,0,MPI_COMM_WORLD,STATUS)
##ENDIF
      RETURN
      END
C
      SUBROUTINE PSND8(ARRAY, LENGTH)
C-----------------------------------------------------------------------
C This routine performs on node 0 broadcast to all other nodes
C and receive from node 0 on all other nodes.
C Usually called after read on node 0. For Real*8 precision arrays.
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
      REAL*8 ARRAY(*)
      INTEGER LENGTH
C
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
C
##IF CMPI
      INTEGER MPI_COMM_WORLD,MPI_BYTE
      PARAMETER(MPI_COMM_WORLD=0,MPI_BYTE=2)
##ELIF IBMSP
      INCLUDE  'mpif.h'
##ELSE
      INCLUDE 'mpi/mpif.h'
##ENDIF
      INTEGER STATUS
C
      IF(NUMNOD.EQ.1) RETURN
##IF CMPI
      CALL CMPI_BCAST(ARRAY,8*LENGTH,MPI_BYTE,0,MPI_COMM_WORLD,STATUS)
##ELIF TERRA
      CALL TCOMMLIBBROADCAST(0,ARRAY,ARRAY,8*LENGTH,STATUS)
##ELSE
      CALL MPI_BCAST(ARRAY,8*LENGTH,MPI_BYTE,0,MPI_COMM_WORLD,STATUS)
##ENDIF
      RETURN
      END
C
      SUBROUTINE BALPRN()
C-----------------------------------------------------------------------
C
C     This routine prints load balance information at the end
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
##INCLUDE '~/charmm_fcm/number.fcm'
##INCLUDE '~/charmm_fcm/stream.fcm'
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
C
      INTEGER I,J
      REAL*8 TSUMT
      REAL*8 TIMONE(MAXNODE),TIMALL(MAXNODE,TMGNUM)
C
C     Subtract external energy times from total sum
c      TMERI(TIMINTE)=TMERI(TIMINTE)-TMERI(TIMEXTE)
      IF(PRNLEV.GE.2) THEN
         WRITE(OUTU,*)'Parallel load balance (sec.):'
         WRITE(OUTU,*)
     $     'Node Eext      Eint   Wait    Comm    List   Integ   Total'
         TSUMT=ZERO
         DO I = 1, TMGNUM
            TSUMT=TSUMT+TMERI(I)
         ENDDO
         WRITE(OUTU,1510)0, (TMERI(I), I=1,TMGNUM),TSUMT
      ENDIF
C
      DO I = 1, TMGNUM
         TSUMI(I) = TMERI(I)
      ENDDO
C
      IF (NUMNOD.GT.1) THEN
C
         DO I = 1, TMGNUM
Cmh..09-FEB-98 fix for T3E port
            DO J = 1, NUMNOD
               TIMONE(J) = ZERO
            ENDDO
            TIMONE(MYNODP) = TMERI(I)
            CALL GCOMB(TIMONE,NUMNOD)
            DO J = 2, NUMNOD
               TIMALL(J,I) = TIMONE(J)
            ENDDO
         ENDDO
         DO I = 2, NUMNOD
            TSUMT=ZERO
            DO J = 1, TMGNUM
               TSUMT=TSUMT+TIMALL(I,J)
            ENDDO
            IF(PRNLEV.GE.2) WRITE(OUTU,1510) I-1,
     $           (TIMALL(I,J), J=1, TMGNUM),TSUMT
            DO J = 1, TMGNUM
               TSUMI(J) = TSUMI(J) + TIMALL(I,J)
            ENDDO
         ENDDO
         IF(PRNLEV.GE.2) THEN
            DO I = 1, TMGNUM
               TSUMI(I) = TSUMI(I)/NUMNOD
            ENDDO
            TSUMT=ZERO
            DO I = 1, TMGNUM
               TSUMT=TSUMT+TSUMI(I)
            ENDDO
            WRITE(OUTU,'(A)') 'PARALLEL> Average timing for all nodes:'
            WRITE(OUTU,1510)
     &           NUMNOD,(TSUMI(I), I = 1, TMGNUM),TSUMT
         ENDIF
      ENDIF
 1510 FORMAT(I4,13F8.1)
C
      RETURN
      END
C

##IF TIMESTAMP (test_nocommunication)
      SUBROUTINE TSTNCM(COMLYN,COMLEN)
C-----------------------------------------------------------------------
C
C     Turns off the interprocessor communications in parallel.
C     Syntax:
C
C     TEST NOCOmmunication READ|WRITE UNIT <int> STEP <int> MEMO <int>
C     TEST CLOSe
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
##INCLUDE '~/charmm_fcm/dimens.fcm'
##INCLUDE '~/charmm_fcm/psf.fcm'
##INCLUDE '~/charmm_fcm/heap.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
##INCLUDE '~/charmm_fcm/timer.fcm'
C
CCC      INCLUDE 'mpif.h'
      INTEGER STATUS
C
C . Passed variables.
      CHARACTER COMLYN*(*)
      INTEGER   COMLEN
C . Local variables.
      CHARACTER*4  WRD
      INTEGER LENMEM,CTAG,CSIZE,CDAT
C
      WRD=NEXTA4(COMLYN,COMLEN)
C . Branch on the option.
C=======================================================================
      IF(WRD.EQ.'READ') THEN
         ICMUN=GTRMI(COMLYN,COMLEN,'UNIT',17)
         ICMST=GTRMI(COMLYN,COMLEN,'STEP',1)
         ICMRW=1
         ICMHPX=GTRMI(COMLYN,COMLEN,'MEMO',0)
         CALL PSND4(ICMST,1)
         CALL PSND4(ICMHPX,1)
         CALL PSND4(ICMRW,1)
         IF(ICMHPX.LT.NATOM*6*ICMST) THEN
            CALL WRNDIE(-5,'<TSTNCM>','Not enough memory.')
         ENDIF
C        Allocate enough memory for all steps (coor + force + fft)
         ICMHP=ALLHP(INTEG4(ICMHPX))
         ICMHPO=0
C        Read the data from the file into HEAP
         IF(MYNOD.EQ.0) CALL RFPDIM(LENMEM)
         ICMHPO=0
         ICMHPL=LENMEM
         CALL PSND4(ICMHPL,1)
C        Broadcast to the net
CC         write(*,*)'TEST NOCOM> About to send HEAP of length =',ICMHPL
         CALL PSND4(HEAP(ICMHP),ICMHPL)
CC         write(*,*)'TEST NOCOM> PSND4 OK.'
C
C        call barrier to synchronize everybody!
C
CCC         CALL MPI_BARRIER(MPI_COMM_WORLD,STATUS)
         CALL PSYNC()
      ELSE IF(WRD.EQ.'WRIT') THEN
         ICMUN=GTRMI(COMLYN,COMLEN,'UNIT',17)
         ICMST=GTRMI(COMLYN,COMLEN,'STEP',1)
         ICMHPX=GTRMI(COMLYN,COMLEN,'MEMO',0)
         ICMRW=2
         IF(NUMNOD.GT.1) CALL WRNDIE(-5,'<TSTNCM>',
     $        'TEST NOCO WRITE must be performed on a single node.')
      ELSE IF(WRD.EQ.'CLOS') THEN
CCC         CALL MPI_BARRIER(MPI_COMM_WORLD,STATUS)
         CALL PSYNC()
         IF(ICMRW.EQ.1) CALL FREHP(ICMHP,ICMHPL)
         IF(ICMRW.EQ.2) THEN
            write(*,*)'TEST NOCOM> Closing the file.'
            CTAG=0 
            CSIZE=0
            CDAT=0
            CALL WRTFMI(CDAT,'ENDOFFILE',CTAG,CSIZE)
         ENDIF
C        Disable further simulation of the communication:
         ICMRW=0
      ELSE
         CALL WRNDIE(-5,'<TSTNCM>',
     $        'TEST NOCOmmunication needs READ, WRITE, CLOSE.')
      ENDIF
      RETURN
      END
C
      SUBROUTINE WRTFMR(DATA,NAME,TAG,SIZE)
C-----------------------------------------------------------------------
C
C     [WRTFMR = WRiTe From Memory Real*8]
C
C This routine write the data array into the file and
C prepends the header to it so later RDMEM can put data
C in the right place into a HEAP storage
C
      REAL*8 DATA(*)
      CHARACTER*9 NAME
      INTEGER TAG,SIZE
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
##INCLUDE '~/charmm_fcm/timer.fcm'
C
      INTEGER I
C
      TAG = ICMHPO
      write(*,*)'TEST NOCOM WRITE>name,tag,size/real=',name,tag,size
      WRITE(ICMUN)NAME,'R8',TAG,SIZE
      WRITE(ICMUN)(DATA(I),I=1,SIZE)
      ICMHPO=TAG+IREAL8(SIZE)
C
      RETURN
      END
C
      SUBROUTINE WRTFMI(DATA,NAME,TAG,SIZE)
C-----------------------------------------------------------------------
C
C     [WRTFMR = WRiTe From Memory Integer]
C
C This routine write the data array into the file and
C prepends the header to it so later RDMEM can put data
C in the right place into a HEAP storage
C
      INTEGER DATA(*)
      CHARACTER*9 NAME
      INTEGER TAG,SIZE
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
##INCLUDE '~/charmm_fcm/timer.fcm'
C
      INTEGER I
C
      TAG = ICMHPO
      write(*,*)'TEST NOCOM WRITE>name,tag,size/int=',name,tag,size
      WRITE(ICMUN)NAME,'I4',TAG,SIZE
      IF(SIZE.GT.0)WRITE(ICMUN)(DATA(I),I=1,SIZE)
      ICMHPO=TAG+INTEG4(SIZE)
      IF(NAME.EQ.'ENDOFFILE') THEN
         ICMHPO=0
         CLOSE(ICMUN)
      ENDIF
C
      RETURN
      END
C
      SUBROUTINE RFPDIM(LENMEM)
C-----------------------------------------------------------------------
C
C     [RFPDIM = Read File and Put Data Into Memory]
C
C This routine reads the entire file made by WRTFM{I,R} and stores 
C all the data contiguously into the HEAP. Pointers are in parallel.fcm
C
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
      INTEGER LENMEM
C
##INCLUDE '~/charmm_fcm/exfunc.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
##INCLUDE '~/charmm_fcm/heap.fcm'
##INCLUDE '~/charmm_fcm/timer.fcm'
C
      CHARACTER*9 NAME
      CHARACTER*2 TYPE
      INTEGER TAG,SIZE,I,N
C
      ICMHPO=0
C
 100  CONTINUE
      READ(ICMUN)NAME,TYPE,TAG,SIZE
      write(*,*)'TEST NOCOM READ>name,type,tag,size=',
     $     name,type,tag,size
C
      IF(NAME.EQ.'ENDOFFILE') GOTO 999
      IF(TYPE.EQ.'R8') THEN
         N=IREAL8(SIZE)
         ICMHPO=ICMHPO+N
      ELSE IF(TYPE.EQ.'I4') THEN
         N=INTEG4(SIZE)
         ICMHPO=ICMHPO+N
      ELSE
         CALL WRNDIE(-5,'<RFPDIM>','Error in type spec.')
      ENDIF
      IF(ICMHPO.GT.ICMHPX) THEN
         WRITE(*,*)'TEST NOCOM> Max memory specified =',ICMHPX
         WRITE(*,*)'TEST NOCOM> Memory used =',ICMHPO
         CALL WRNDIE(-5,'<TSTNCM>','Not enough memory.')
      ENDIF
      READ(ICMUN)(HEAP(ICMHP+ICMHPO-N+I-1),I=1,N)
      GOTO 100
 999  CONTINUE
C     Print some statistics:
C
      write(*,*)'TEST NOCOM> Memory used =',ICMHPO
C
      LENMEM=ICMHPO
C
      RETURN
      END
C
      SUBROUTINE GRDFMM(X,SIZE,NAME)
C-----------------------------------------------------------------------
C
C     [GRDFMM = Get Real*8 Data From MeMory]
C
C This routine puts the data from contiguous HEAP storage
C into a data array. Real case. Offset is stored in parallel.fcm
C It replaces getting data over the network - instead we get
C it from the memory.
C
      REAL*8 X(*)
      INTEGER SIZE,I
      CHARACTER*9 NAME
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
##INCLUDE '~/charmm_fcm/heap.fcm'
##INCLUDE '~/charmm_fcm/timer.fcm'
C
CC      write(*,*)'TEST NOCOM GRDFMM>me,name,offset,size=',
CC     $     mynod,name,icmhpo,size
      CALL COPYR8(HEAP(ICMHP+ICMHPO),X,SIZE)
      ICMHPO=ICMHPO+IREAL8(SIZE)
C
      RETURN
      END
C
      SUBROUTINE GIDFMM(X,SIZE,NAME)
C-----------------------------------------------------------------------
C
C     [GIDFMM = Get Integer Data Form MeMory]
C
C This routine puts the data from contiguous HEAP storage
C into a data array. Integer case. Offset is stored in COMMON
C It replaces getting data over the network - instead we get
C it from the memory.
C
      INTEGER X(*)
      INTEGER SIZE
      CHARACTER*9 NAME
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/exfunc.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
##INCLUDE '~/charmm_fcm/heap.fcm'
##INCLUDE '~/charmm_fcm/timer.fcm'
C
CC      write(*,*)'TEST NOCOM GIDFMM>me,name,offset,size=',
CC     $     mynod,name,icmhpo,size
      CALL COPYI4(HEAP(ICMHP+ICMHPO),X,SIZE)
      ICMHPO=ICMHPO+INTEG4(SIZE)
C
      RETURN
      END
C
##ENDIF (test_nocommunication)
C
      SUBROUTINE PARFIN
C-----------------------------------------------------------------------
C
C Finishing routine - required for some platforms or libraries
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
##INCLUDE '~/charmm_fcm/parallel.fcm'
##INCLUDE '~/charmm_fcm/heap.fcm'
##IF TIMESTAMP
##INCLUDE '~/charmm_fcm/timer.fcm'
##ENDIF
C
      INTEGER STATUS     !##.NOT.INTEGER8
      INTEGER*4 STATUS   !##INTEGER8
##IF GAMESSUK
      INTEGER ISTATUS
##ENDIF
##IF SOCKET
      INTEGER J,LENH,LUNASS
      EXTERNAL LUNASS
      CHARACTER*80 ENV0
##ENDIF
C
##IF CMPI
##IF MPI
##IF TIMESTAMP
      IF(ISTMPL.NE.0) THEN
         call time_stamp(mynod,
     $        'Finished CHARMM time stamp'//char(0),NUMNOD)
         call flush_time_stamp(mynod,numnod)
      ENDIF
##ENDIF

##IF GAMESSUK
      CALL PARALLEL_END(ISTATUS)
##ELSE
      CALL MPI_FINALIZE(STATUS)
##ENDIF
##ELSE
      STATUS=0
##ENDIF
##ELSE
      CALL MPI_FINALIZE(STATUS)
##ENDIF
##IF SOCKET
C LNI March 1998.
C Cleanup
C  It seems that this job has finished OK. The
C  script for slave removal will not be needed so we delete it
C  unless its name starts with 'KEEP' (for debugging)
C
      IF(MYNOD.EQ.0 .AND. NUMNOD.GT.1)THEN
         CALL FGETENV(ENV0,LENH,'REMOVE_CH_SLAVES'//CHAR(0))
         IF(LENH.LE.0) RETURN
         IF(ENV0(1:4).EQ.'KEEP') RETURN
         J=LUNASS(90)
         OPEN(J,NAME=ENV0(1:LENH),ERR=999,STATUS='UNKNOWN')
         CLOSE(J,STATUS='DELETE',ERR=999)
      ENDIF
999   RETURN
##ELSE
      RETURN
##ENDIF
      END
C
      REAL*8 FUNCTION ECLOCK()
C-----------------------------------------------------------------------
C
C Timing routine which ususally returns elapsed time
C
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
##IF IBMSP1
Cln...23-Jul-96
      REAL*8 RTC
      ECLOCK=RTC()
C
##ELIF T3D
      ECLOCK=RTC() * 6.67 E-9
##ELIF T3E
Clb...ECLOCK=RTC() * 3.335 E-9
      integer get_clk
      ECLOCK=RTC() * 1.0d0 / ( get_clk() )
##ELIF CM5
      INCLUDE '/usr/include/cm/cmmd_fort.h'
      INTEGER ITMP
      ITMP=CMMD_NODE_TIMER_STOP(1)
      ECLOCK=CMMD_NODE_TIMER_BUSY(1)
      ITMP=CMMD_NODE_TIMER_START(1)
C
##ELIF TERRA
CC      REAL*8 Ttime_Wctime
      INCLUDE 'Ttimelib_f.h'
      real timeval(2),tt
      tt = ttime_processtime(timeval)
      eclock = timeval(1) + timeval(2)
##ELIF ULTRA
      REAL*8 TIM
      real tim2(2),dtime
      TIM=dtime(tim2)
      ECLOCK=TIM
##ELSE
      REAL*8 TIM,DCLOCK
      TIM=DCLOCK()
      ECLOCK=TIM
##ENDIF
C
      RETURN
      END
C
      SUBROUTINE CUBE(MYNOD,NUMNOD,IPPMAP)
C-----------------------------------------------------------------------
##INCLUDE '~/charmm_fcm/impnon.fcm'
C
C     Defines IPPMAP for standard cube network node numbers
C
      INTEGER MYNOD, NUMNOD, IPPMAP(0:*)
      INTEGER I
C
      DO I = 0, NUMNOD
         IPPMAP(I)=I
      ENDDO
      RETURN
##ENDIF ! PARALLEL
      END
