! =====================================================================
! TREECODE: Fortran77 version of hierarchical N-body code (O(N log N)).

!     This program evolves a self-gravitating N-body system, using
!     a hierarchical O(N log N) algorithm to compute gravitational
!     forces (see Barnes, J.E. & Hut, P. 1986, Nature, 324, 446).
!     The computational system of units is determined by the input
!     data, with the assumption that G = 1.  Particles are not
!     required to have identical masses.

!     This code is based on a Fortran implementation of the
!     hierarchical algorithm developed and kindly provided by
!     Lars Hernquist.  Several additional refinments by Jun Makino
!     have been incorporated as well.

!     The present version of the code includes quadrupole moments
!     when calculating the gravitational fields of cells.  It also
!     incorporates an improved opening criterion which increases
!     the accuracy of the forces in certain exceptional situations.

! ==============================================================
! TREECD: top-level N-body evolution program.  Its tasks are to:
!  *  read input parameters and initial conditions;
!  *  perform an initial force calculation, reporting
!     system diagnostics and time required;
!  *  advance the system for NSTEPS time-steps, periodically
!     saving a snapshot of the system;
!  *  terminate the simulation and close data files.
! ==============================================================

         PROGRAM TREECD

	     INCLUDE 'treedefs.f90'
	     INTEGER NDIAG, N
	     data iprint /0/

		CALL MPI_INIT(ierr)
		CALL MPI_COMM_RANK(MPI_COMM_WORLD,my_rank,ierr)
		CALL MPI_COMM_SIZE(MPI_COMM_WORLD,nodes,ierr)
		
		Workers = nodes - 1
		IF (my_rank .NE. MASTER) CALL NODE_MAIN
		
!       ---------------------
!       Initialize cpu timer.
!       ---------------------
!          CALL SECOND(CPUT0)
!          CALL SECNDS(CPUT0)
! ************* CHANGE *******************
          CPUT0=MPI_WTIME()
!	----------------------------------------------------------------
!	Open data files, read input parameters and initial system state,
!       initialize system parameters, and start output.
!	----------------------------------------------------------------
          CALL INPARM
          CALL GETBDS
          CALL INITCD
          CALL BEGOUT
!	----------------------------------------------------------------
!	Distribute Parameters to nodes         
!	----------------------------------------------------------------
	CALL SEND_PARAMETERS
	
!   -----------------------------------------------
!	Generate tree and compute initial acceleration.
!	-----------------------------------------------
          CALL MKTREE
	      CALL ACCEL
!	------------------------------------------
!	Output timing data and system diagnostics.
!	------------------------------------------
          CALL OUTTIM
          CALL OUTLOG
          CALL DIAGS
!       -----------------------------------------------
!       Set number of steps between diagnostic outputs.
!       -----------------------------------------------
	      NDIAG = MAX(1, NSTEPS / 128)
!	=================================================================
!	Primary loop to advance system state for a given number of steps.
!	=================================================================
          DO N = 1, NSTEPS
!	  ------------------------------------------------------
!	  Advance velocity by 1/2 step, positions by 1 timestep.
!	  ------------------------------------------------------
          CALL ADVVEL(DTIME / 2.0)
          CALL ADVPOS(DTIME)
!	  ----------------------------
!	  Create tree, compute forces.
!	  ----------------------------
          CALL MKTREE
          CALL ACCEL
!	  ----------------------------------------------------
!	  Advance velocity 1/2 step to synchronize tvel, tpos.
!	  ----------------------------------------------------
          CALL ADVVEL(DTIME / 2.0)
!	  ----------------------------------------------------
!	  At appropriate intervals, output system diagnostics.
!	  ----------------------------------------------------
          IF (MOD(N, NDIAG) .EQ. 0) THEN
            CALL OUTLOG
            CALL DIAGS
	      ENDIF
!         ----------------------------------------------
!	  At appropriate intervals, output system state.
!         ----------------------------------------------
          IF (MOD(N, NOUT) .EQ. 0) &
           CALL PUTBDS(.FALSE.)
          END DO
!	----------------------------------
!	Output final snapshot with masses.
!	----------------------------------
          CALL PUTBDS(.TRUE.)
!	------------------------------------------------------------------
!	Stop timing, write timing data, close files, terminate simulation.
!	------------------------------------------------------------------
          CALL OUTCPU
          CALL STPOUT
          CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
          CALL MPI_FINALIZE(ierr)
          END

! ----------------------------------------------------------
! INITCD: initialize system parameters that depend on either
! the input data or defined PARAMETERS.
! ----------------------------------------------------------

        SUBROUTINE INITCD
 
	    INCLUDE 'treedefs.f90'
 
!	---------------------------------------
!	Initialize position and velocity times.
!	---------------------------------------
        TPOS = TNOW
        TVEL = TNOW
!	---------------------------------
!	Initialize box coordinate system.
!	---------------------------------
	    RSIZE = 4.0
        END

! ---------------------------------------------------------------
! ADVPOS: advance the positions of the bodies for a timestep DT.
! ---------------------------------------------------------------
 
        SUBROUTINE ADVPOS(DT)
        DOUBLE PRECISION::  DT

	    INCLUDE 'treedefs.f90'
        INTEGER K, P

!	------------------------------------------------- 
!	Loop over all spatial coordinates for all bodies.
!	------------------------------------------------- 
        DO K = 1, NDIM
          DO P = 1, NBODY
            POS(P,K) = POS(P,K) + DT*VEL(P,K)
          END DO
        END DO
!	----------------------------------
!	Update position time, system time.
!	----------------------------------
        TPOS = TPOS+DT
        TNOW = TPOS
        END

! ----------------------------------------------------------------
! ADVVEL: advance the velocities of the bodies for a timestep dt.
! ----------------------------------------------------------------
 
        SUBROUTINE ADVVEL(DT)
        REAL*8 DT

	    INCLUDE 'treedefs.f90'
        INTEGER K, P

!	------------------------------------------------- 
!	Loop over all spatial coordinates for all bodies.
!	------------------------------------------------- 
        DO K = 1, NDIM
          DO P = 1, NBODY
            VEL(P,K) = VEL(P,K) + DT*ACC(P,K)
          END DO
        END DO
!	----------------------------------
!	Update velocity time, system time.
!	----------------------------------
        TVEL = TVEL+DT
        TNOW = TVEL
        END

! ----------------------------------------------------------------
! DIAGS: compute diagnostics for the system: total energy, total
! kinetic energy, total potential energy, angular momentum, center
! of mass coordinates, and center of mass velocity.
! ----------------------------------------------------------------
 
        SUBROUTINE DIAGS

	    INCLUDE 'treedefs.f90'
        INTEGER P, K
 
!       --------------------------------------------------------------
!       Check that tpos and tvel are synchronized to one part in 1000.
!       --------------------------------------------------------------
        IF (ABS(TPOS - TVEL) .GT. DTIME/1.0E3) &
          CALL TERROR(' DIAGS: TIMES NOT SYNCHRONIZED')
!       ---------------------------------------------
!       Zero the accumulators for system diagnostics.
!       ---------------------------------------------
        MTOT = 0.0
        EKTOT = 0.0
        EPTOT = 0.0
        DO K = 1, NDIM
          CMPOS(K) = 0.0
          CMVEL(K) = 0.0
          AMVEC(K) = 0.0
        END DO
!       ----------------------------------------------------------
!       Loop over bodies to compute net mass and potential energy.
!       ----------------------------------------------------------
        DO P = 1, NBODY
          MTOT = MTOT + MASS(P)
          EPTOT = EPTOT + 0.5 * MASS(P) * PHI(P)
        END DO
!       -----------------------------------------------------------------
!       Compute net kinetic energy, center of mass position and velocity.
!       -----------------------------------------------------------------
        DO K = 1, NDIM
          DO P = 1, NBODY
            EKTOT = EKTOT + 0.5 * MASS(P) * VEL(P,K)**2
            CMPOS(K) = CMPOS(K) + MASS(P) * POS(P,K)
            CMVEL(K) = CMVEL(K) + MASS(P) * VEL(P,K)
          END DO
          CMVEL(K) = CMVEL(K) / MTOT
          CMPOS(K) = CMPOS(K) / MTOT
        END DO
!       ----------------------------
!       Compute total system energy.
!       ----------------------------
        ETOT = EKTOT + EPTOT
!       ---------------------------------------
!       Compute angular momentum of the system.
!       ---------------------------------------
        DO P = 1, NBODY
          AMVEC(1) = AMVEC(1) + MASS(P) * &
                      (POS(P,2)*VEL(P,3) - POS(P,3)*VEL(P,2))
          AMVEC(2) = AMVEC(2) + MASS(P) * &
                      (POS(P,3)*VEL(P,1) - POS(P,1)*VEL(P,3))
          AMVEC(3) = AMVEC(3) + MASS(P) * &
                      (POS(P,1)*VEL(P,2) - POS(P,2)*VEL(P,1))
        END DO
!       ----------------------------------
!       Write diagnostics to the log file.
!       ----------------------------------
        CALL OUTDAG
        END

! -------------------------------------------------------------
! TERROR: terminate the program as the result of a fatal error,
! closing the output files and writing timing information.
! -------------------------------------------------------------

        SUBROUTINE TERROR(MSG)
        CHARACTER*(*) MSG

	    INCLUDE 'treedefs.f90'

!       ------------------------------------
!       Write error message to the log file.
!       ------------------------------------
        CALL OUTERR(MSG)
!       ---------------------------------------------------
!       Output timing data, close files, terminate the run.
!       ---------------------------------------------------
        CALL OUTCPU
        CALL STPOUT
        END
        
!	----------------------------------------------------------------
!	NODE MAIN: Node Runtime area
!	----------------------------------------------------------------

	SUBROUTINE NODE_MAIN
	IMPLICIT NONE
	INCLUDE 'treedefs.f90'
	INTEGER :: P,K
	
	PRINT *, "Process ", my_rank, " has entered Node Main Subroutine"
	
	!	----------------------------------------------------------------
	!	Recieve Parameters from Master
	!	----------------------------------------------------------------
		CALL RECV_PARAMETERS
		
		
! N steps +   because there is the initial accel
! before master enters primary do loop

DO K = 1, NSTEPS + 1
	PRINT *, "The number of steps is : ", NSTEPS
	PRINT *, "On itertation ", K
	!	----------------------------------------------------------------
	!   Recieve Body Data from Master
	!	----------------------------------------------------------------
		CALL RECV_RAW_DATA
	
	!	----------------------------------------------------------------
	!	Reset Force Calculation Diagnostic Values
	!	----------------------------------------------------------------
			NCTOT = 0
			NBTOT = 0
			NTMAX = 0
			
	!	----------------------------------------------------------------
	!   Do each node's portion of the force calculations
	!	----------------------------------------------------------------

		DO P = mpi_index, mpi_index + workload(my_rank)
			CALL TRWALK(P)
		END DO
		
		
	!	----------------------------------------------------------------
	!   Return Processed Data to Master
	!	----------------------------------------------------------------
		CALL SEND_PROCESSED_DATA
END DO

CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
CALL MPI_FINALIZE(ierr)

STOP "NODE CALCULATION TERMINATED"

END SUBROUTINE NODE_MAIN