
	program NBODY

	include '../include/express.h'

	include '../include/nbody.h'

	include '../include/sysdep.h'

	data 	allNode,npid,gtype,stype,intSize,realSize,cmplxSize 
     &			/ -1,0,999,0,4,4,8 /

	integer timeBegin, timeEnd, mclock
	real dx(0:NN-1), dy(0:NN-1), dz(0:NN-1), dist(0:NN-1), sq(0:NN-1)
        real fac(0:NN-1), tx(0:NN-1), ty(0:NN-1), tz(0:NN-1)
	real p(0:NN-1,0:6), q(0:NN-1,0:6) 
	real buf1(0:NN*7), buf2(0:NN*7) 
	integer i,j,k,npts,dest,t,tp

        integer status,src,dummy,dest

        CALL KXINIT
        status = KXPARA(ENV)

        thisNode = ENV(1)
        hostNode = IHOST

c	the total number of points is received by each node.
        status = KXBROD(npts, IHOST, intSize, IALNOD, dummy, gtype)

	numNode	= ENV(2)
	numDim	= nodedim()

c       numLoc gives the number of points allocated to each node.
c       The optimum decomposition is such that equal number of particles are
c       stored in every node.
	numLoc  = npts/numNode

	src = NOCARE
	status = KXREAD(buf1, numLoc*7*realSize, src, gtype)

	minLoc	= thisNode*numLoc
	maxLoc 	= minLoc + numLoc - 1
	if (maxLoc .GT. npts-1) then
	  maxLoc = npts-1
	  numLoc = maxLoc - minLoc + 1
	endif
	do i = 0, numLoc-1
	  do j = 0, 6
  	    p(i,j) = buf1(i*7+j)
	  end do
	end do

        timeBegin = mclock()
	
        do i = 0, numLoc-1
          do j = 0, 6
            q(i,j) = p(i,j)
          end do
        end do

        k = npts/2 - 1

        l = k / numLoc
        if (MOD(k, numLoc) .NE. 0)  l = l + 1

	do i = 0, numLoc-1
	  do j = i+1, numLoc-1
  	    dx(i) = p(i,PX) - q(j,PX)
	    dy(i) = p(i,PY) - q(j,PY)
	    dz(i) = p(i,PZ) - q(j,PZ)
	    sq(i) = dx(i)**2+dy(i)**2+dz(i)**2
	    dist(i) = SQRT(sq(i))
	    fac(i) = p(i,MM) * q(j,MM) / (dist(i) * sq(i))
	    tx(i) = fac(i) * dx(i)
	    ty(i) = fac(i) * dy(i)
	    tz(i) = fac(i) * dz(i)
	    p(i,FX) = p(i,FX)-tx(i)
	    q(j,FX) = q(j,FX)+tx(i)
	    p(i,FY) = p(i,FY)-ty(i)
	    q(j,FY) = q(j,FY)+ty(i)
	    p(i,FZ) = p(i,FZ)-tz(i)
	    q(j,FZ) = q(j,FZ)+tz(i)
          end do
        end do

        do i = 1, l - 1
          call CSHIFT(q, 1, numLoc, numLoc, 7, buf1, buf2, i)
          t = MIN(i+k-i*numLoc, numLoc-1)
          do j = 0, numLoc-1
            do j1 = 0, t
  	      dx(j) = p(j,PX) - q(j1,PX)
	      dy(j) = p(j,PY) - q(j1,PY)
	      dz(j) = p(j,PZ) - q(j1,PZ)
	      sq(j) = dx(j)**2+dy(j)**2+dz(j)**2
	      dist(j) = SQRT(sq(j))
	      fac(j) = p(j,MM) * q(j1,MM) / (dist(j) * sq(j))
	      tx(j) = fac(j) * dx(j)
	      ty(j) = fac(j) * dy(j)
	      tz(j) = fac(j) * dz(j)
	      p(j,FX) = p(j,FX)-tx(j)
	      q(j1,FX) = q(j1,FX)+tx(j)
	      p(j,FY) = p(j,FY)-ty(j)
	      q(j1,FY) = q(j1,FY)+ty(j)
	      p(j,FZ) = p(j,FZ)-tz(j)
	      q(j1,FZ) = q(j1,FZ)+tz(j)
            end do
          end do
        end do

	if (numNode .GT. 1) then
          tp = MOD(k, numLoc)
          if (tp .NE. 0) then
            call CSHIFT(q, 1, tp, numLoc, 7, buf1, buf2, i)
            do j = numLoc-tp, numLoc-1
              do j1 = numLoc-tp, j
      	        dx(j) = p(j,PX) - q(j1,PX)
	        dy(j) = p(j,PY) - q(j1,PY)
	        dz(j) = p(j,PZ) - q(j1,PZ)
  	        sq(j) = dx(j)**2+dy(j)**2+dz(j)**2
	        dist(j) = SQRT(sq(j))
	        fac(j) = p(j,MM) * q(j1,MM) / (dist(j) * sq(j))
	        tx(j) = fac(j) * dx(j)
	        ty(j) = fac(j) * dy(j)
	        tz(j) = fac(j) * dz(j)
	        p(j,FX) = p(j,FX)-tx(j)
	        q(j1,FX) = q(j1,FX)+tx(j)
	        p(j,FY) = p(j,FY)-ty(j)
	        q(j1,FY) = q(j1,FY)+ty(j)
	        p(j,FZ) = p(j,FZ)-tz(j)
	        q(j1,FZ) = q(j1,FZ)+tz(j)
              end do
            end do
          end if

          call CSHIFT(q, 1, 1, numLoc, 7, buf1, buf2, i+1)
          tp = MIN(k-minLoc, numLoc-1)
          do j = 0, tp
            dx(j) = p(j,PX) - q(j,PX)
	    dy(j) = p(j,PY) - q(j,PY)
	    dz(j) = p(j,PZ) - q(j,PZ)
  	    sq(j) = dx(j)**2+dy(j)**2+dz(j)**2
	    dist(j) = SQRT(sq(j))
	    fac(j) = p(j,MM) * q(j,MM) / (dist(j) * sq(j))
	    tx(j) = fac(j) * dx(j)
	    ty(j) = fac(j) * dy(j)
	    tz(j) = fac(j) * dz(j)
	    p(j,FX) = p(j,FX)-tx(j)
	    q(j,FX) = q(j,FX)+tx(j)
	    p(j,FY) = p(j,FY)-ty(j)
	    q(j,FY) = q(j,FY)+ty(j)
	    p(j,FZ) = p(j,FZ)-tz(j)
	    q(j,FZ) = q(j,FZ)+tz(j)
          end do

          call CSHIFT(q, 1, npts/2, numLoc, 7, buf1, buf2, i+2)
        end if
        do j = 0, numLoc-1
	  p(FX,i) = p(FX,i) + q(FX,i)
	  p(FY,i) = p(FY,i) + q(FY,i)
	  p(FZ,i) = p(FZ,i) + q(FZ,i)
	end do

        timeEnd = mclock()
        print *,'Node', thisNode,'  Elapsed time:',
     &          timeEnd-timeBegin,'mS'

        k = 0
        do i = 0, numLoc-1
          do j = 0, 6
            buf1(k) = p(i,j)
            k = k + 1
          end do
        end do
        status = KXWRIT(buf1,realSize*7*npts/numNode,IHOST,gtype)

	stop
	end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C									C
C			CSHIFT						C
C									C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

	subroutine CSHIFT(a, dim, sh, n1, n2, buf1, buf2, indx)
	integer dim, sh, n1, n2, indx

	include '../include/express.h'

	include '../include/nbody.h'

	real a(0:NN-1,0:n2-1), buf1(0:n1*n2-1), buf2(0:n1*n2-1)
        
	include '../include/sysdep.h'


	integer i, j, k, desPE1, desPE2, srcPE1, srcPE2, offset1, offset2

C       send, DIM=1
	if ((maxLoc-sh) .LT. 0) then
	    desPE1 = (maxLoc-sh+npts)/n1
	else
	    desPE1 = (maxLoc-sh)/n1
	endif
        if (desPE1 .NE. thisNode) then
	  k = 0
	  do i = 0, n1-1
	    do j = 0, n2-1
	      buf1(k) = a(i,j)
	      k = k + 1
	    end do
	  end do
	 status= KXWRIT(buf1,realSize*n1*n2,desPE1,stype+indx*2)
        endif

	if ((minLoc-sh) .LT. 0) then
	    desPE2 = (minLoc-sh+npts)/n1
	else
            desPE2 = (minLoc-sh)/n1
	endif
        if ((desPE2 .NE. thisNode) .AND. (desPE2 .NE. desPE1)) then
	  k = 0
	  do i = 0, n1-1
	    do j = 0, n2-1
	      buf2(k) = a(i,j)
	      k = k + 1
	    end do
	  end do
	 status = KXWRIT(buf2,realSize*n1*n2,desPE2,stype+indx*2+1)

        endif

C       receive
	if ((minLoc+sh) .GE. npts) then
	    srcPE1 = (minLoc+sh-npts)/n1
	else
	    srcPE1 = (minLoc+sh)/n1
	endif
        if (srcPE1 .NE. thisNode) then
	  src = NOCARE
	  status = KXREAD(buf1,realSize*n1*n2,src,stype+indx*2)
          offset1 = (srcPE1-thisNode)*n1
	  if (offset1 .LT. 0) then
	    offset1 = offset1 + npts
	  endif
        endif

	if ((maxLoc+sh) .GE. npts) then
	    srcPE2 = (maxLoc+sh-npts)/n1
	else
            srcPE2 = (maxLoc+sh)/n1
	endif
        if ((srcPE2 .NE. thisNode) .AND. (srcPE2 .NE. srcPE1)) then
	  src = NOCARE
	  status = KXREAD(buf2,realSize*n1*n2,src,stype+indx*2+1)
          offset2 = (srcPE2-thisNode)*n1
	  if (offset2 .LT. 0) then
	    offset2 = offset2 + npts
	  endif
        endif

	if (numNode .LE. 1) then
	  k = 0
	  do i = 0, n1-1
	    do j = 0, n2-1
		buf1(k) = a(i,j)
		k = k + 1
	    end do
	  end do
	  do i = 0, n1-1
	    do j = 0, n2-1
  	      a(i,j) = buf1(MOD(i+sh,n1)*n2+j)
	    end do
	  end do
	else
          do i = 0, n1-1
	    if ((i+sh) .LE. n1-1) then
	      do j = 0, n2-1
	        a(i,j) = a(i+sh,j)
	      end do
            else
              if (MOD(INT((i+minLoc+sh)/n1),numNode) .EQ. srcPE1) then
	        do j = 0, n2-1
  	          a(i,j) = buf1((i+sh-offset1)*n2+j)
	        end do
              else
	        do j = 0, n2-1
  	          a(i,j) = buf2((i+sh-offset2)*n2+j)
	        end do
	      endif
	    endif
	  end do
	endif

	return
	end

