program Fortran_Figure5_2_b

implicit none

save

 integer, parameter :: long = selected_real_kind(15,50)
 integer, parameter :: numbermesh=60
 integer, parameter :: noreal=1000000
 integer, parameter :: nrn=100000 
 real(kind=long), dimension(nrn) :: RanNum  
 real(kind=long),  dimension(numbermesh) :: histC 
 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
 real(kind=long) :: reaction1,reaction4,pom
 integer :: ii,ir,jj
 integer size, seed(2)
 data seed /10,10/   
      size = 2
 
 call random_seed(SIZE=size)
 call random_seed(PUT=seed(1:size))
 call random_number(RanNum)
  
 ir=0
 finaltime=50.0_long
 reaction1=0.0_long
 reaction4=0.0_long
  
 do ii=1,numbermesh 
    histC(ii)=0.0_long
 end do
 
 do jj=1,noreal 
    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
	  reaction1=reaction1+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
	           reaction4=reaction4+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			
    end do
    histC(nint(C)+1)=histC(nint(C)+1)+1.0_long
 end do
 
 pom=reaction1+reaction4
 reaction1=reaction1/pom
 reaction4=reaction4/pom

 write(*,*) 'Reaction ratios:',reaction1,reaction4
 	          
 open(unit=1, file='data_Figure5_2_b.dat',status='replace')

100 format(e15.6)
 
    do ii=1,numbermesh 
        write(1,100) histC(ii)/real(noreal)
    end do

 close(unit=1)
 

end program Fortran_Figure5_2_b