Module Nbody        ! Not the best code, but a very short one
integer,  parameter               :: Np =100
Real*8, DIMENSION(3,Np)  :: X,V,G
Real*8, DIMENSION(Np)     :: Mass
Real*8, Parameter                :: eps =1.d-3, eps2=eps**2, dt =1.d-4, t_end =1.
Real*8                                     :: Ekin, Epot,t=0.
Contains
Subroutine Accel              !-------------------------------------
   G = 0.
   Do i=1,Np       
   Do j=1,Np
        G(:,i) =G(:,i) +Mass(j)*(X(:,j)-X(:,i))/sqrt(  SUM( (X(:,j)-X(:,i))**2+eps2)**3  )
    EndDo
    EndDo
end Subroutine Accel
Subroutine Energy              !-------------------------------------
   Epot =0.
   Ekin = 0.5*SUM(Mass*SUM(V**2,1))
   Do i=1,Np-1     
   Do j=i+1,Np
        Epot = Epot-Mass(i)*Mass(j)/sqrt(  (SUM( (X(:,j)-X(:,i))**2+eps2)  ))
    EndDo
    EndDo
    write(*,'(4(2x,a,g12.5))') 't=',t,' Ekin=',Ekin,' Epot=',Epot,' Etot=',Ekin+Epot
end Subroutine Energy
end Module Nbody         
Program MakeItSimple     !-------------------------------------
Use Nbody
   Mass =1./Np
   Call random_number(X)
   Call random_number(V)
   X = 2.*X-1.;    V = (V-0.5)/3.
    Do
         If(mod(INT(t/dt),100)==0)Call Energy
         Call Accel
         V = V+G*dt ;    X = X+V*dt;   t = t +dt
         If(t> t_end)exit
     EndDo
End Program MakeItSimple