program Fortran_Figure5_1_a_b_c

implicit none

save

 integer, parameter :: long = selected_real_kind(15,50)
 integer, parameter :: numbertimes=510
 integer, parameter :: noreal=5
 integer, parameter :: nrn=100000 
 real(kind=long), dimension(nrn) :: RanNum  
 real(kind=long),  dimension(numbertimes,noreal) :: Asave,Bsave,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,A,B,C
 integer :: ii,ir,jj,kk
 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
 finaltime=60.0_long
  
 do ii=1,numbertimes 
    do jj=1,noreal 
    Asave(ii,jj)=0.0_long
    Bsave(ii,jj)=0.0_long
    Csave(ii,jj)=0.0_long
    timesave(ii,jj)=0.0_long
    end do
 end do
  
 do jj=1,noreal 
    ii=1
    A=0.0_long
    B=0.0_long
    C=0.0_long
    ir=ir+1
    if (ir>nrn) then
       call random_number(RanNum)
       ir=1
    end if
    a0=k1+k2*A+(k3+k4)*B
    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
	  A=A+1.0_long
       else
          ss=ss+k2*A
          if (ss>RanNum(ir)*a0) then
             A=A-1.0_long
	     B=B+1.0_long
	  else
	     ss=ss+k3*B
             if (ss>RanNum(ir)*a0) then
	        B=B-1.0_long
	        A=A+1.0_long
	     else
	        ss=ss+k4*B
                if (ss>RanNum(ir)*a0) then
   	           B=B-1.0_long
		   C=C+1.0_long
		   ii=ii+1
                   Asave(ii,jj)=A
                   Bsave(ii,jj)=B
                   Csave(ii,jj)=C
                   timesave(ii,jj)=time
	        else	   
	           write(*,*) 'Error.'
		   STOP
	        end if	   	   
	     end if
          end if
       end if
       ir=ir+1
       if (ir>nrn) then
          call random_number(RanNum)
          ir=1
       end if
       a0=k1+k2*A+(k3+k4)*B
       tau=(1/a0)*log(1/RanNum(ir))
       time=time+tau			
       if (int(8.0_long*time)-int(8.0_long*(time-tau))>0.5_long) then
	  ii=ii+1
          Asave(ii,jj)=A
          Bsave(ii,jj)=B             
	  Csave(ii,jj)=C
          timesave(ii,jj)=time       
       end if
    end do
    do kk=ii+1,numbertimes
       Asave(kk,jj)=A
       Bsave(kk,jj)=B
       Csave(kk,jj)=C
       timesave(kk,jj)=time
    end do          
 end do
  	          
 open(unit=1, file='data_Figure5_1_a_b_c.dat',status='replace')

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

 close(unit=1)
 

end program Fortran_Figure5_1_a_b_c