program Fortran_Figure5_1_d

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) :: histA,histB,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,A,B,C
 real(kind=long) :: reaction1,reaction2,reaction3,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
 reaction2=0.0_long
 reaction3=0.0_long
 reaction4=0.0_long
  
 do ii=1,numbermesh 
    histA(ii)=0.0_long
    histB(ii)=0.0_long
    histC(ii)=0.0_long
 end do
 
 do jj=1,noreal 
    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
	  reaction1=reaction1+1.0_long
       else
          ss=ss+k2*A
          if (ss>RanNum(ir)*a0) then
             A=A-1.0_long
	     B=B+1.0_long
	     reaction2=reaction2+1.0_long
	  else
	     ss=ss+k3*B
             if (ss>RanNum(ir)*a0) then
	        B=B-1.0_long
	        A=A+1.0_long
	        reaction3=reaction3+1.0_long
	     else
	        ss=ss+k4*B
                if (ss>RanNum(ir)*a0) then
   	           B=B-1.0_long
		   C=C+1.0_long
	           reaction4=reaction4+1.0_long
	        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			
    end do
    histA(nint(A)+1)=histA(nint(A)+1)+1.0_long
    histB(nint(B)+1)=histB(nint(B)+1)+1.0_long
    histC(nint(C)+1)=histC(nint(C)+1)+1.0_long
 end do
 
 pom=reaction1+reaction2+reaction3+reaction4
 reaction1=reaction1/pom
 reaction2=reaction2/pom
 reaction3=reaction3/pom
 reaction4=reaction4/pom

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

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

 close(unit=1)
 

end program Fortran_Figure5_1_d