program Fortran_Figure3_5a

implicit none

save

 integer, parameter :: long = selected_real_kind(15,50)
 integer, parameter :: ranpom=1000000
 
 real(kind=long), parameter :: k1=0.18_long
 real(kind=long), parameter :: k2=0.00025_long
 real(kind=long), parameter :: k3=37.5_long
 real(kind=long), parameter :: k4=2200.0_long
 real(kind=long), dimension(2,ranpom) :: RanNumUnif 
 real(kind=long), dimension(600) :: histogram 
 real(kind=long) :: X,ss,a0,time,tau,nsave
 integer :: ii,iru,tt,oldtt
 
 integer size, seed(2)
 data seed /100, 100/   !(you can choose different values here)
      size = 2

 call random_seed(SIZE=size)
 call random_seed(PUT=seed(1:size))
 call random_number(RanNumUnif)
 
 iru=0
 time=0.0_long
 nsave=0.0_long
 X=100.0_long
   
 do ii=1,600
    histogram(ii)=0.0_long
 end do  
 
 oldtt=0
   
 do while (time<10000000.0_long)
          
       a0=k1*X*(X-1)+k2*X*(X-1)*(X-2)+k3*X+k4
       iru=iru+1
       if (iru>ranpom) then
           call random_number(RanNumUnif)
	   iru=1
       endif
       tau=(1/a0)*log(1/RanNumUnif(1,iru))
       time=time+tau
       ss=k1*X*(X-1)
       if (ss>RanNumUnif(2,iru)*a0) then
          X=X+1
       else
          ss=ss+k2*X*(X-1)*(X-2)
          if (ss>RanNumUnif(2,iru)*a0) then
             X=X-1
          else
             ss=ss+k3*X     
             if (ss>RanNumUnif(2,iru)*a0) then
                 X=X-1
             else
                 ss=ss+k4
                 if (ss>RanNumUnif(2,iru)*a0) then
                    X=X+1
                 else
                    STOP
                 endif
             endif    
          endif       
       endif
       tt=nint(10.0*time)
       if (tt>oldtt) then
          ii=nint(X)
          histogram(ii)=histogram(ii)+1.0_long           
          nsave=nsave+1.0_long
	  oldtt=tt
       end if
 end do   

 do ii=1,600
    histogram(ii)=histogram(ii)/nsave
 end do  
 
 open(unit=1, file='data_Figure3_5a.dat', status='replace')

100 format(e14.8)
 
 do ii=1,600
    write(1,100) histogram(ii)
 end do

 close(unit=1)
 
end program Fortran_Figure3_5a
