C     ==========================================================
C     NBODY PROBLEM
C     ==========================================================
      PROGRAM NBODY
      IMPLICIT NONE

      INTERFACE

        SUBROUTINE TIMER(T)
          DOUBLE PRECISION, INTENT(OUT) :: T
        END SUBROUTINE TIMER

        EXTRINSIC(HPF_LOCAL) SUBROUTINE L_TIMER(R_TIME,I_TIME) 
              REAL, DIMENSION(:), INTENT(IN)  :: I_TIME
              REAL, DIMENSION(:), INTENT(OUT) :: R_TIME
!HPF$         DISTRIBUTE *(BLOCK)             :: I_TIME
!HPF$         ALIGN R_TIME WITH I_TIME
        END SUBROUTINE L_TIMER

      END INTERFACE   
   
C     ==========================================================
C     VARIABLES , ALIGNED WRT THE PROBLEM'S SPACE
C     ==========================================================

      integer, parameter :: MM=0, PX=1, PY=2, PZ=3, FX=4, FY=5, FZ=6
!     MM:mass, PX:position in x coordinate, PY:position in y coordinate
!     FX: force in x coordinate, FY: force in y coordinate

      real, parameter    :: G = 1.0

      integer, parameter          :: NN = 1024 
      real, dimension(0:NN-1)     :: x,t,dx,dy,dz,tx,ty,tz,sq,dist,fac
      real, dimension(0:NN-1,0:6) :: p,q

!   Specify the HPf directives

!HPF$ DISTRIBUTE x(block)
!HPF$ ALIGN (:)   WITH x(:) :: t,dx,dy,dz,tx,ty,tz,sq,dist,fac
!HPF$ ALIGN (:,*) WITH x(:) :: p,q

      integer           :: i,j,k,npts
      logical           :: odd, IsOdd
 
		
C     ==========================================================
C     DEBUG AND PRINTING VARIABLES
C     ==========================================================

      DOUBLE PRECISION                       :: START_TIME
      DOUBLE PRECISION                       :: END_TIME
      REAL, ALLOCATABLE,DIMENSION(:)         :: ti, tf
!HPF$ DISTRIBUTE ti(BLOCK) 
!HPF$ DISTRIBUTE tf(BLOCK) 

      INTEGER, PARAMETER                     :: LEVEL0 = 0
      INTEGER, PARAMETER                     :: LEVEL1 = 1
      INTEGER, PARAMETER                     :: LEVEL2 = 2 
      INTEGER, PARAMETER                     :: DEBUG  = LEVEL0
 
C     ==========================================================
C     INITIALIZATION
C     ==========================================================

        print *,'********************************************'
        print *,'             NBODY PROBLEM'
        print *,'NUMBER OF POINTS ' , NN
        print *,'NUMBER OF PROCESSORS ', NUMBER_OF_PROCESSORS()
 
        ALLOCATE(ti(NUMBER_OF_PROCESSORS()))
        ALLOCATE(tf(NUMBER_OF_PROCESSORS()))

!       Start Timer

        CALL TIMER(START_TIME)
        call l_timer(tf, ti)

        call RANDOM_SEED
        call RANDOM_NUMBER(x)
        p(:,MM) = x(:)
        call RANDOM_NUMBER(x)
        p(:,PX) = x(:)
        call RANDOM_NUMBER(x)
        p(:,PY) = x(:)
        call RANDOM_NUMBER(x)
        p(:,PZ) = x(:)
        p(:,FX) = 0.0
        p(:,FY) = 0.0
        p(:,FZ) = 0.0

! Stop timer
        call l_timer(ti, tf)
        CALL TIMER(END_TIME)

        print *,'time taken for initialization is ',ti,'seconds'
        print *,'INITIALIZATION ELAPSED TIME ', END_TIME-START_TIME

C     ==========================================================
C     PROCESSING
C     ==========================================================

!       Start Timer
        ti = 0.0
        tf = 0.0
        CALL TIMER(START_TIME)
        call l_timer(tf, ti)

        npts = NN

        q = p
        k = npts/2 - 1
        do i = 0, k - 1
          q = CSHIFT(q,DIM=1,SHIFT=1)
          dx(0:npts-1) = p(0:npts-1,PX) - q(0:npts-1,PX)
          dy(0:npts-1) = p(0:npts-1,PY) - q(0:npts-1,PY)
          dz(0:npts-1) = p(0:npts-1,PZ) - q(0:npts-1,PZ)
          sq(0:npts-1) = dx(0:npts-1)**2+dy(0:npts-1)**2+
     &        dz(0:npts-1)**2
          dist(0:npts-1) = SQRT(sq(0:npts-1))
          fac(0:npts-1) = p(0:npts-1,MM) * q(0:npts-1,MM) 
     &                   / (dist(0:npts-1) * sq(0:npts-1))
          tx(0:npts-1) = fac(0:npts-1) * dx(0:npts-1)
          ty(0:npts-1) = fac(0:npts-1) * dy(0:npts-1)
          tz(0:npts-1) = fac(0:npts-1) * dz(0:npts-1)
          p(0:npts-1,FX) = p(0:npts-1,FX)-tx(0:npts-1)
          q(0:npts-1,FX) = q(0:npts-1,FX)+tx(0:npts-1)
          p(0:npts-1,FY) = p(0:npts-1,FY)-ty(0:npts-1)
          q(0:npts-1,FY) = q(0:npts-1,FY)+ty(0:npts-1)
          p(0:npts-1,FZ) = p(0:npts-1,FZ)-tz(0:npts-1)
          q(0:npts-1,FZ) = q(0:npts-1,FZ)+tz(0:npts-1)
        end do

        q = CSHIFT(q,DIM=1,SHIFT=1)
        k = npts/2
        dx(0:k-1) = p(0:k-1,PX) - q(0:k-1,PX)
        dy(0:k-1) = p(0:k-1,PY) - q(0:k-1,PY)
        dz(0:k-1) = p(0:k-1,PZ) - q(0:k-1,PZ)
        sq(0:k-1) = dx(0:k-1)**2+dy(0:k-1)**2+dz(0:k-1)**2
        dist(0:k-1) = SQRT(sq(0:k-1))
        fac(0:k-1) = p(0:k-1,MM) * q(0:k-1,MM) 
     &              / (dist(0:k-1) * sq(0:k-1))
        tx(0:k-1) = fac(0:k-1) * dx(0:k-1)
        ty(0:k-1) = fac(0:k-1) * dy(0:k-1)
        tz(0:k-1) = fac(0:k-1) * dz(0:k-1)
        p(0:k-1,FX) = p(0:k-1,FX)-tx(0:k-1)
        q(0:k-1,FX) = q(0:k-1,FX)+tx(0:k-1)
        p(0:k-1,FY) = p(0:k-1,FY)-ty(0:k-1)
        q(0:k-1,FY) = q(0:k-1,FY)+ty(0:k-1)
        p(0:k-1,FZ) = p(0:k-1,FZ)-tz(0:k-1)
        q(0:k-1,FZ) = q(0:k-1,FZ)+tz(0:k-1)

        q = CSHIFT(q,DIM=1,SHIFT=npts/2)
        p(:,FX) = p(:,FX) + q(:,FX)
        p(:,FY) = p(:,FY) + q(:,FY)
        p(:,FZ) = p(:,FZ) + q(:,FZ)
        
! Stop timer
        call l_timer(ti, tf)
        CALL TIMER(END_TIME)

        print *,'time taken is ',ti,'seconds'
        print *,'ELAPSED TIME PROCESSING ', END_TIME-START_TIME

C     ==========================================================
C     PRINT RESULTS
C     ==========================================================
 
        print *, 'FX', p(0:6,FX)
        print *, 'FY', p(0:6,FY)
        print *, 'FZ', p(0:6,FZ)

      STOP
      END PROGRAM NBODY


      SUBROUTINE TIMER(T)
      DOUBLE PRECISION, INTENT(OUT) :: T
      
      INTEGER, DIMENSION(8) :: VAL

        CALL DATE_AND_TIME(VALUES=VAL) 
        T = (DBLE(VAL(8)) * 1D-3) + (DBLE(VAL(7))) +
     &      (DBLE(VAL(6)) * 60D0) + (DBLE(VAL(5)) * 3600D0)

      RETURN
      END SUBROUTINE TIMER

C     ==========================================================
C     L_TIMER SUBROUTINE
C        RETURN NUMBER OF SECONDS IN EACH PROCESSOR
C     ==========================================================

      extrinsic(hpf_local) subroutine l_timer(r_time,i_time)
        real, dimension(:), intent(in)  :: i_time
        real, dimension(:), intent(out) :: r_time
        
        r_time(lbound(r_time, 1)) =
     &  secnds(i_time(lbound(i_time, 1)))

      end


