 module global 
    integer N, NB, switch_ini
    integer, parameter :: NX = 100000
    integer, parameter :: NBX = 1000
    double precision xmin, x_aver, Dx, Dxh, Dx_min, Dx_max
    double precision, dimension(NX) :: x
    double precision, dimension(NBX) :: P
 end module
!===================================================
 program main 
 use global
 integer N_ITER, NBN, ic
 double precision S
 character (len=4) iterstr

!====================================================
!  This code perform a microscopic energy-conservinhg 
!  exchange process in a system of N agents.
!
!  x = array, state of the system  
!  N = number of agents
!  NB = Number of bins for computing the entropy 
!  xmin = minimum allowed quantity of "substance"
!         per agent
!  x_aver = Fixed average value of "substance" in 
!           the system 
!  Dx = quantity of "substance" exchanged at each 
!       step (may be different at each step)
!  switch_ini : switch for initial condition 
!               switch_ini = 0 ... uniform IC
!                            x(i) = x_aver
!               switch_ini = 1 ... random IC
!                            x(i) = 2*x_aver*RAND()
!====================================================

 N = 10000
 xmin = 0.
 x_aver = 1.
 Dx_min = 0.1 
 Dx_max = 0.2
 Dx  = 0.1          
 DXh = 0.1
 NBN = 4*x_aver/DXh

 print*, ' # of bins necessary at most :', NBN
 if ( NBN > NBX ) then 
     print*, ' Max number of bins allowed  :', NBX
     print*, ' Please resize '
     STOP
 end if 

 N_ITER = 50000
 switch_ini = 0

 !==========================================
 ! Relaxation to equilibrium 

 call INI_COND()
 call entropy(S)
 open (2,file='entropy.dat')
 write (2,*) 0, S
 ic = 0
 do i=1, N_ITER
     call one_exch()
     if ( mod(i,1000) == 0) then 
         ic = ic + 1
         call entropy(S)
         write (2,*) i, S
         !**********************************
         write (iterstr,'(i4.4)') ic
         open (1,file='profile_'//iterstr)
             do ib=1, NB
                write (1,*) (ib-1+0.5)*DXh,&
                            P(ib)/DXh
             end do
         close (1)
         !**********************************
     end if
 end do
 close (2)
 !==========================================

 open (1,file='xfin.dat')
 do i=1, N
    write (1,*) x(i)
 end do
 close (1)

 end program 

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 subroutine INI_COND()
 use global 
 
 if ( switch_ini == 0 ) then 
    x = x_aver
 else 
    do i=1, N
       x(i) = 2*x_aver*RAND()
    end do
 end if
 
 end subroutine
 
!--------------------------------------------------
 subroutine one_exch()
 use global 
 integer loser, winner

 ! 1. select loser (i) and winner (j)
 
 loser = int( N*RAND() ) + 1
 
 ! perform x-conserving exchange if loser 
 ! has enough money
 ! 
 call gen_Dx_rand()  ! Comment for fixed Dx 
                     ! trades

 if ( x(loser) - xmin > Dx ) then
 10 winner = int( N*RAND() ) + 1 
   if ( winner == loser )  go to 10
 	x(loser)  = x(loser) - Dx
 	x(winner) = x(winner) + Dx
 end if	
 
 end subroutine
 !--------------------------------------------------
 subroutine gen_Dx_rand()
 use global
 
 Dx = RAND() * x_aver
 
 end subroutine
 
!=================================================== 
 
 subroutine entropy(S)
 use global 
 double precision S 

 ! Binning 
 NB = 0
 P = 0.d0
 do i=1, N
    ib = int((x(i) - xmin)/Dxh) + 1
    if ( ib > NB ) NB = ib
    P(ib) = P(ib) + 1.d0
 end do
 P = P/N

 ! Compute entropy 
 S = 0.d0
 do k=1, NB
    if ( P(k) > 0 ) S = S - P(k)*log(P(k))
 end do 

 end subroutine 




  


