program Fortran_Figure6_3

implicit none

save

 integer, parameter :: long = selected_real_kind(15,50)
 integer, parameter :: numbermesh=40
 integer, parameter :: nrn=100000 
 real(kind=long), dimension(nrn) :: RanNum  
 real(kind=long),  dimension(numbermesh) :: A,B
 real(kind=long),  dimension(numbermesh) :: k1S,k2S,k5S,k6S,DhS 
 real(kind=long), parameter :: DiffConst=1.0_long ! diffusion constant
 real(kind=long), parameter :: k1c=0.03_long ! A+A->0
 real(kind=long), parameter :: k2c=0.3_long ! A+B->0
 real(kind=long), parameter :: k3c=0.0001_long ! A->0
 real(kind=long), parameter :: k4c=0.0001_long ! B->0  
 real(kind=long), parameter :: k5c=0.0000001_long ! 0->A
 real(kind=long), parameter :: k6c=0.000001_long ! 0->B
 real(kind=long) :: time,Dh,h,ss,a0,finaltime,ate,tau
 real(kind=long) :: k1,k2,k3,k4,k5,k6
 
 integer :: ii,ir,ns
 integer size, seed(2)
 data seed /2,2/ 
      size = 2
  
 h=25.0_long ! h=1000/real(numbermesh)
 Dh=DiffConst/(h*h)
 time=0_long
 ir=0
 a0=0.0_long
 finaltime=200.0_long*60.0_long
 
 k1=k1c/(h*h*h)
 k2=k2c/(h*h*h)
 k3=k3c
 k4=k4c
 k5=k5c*(h*h*h)
 k6=k6c*(h*h*h)
 
 call random_seed(SIZE=size)
 call random_seed(PUT=seed(1:size))
 call random_number(RanNum)
 
 do ii=1,numbermesh
    A(ii)=0.0_long
 end do
    
 do ii=2,numbermesh-1
    DhS(ii)=2.0_long*Dh
 end do
 ii=1
    DhS(ii)=Dh
 ii=numbermesh
    DhS(ii)=Dh
 
 do ii=1,numbermesh
    k1S(ii)=k1*A(ii)*(A(ii)-1)
    k2S(ii)=k2*A(ii)*B(ii)
    k5S(ii)=k5
    k6S(ii)=0.0_long
 end do
    
 do ii=3*numbermesh/5+1,numbermesh   
    k6S(ii)=k6
 end do
 
 do ii=1,numbermesh    
    a0=a0+DhS(ii)*A(ii)+DhS(ii)*B(ii)+k1S(ii)+k2S(ii)+k3*A(ii)+k4*B(ii)+k5S(ii)+k6S(ii)
 end do
 
 do while (time<finaltime)
    ir=ir+1
    if (ir>nrn) then
        call random_number(RanNum)
        ir=1
    end if
    tau=(1/a0)*log(1/RanNum(ir))
    ir=ir+1
    if (ir>nrn) then
        call random_number(RanNum)
        ir=1
    end if
    ss=0_long
    ii=0
    do while ((ss<RanNum(ir)*a0).AND.(ii<(numbermesh)))
       ii=ii+1
       ss=ss+k1S(ii)
    end do
    if (ss>RanNum(ir)*a0) then
        a0=a0-k1S(ii)-k2S(ii)
	A(ii)=A(ii)-2.0_long
        k1S(ii)=k1*A(ii)*(A(ii)-1)
        k2S(ii)=k2*A(ii)*B(ii)
        a0=a0+k1S(ii)+k2S(ii)-2.0_long*k3-2.0_long*DhS(ii)
    else
        ii=0
        do while ((ss<RanNum(ir)*a0).AND.(ii<(numbermesh)))
           ii=ii+1
	   ss=ss+k2S(ii)
	end do
        if (ss>RanNum(ir)*a0) then
           a0=a0-k1S(ii)-k2S(ii)
           A(ii)=A(ii)-1_long
           B(ii)=B(ii)-1_long
           k1S(ii)=k1*A(ii)*(A(ii)-1)
           k2S(ii)=k2*A(ii)*B(ii)
           a0=a0+k1S(ii)+k2S(ii)-k3-k4-2.0_long*DhS(ii)
	else
           ii=0
           do while ((ss<RanNum(ir)*a0).AND.(ii<(numbermesh)))
              ii=ii+1
	      ss=ss+k5S(ii)
	   end do
           if (ss>RanNum(ir)*a0) then
               a0=a0-k1S(ii)-k2S(ii)
               A(ii)=A(ii)+1_long
               k1S(ii)=k1*A(ii)*(A(ii)-1)
               k2S(ii)=k2*A(ii)*B(ii)
               a0=a0+k1S(ii)+k2S(ii)+k3+DhS(ii)
	   else
               ii=0
               do while ((ss<RanNum(ir)*a0).AND.(ii<(numbermesh)))
                  ii=ii+1
	          ss=ss+k6S(ii)
	       end do
               if (ss>RanNum(ir)*a0) then
                   a0=a0-k2S(ii)
                   B(ii)=B(ii)+1_long
                   k2S(ii)=k2*A(ii)*B(ii)
                   a0=a0+k2S(ii)+k4+DhS(ii)
 	       else
                   ii=0
                   do while ((ss<RanNum(ir)*a0).AND.(ii<(numbermesh-1)))
                      ii=ii+1
	              ss=ss+Dh*A(ii)
                   end do
                   if (ss>RanNum(ir)*a0) then
                       a0=a0-k1S(ii)-k2S(ii)-k1S(ii+1)-k2S(ii+1)
	               A(ii)=A(ii)-1_long
	               A(ii+1)=A(ii+1)+1_long
                       k1S(ii)=k1*A(ii)*(A(ii)-1)
                       k2S(ii)=k2*A(ii)*B(ii)
                       k1S(ii+1)=k1*A(ii+1)*(A(ii+1)-1)
                       k2S(ii+1)=k2*A(ii+1)*B(ii+1)
                       a0=a0+k1S(ii)+k2S(ii)+k1S(ii+1)+k2S(ii+1)
	               a0=a0+DhS(ii+1)-DhS(ii)
	           else
                       ii=1
                       do while ((ss<RanNum(ir)*a0).AND.(ii<(numbermesh)))
                          ii=ii+1
	                  ss=ss+Dh*A(ii)
                       end do
                       if (ss>RanNum(ir)*a0) then
                           a0=a0-k1S(ii)-k2S(ii)-k1S(ii-1)-k2S(ii-1)
		           A(ii)=A(ii)-1_long
	                   A(ii-1)=A(ii-1)+1_long   
                           k1S(ii)=k1*A(ii)*(A(ii)-1)
                           k2S(ii)=k2*A(ii)*B(ii)
                           k1S(ii-1)=k1*A(ii-1)*(A(ii-1)-1)
                           k2S(ii-1)=k2*A(ii-1)*B(ii-1)
                           a0=a0+k1S(ii)+k2S(ii)+k1S(ii-1)+k2S(ii-1)
		           a0=a0+DhS(ii-1)-DhS(ii)
		       else	   
                           ii=0
                           do while ((ss<RanNum(ir)*a0).AND.(ii<(numbermesh-1)))
                              ii=ii+1
	                      ss=ss+Dh*B(ii)
                           end do
                           if (ss>RanNum(ir)*a0) then
                               a0=a0-k2S(ii)-k2S(ii+1)
	                       B(ii)=B(ii)-1_long
	                       B(ii+1)=B(ii+1)+1_long
                               k2S(ii)=k2*A(ii)*B(ii)
                               k2S(ii+1)=k2*A(ii+1)*B(ii+1)
                               a0=a0+k2S(ii)+k2S(ii+1)+DhS(ii+1)-DhS(ii)
	                   else
                               ii=1
                               do while ((ss<RanNum(ir)*a0).AND.(ii<(numbermesh)))
                                  ii=ii+1
	                          ss=ss+Dh*B(ii)
                               end do
                               if (ss>RanNum(ir)*a0) then
                                   a0=a0-k2S(ii)-k2S(ii-1)
	                           B(ii)=B(ii)-1_long
	                           B(ii-1)=B(ii-1)+1_long
                                   k2S(ii)=k2*A(ii)*B(ii)
                                   k2S(ii-1)=k2*A(ii-1)*B(ii-1)
                                   a0=a0+k2S(ii)+k2S(ii-1)+DhS(ii-1)-DhS(ii)
	                       else
                                   ii=0
                                   do while ((ss<RanNum(ir)*a0).AND.(ii<(numbermesh)))
                                      ii=ii+1
	                              ss=ss+k3*A(ii)
	                           end do
                                   if (ss>RanNum(ir)*a0) then
                                       a0=a0-k1S(ii)-k2S(ii)
                                       A(ii)=A(ii)-1_long
                                       k1S(ii)=k1*A(ii)*(A(ii)-1)
                                       k2S(ii)=k2*A(ii)*B(ii)
                                       a0=a0+k1S(ii)+k2S(ii)-k3-DhS(ii)
	                           else
                                       ii=0
                                       do while ((ss<RanNum(ir)*a0).AND.(ii<(numbermesh)))
                                          ii=ii+1
	                                  ss=ss+k4*B(ii)
	                               end do
                                       if (ss>RanNum(ir)*a0) then
                                          a0=a0-k2S(ii)
                                          B(ii)=B(ii)-1_long
                                          k2S(ii)=k2*A(ii)*B(ii)
                                          a0=a0+k2S(ii)-k4-DhS(ii)
 	                               else
			                  write(*,*) 'Error.'
				       end if
				   end if   	  
	                       end if
			   end if
		       end if
		   end if    	       
	       end if	   	   
	   end if
        end if
    end if
    time=time+tau			
    if (int(0.01_long*time)-int(0.01_long*(time-tau))>0.5_long) then
        write(*,*) 'Time computed:',nint(time),'seconds (final time 18000 sec)'
    end if
 end do
 
	          
 open(unit=1, file='data_Figure6_3.dat', status='replace')

100 format(i6,X,i6)
 
    do ii=1,numbermesh 
        write(1,100) nint(A(ii)),nint(B(ii))
    end do

 close(unit=1)
 
end program Fortran_Figure6_3