C  N-body calculation using the modified Euler method,
C  which is globally 2nd-order.
C  Uses a ring organization to calculate accelerations.
C
C  Data arrays (position, velocity, acceleration, mass)
C  are organized:
C
C          1	     Particle number (i) -->	       N
C	  +--------------------------------------------+
C   |  x  |					       |
C  (c) y  |					       |
C   v  z  |					       |
C	  +--------------------------------------------+
C
C  The parameter N and common blocks /NBODYFE/, /NBODYPE/ with arrays
C
C	X0, X1	- Current, next positions
C 	V0, V1	- velocities
C	 M, GM	- mass, G*mass (redundant)
C	 c,  i	- array coordinates
C  and
C	    NB	- integer <= N, actual no. of bodies
C	     G	- real, the gravitational constant
C
C  are defined in the file nbmode.h
C
C  File:  nbmode.fcm		Header:  nbmode.h


C------- setup() ---------------------------------------

      subroutine setup()

      INCLUDE 'nbmode.h'

C  Initialize coordinates
      c = SPREAD( [1:3],  2, N)
      i = SPREAD([1:N], 1, 3)
      end subroutine setup


C------- init() ----------------------------------------

C  Load initial positions, velocities, also masses,
C  gravitational const.  All are read from logical
C  unit 10.  Replicates masses over x,y,z.

      subroutine init()

      INCLUDE 'nbmode.h'

      integer k
      real mx
      X0 = 0.0		! For unused particles
      V0 = 0.0
      do, k = 1, NB
        read (10, *) X0(1, k), X0(2, k), X0(3, k)
      end do
      do, k = 1, NB
        read (10, *) V0(1, k), V0(2, k), V0(3, k)
      end do
      M = 0.0		! For unused particles
      do, k = 1, NB
        read (10, *) mx
        M(:, k) = mx
      end do
      read (10, *) G
      end subroutine init


C------- mkast(nx, ar, am) -----------------------------

C  Insert nx asteroids at radius ar with total mass am.
C  This assumes that an orbit of radius 1 about origin
C  corresponds to speed 2*pi (e.g., au-earth-yr units).

      subroutine mkast(nx, ar, am)
      integer nx	! Number of new asteroids
      real ar, am	! Orbit radius, total mass

      INCLUDE 'nbmode.h'

      real, parameter:: TWOPI = 6.2831853
      real, array(1:3, 1:N) :: sx, cx, arg
      real spd
      integer nn
      print 21, nx, ar, am
 21   FORMAT (" Making ", I4, " asteroids, radius ",
     *        F10.5, ", mass ", F8.5)
      nn = NB + nx
      spd = TWOPI/sqrt(ar)
      where ( (NB < i) .and. (i <= nn) )
        arg = (TWOPI*(i - NB))/nx
        sx = sin(arg)
        cx = cos(arg)
      end where
      where ( (NB < i) .and. (i <= nn) .and. (c == 1) )
        X0 =  cx
        V0 = -sx
      end where
      where ( (NB < i) .and. (i <= nn) .and. (c == 2) )
        X0 =  sx
        V0 =  cx
      end where
      where ( (NB < i) .and. (i <= nn) .and. (c == 3) )
        X0 = 0.0
        V0 = 0.0
      end where
      where ( (NB < i) .and. (i <= nn) )
        X0 = X0*ar
        V0 = V0*spd
         M = am/nx
      end where
      NB = nn
      end subroutine mkast


C------- print_state(PX) -------------------------------

C  Print current state of system according to PX.
C  See comments for main program for PX settings.

      subroutine print_state(PX)
      integer PX	! Print control

      INCLUDE 'nbmode.h'

      integer k, np
      np = NB
      if ( (PX >= 0) .and. (PX < NB) ) np = PX
      print 21
 21   FORMAT (/ "         x         y         z     ",
     *          "    vx        vy        vz        m")
      do, k = 1, np
        print 22, k, X0(1,k), X0(2,k), X0(3,k),
     *               V0(1,k), V0(2,k), V0(3,k), M(1,k)
      end do
 22   FORMAT (1X, I3, 6F10.5, E10.2)
      end subroutine print_state


C------- print_x(ax, na, nf) -------------------------------

C  Print ax(:, na:nf), nf-na <= 9 (for testing)

      subroutine print_x(ax, na, nf)
      INCLUDE 'nbmode.h'
      real, array(1:3, 1:N) :: ax
CMPF  ONDPU ax
CMF$  LAYOUT ax(,)
      integer na, nf

      integer u, lm
      lm = min(nf, na+9)
      do, u = 1, 3
        print 21, ax(u, na:lm)
 21     FORMAT ( 10(1X : E9.2 :) )
      end do
      end subroutine print_x

C------- print_c(Xf, Xc, Af, Ac, na, nb) -----------------------

C  Print conveyor status in positions na:nb (for testing)

      subroutine print_c(Xf, Xc, Af, Ac, na, nf)
      INCLUDE 'nbmode.h'
      real, array(1:3, 1:N) :: Xf, Xc, Af, Ac
      integer na, nf

      integer kk
      print 22, (kk, kk = na, nf)
 22   FORMAT ( 10(1X : I9 :) )
      print *, "  - Xf ----------"
      call print_x(Xf, na, nf)
      print *, "  - Xc ----------"
      call print_x(Xc, na, nf)
      print *, "  - Af ----------"
      call print_x(Af, na, nf)
      print *, "  - Ac ----------"
      call print_x(Ac, na, nf)
      end subroutine print_c


C------- Grav(Xf) --------------------------------------

C  Calculate accelerations (Newtonian gravitation)
C  Uses systolic ring calculation, taking advantage
C  of anti-symmetric force.  Data being circulated
C  on `conveyor' are marked `c', fixed items are `f'.

      function Grav(Xf)

      INCLUDE 'nbmode.h'

      real, array(1:3, 1:N) :: Xf, Grav
CMPF  ONDPU Xf
CMF$  LAYOUT Xf(,), Grav(,)
      real, array(1:3, 1:N) :: Af, Ac, D, R, Xc, Mc
      integer k, lm

      Af = 0.0		! Fixed accelerations
      Ac = 0.0		! Circulating accelerations
      Xc = Xf		! Circulating positions
      Mc = GM		! Circulating scaled masses
      lm = min(NB, N/2) - 1	! For loop control
      do, k = 1, lm		! 2-way accel. steps
	Xc = CSHIFT(Xc, SHIFT=-1, DIM=2)	! Step conveyor
	Mc = CSHIFT(Mc, SHIFT=-1, DIM=2)
	if ( k > 1 ) then
	  Ac = CSHIFT(Ac, SHIFT=-1, DIM=2)
	end if
C  Xc(:,i) = Xf(:,i-k), Similarly Mc and MG
C  Ac(:,i) contributes to particle i-k
	D = Xc - Xf
	R = sqrt(SPREAD(SUM(D*D, 1), 1, 3))
	WHERE ( GM > 0.0 .and. Mc > 0.0 )
	  D = D/R**3
	  Af = Af + Mc*D
	  Ac = Ac - GM*D
	end WHERE
      end do			! for 2-way accel.
C  Ac(i) corresponds to Af(i-lm)
      if ( NB > N/2 ) then	! Final 1-way accel.
	Xc = CSHIFT(Xc, SHIFT=-1, DIM=2)	! Step conveyor
	Mc = CSHIFT(Mc, SHIFT=-1, DIM=2)
	Ac = CSHIFT(Ac, SHIFT=-1, DIM=2)
	lm = N/2		! lm gives shift of Ac
	D = Xc - Xf
	R = sqrt(SPREAD(SUM(D*D, 1), 1, 3))
	WHERE ( GM > 0.0 .and. Mc > 0.0 )
	  D = D/R**3
	  Af = Af + Mc*D
	end WHERE
      end if			! for final accel.
D     call print_c(Xf, Xc, Af, Ac, 29, 35)
      Ac = CSHIFT(Ac, SHIFT=lm, DIM=2)    ! Back to correct position
D     call print_c(Xf, Xc, Af, Ac, 1, 7)
C   Combine accels for final result.
      Grav = Af + Ac
      end function Grav


C------- modeul(h, ns)  --------------------------------

C  Perform ns steps of modified Euler with time-step h.
C  X0, V0 are taken to be current state, and are updated
C  with final state.  X1, V1 are used in the process.

      subroutine modeul(h, ns)
      real h
      integer ns

      INCLUDE 'nbmode.h'

!       INTERFACE
!	  function Grav(Xf)
!	  !! INCLUDE 'nbmode.h'
!	  !real, array(1:3, 1:N) :: Xf, Grav
!CMPF    ONDPU Xf
!CMF$	LAYOUT Xf(,), Grav(,)
!       END ! INTERFACE

      real h2
      integer k

      h2 = 0.5*h	! h/2, for convenience
      do k = 1, ns	! ns steps
	X1 = X0 + h2*V0		! Half-way positions
	V1 = V0 + h2*Grav(X0)	! Half-way velocities
	X0 = X0 + h*V1		! Next positions
	V0 = V0 + h*Grav(X1)	! Next velocities
      end do
      end subroutine modeul


C------- nbmode  -------------------------------

C  Main program:  reads problem and control information
C  from data file named by user, performs and times
C  calculation.  Print control:
C	  0 - just performance data
C	< 0 - all
C	> 0 - number given

      program nbmode

      INCLUDE 'nbmode.h'

      integer ns, np, nv, k, PX
      real dt, ar, am
      character*50 fnm

      print '(" Data file> "$)'    ! Request data file
      read *, fnm
      print '(" Print ctl> "$)'	! Request print control
      read *, PX
C   Begin initialization
      open(10, file=fnm, status = "OLD")
      read (10, *) NB
      call setup()     ! Setup geometric structure
      call init()      ! Read data for individual bodies
C   Read "asteroid" descriptions from file
      read (10, *) ns
      do while ( ns > 0 )
        read (10, *) ar, am
        call mkast(ns, ar, am)
        read (10, *) ns
      end do
      read (10, *) ns, np, dt	! Integration parameters
      close(10)
      print 21, NB, N
 21   FORMAT (/, X, I4, " bodies, geometry allows ", I4)
      print 22, G, dt
 22   FORMAT ("    G = ", E13.6, "  delta-t = ", E12.5)
      print *, "  simulated time = 0"
      call print_state(PX)
C   Scale masses
      GM = G*m
C   np runs of ns steps, timestep dt
      do, k = 1, np
!        call CM_timer_clear(1)
!	call CM_timer_start(1)
!        call modeul(dt, ns)
!	call CM_timer_stop(1)
!	call CM_timer_print(1)
        print 23, k*ns*dt
 23     FORMAT (/ " simulated time = ", F9.3)
	call print_state(PX)
      end do
      end program nbmode

C*******************************************************
C  Data files:
C
C   N			No. of explicit bodies
C
C   x y z		Position of body 1
C     :                      :
C   x y z		Position of body N
C
C   vx vy va		Velocity of body 1
C      :                     :
C   vx vy va		Velocity of body N
C
C   m			Mass of body 1
C   :			      :
C   m			Mass of body N
C
C   G			Gravitational constant
C
C   nx			Number,
C   ar am		radius, mass of asteroid group
C      :                        :     0 or more groups
C   nx			Number,
C   ar am		radius, mass of asteroid group
C   -1			Stopper
C
C   ns			Number of time steps per phase
C   np			Number of phases
C   dt			delta-t (per time step)
C
C  Further data in the file will be ignored.
C*******************************************************
