program Fortran_Figure5_2_c_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 
 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,relaxtime,tau,A,B,C,D,avfull
 integer :: ii,ir,jj,kk
 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
 relaxtime=0.01_long
 avfull=0.0_long
  
 do ii=1,numbermesh 
    histA(ii)=0.0_long
    histB(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-relaxtime))
       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			
    end do
    time=finaltime-relaxtime
    B=real(nint(k2*D/(k3+k2)))
    A=D-B
    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
    kk=0
    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
	        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
       avfull=avfull+1.0_long
       kk=kk+1
    end do
    histA(nint(A)+1)=histA(nint(A)+1)+1.0_long
    histB(nint(B)+1)=histB(nint(B)+1)+1.0_long
 end do
 
 write(*,*) 'Average number of reaction events of the full model:',avfull/real(noreal)
  	          
 open(unit=1, file='data_Figure5_2_c_d.dat',status='replace')

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

 close(unit=1)
 
end program Fortran_Figure5_2_c_d