program Fortran_Figure5_2_a

implicit none

save

 integer, parameter :: long = selected_real_kind(15,50)
 integer, parameter :: numbertimes=200
 integer, parameter :: noreal=5
 integer, parameter :: nrn=100000 
 real(kind=long), dimension(nrn) :: RanNum  
 real(kind=long),  dimension(numbertimes,noreal) :: Csave,timesave 
 real(kind=long), parameter :: k1=1.0_long 
 real(kind=long), parameter :: k2=50.0_long 
 real(kind=long), parameter :: k3=200.0_long 
 real(kind=long), parameter :: k4=0.1_long 
 
 real(kind=long) :: time,ss,a0,finaltime,tau,C,D
 integer :: ii,ir,jj,kk,maxnt
 integer size, seed(2)
 data seed /100,100/   
      size = 2
 
 call random_seed(SIZE=size)
 call random_seed(PUT=seed(1:size))
 call random_number(RanNum)
  
 ir=0
 maxnt=1
 finaltime=60.0_long
  
 do ii=1,numbertimes 
    do jj=1,noreal 
    Csave(ii,jj)=0.0_long
    timesave(ii,jj)=0.0_long
    end do
 end do
   
 do jj=1,noreal 
    ii=1
    C=0.0_long
    D=0.0_long
    ir=ir+1
    if (ir>nrn) then
       call random_number(RanNum)
       ir=1
    end if
    a0=k1+k4*k2*D/(k3+k2)
    tau=(1/a0)*log(1/RanNum(ir))
    time=tau			
    do while (time<finaltime)
       ir=ir+1
       if (ir>nrn) then
          call random_number(RanNum)
          ir=1
       end if
       ss=k1
       if (ss>RanNum(ir)*a0) then
	  D=D+1.0_long
       else
	  ss=ss+k4*k2*D/(k3+k2)
          if (ss>RanNum(ir)*a0) then
   	           D=D-1.0_long
		   C=C+1.0_long
	  else	   
	           write(*,*) 'Error.'
		   STOP
          end if
       end if
       ir=ir+1
       if (ir>nrn) then
          call random_number(RanNum)
          ir=1
       end if
       a0=k1+k4*k2*D/(k3+k2)
       tau=(1/a0)*log(1/RanNum(ir))
       time=time+tau			
       ii=ii+1
       Csave(ii,jj)=C
       timesave(ii,jj)=time       
    end do
    if (ii>maxnt) then
        maxnt=ii
    end if
    do kk=ii+1,numbertimes
       Csave(kk,jj)=C
       timesave(kk,jj)=time
    end do          
 end do
  	          
 open(unit=1, file='data_Figure5_2_a.dat',status='replace')

100 format(e15.6,x,e15.6)
 
    do jj=1,noreal 
       do ii=1,maxnt
          write(1,100) Csave(ii,jj),timesave(ii,jj)
       end do
    end do

 close(unit=1)

end program Fortran_Figure5_2_a