program Fortran_Figure6_4_b

implicit none

save

 integer, parameter :: long = selected_real_kind(15,50)
 integer, parameter :: ranmax=1000000
 integer, parameter :: maxocc=100
 integer, parameter :: maxK=36
 integer, parameter :: numberofrealizations=10000000

 real(kind=long), parameter :: D=10.0_long
 real(kind=long), parameter :: k1=1.0_long 
 real(kind=long), parameter :: k2=8.0_long
 real(kind=long), parameter :: L=1.0_long 
 real(kind=long), dimension(ranmax) :: RanNumUnif 
 integer, dimension(maxocc,4) :: positions ! compartment indices
 real(kind=long), dimension(maxK) :: averageA
 
 real(kind=long) :: time,h,Aone,k1h,diff,a0,savetime,ss,beta
 integer :: ii,iru,maxused,nr,total,K,jj,pom1,pom2,pom3,pp,ll
 logical :: pombool,pombool2
 
 integer size, seed(2)
 data seed /5, 5/ 
      size = 2

 call random_seed(SIZE=size)
 call random_seed(PUT=seed(1:size))
 call random_number(RanNumUnif)
 
 do ii=1,maxK
    averageA(ii)=0.0_long
 end do

 beta=0.275_long
 
 iru=0
 K=1
 
    time=0.0_long 
    savetime=0.01_long
 
    Aone=real(int(sqrt(k2/(2.0_long*k1))))
    
    a0=k1*Aone*(Aone-1.0_long)+k2
    iru=iru+1
    if (iru>ranmax) then
       call random_number(RanNumUnif)
       iru=1
    end if
    time=time+(1.0_long/a0)*log(1/RanNumUnif(iru))
                     
    do nr=1,numberofrealizations
       do while (time<real(nr*savetime))
          iru=iru+1
          if (iru>ranmax) then
             call random_number(RanNumUnif)
             iru=1
          end if
	  if (k2>RanNumUnif(iru)*a0) then
	      Aone=Aone+1.0_long
	  else 	      
	      Aone=Aone-2.0_long
	  end if     
          a0=k1*Aone*(Aone-1.0_long)+k2
          iru=iru+1
          if (iru>ranmax) then
             call random_number(RanNumUnif)
             iru=1
          end if
          time=time+(1.0_long/a0)*log(1/RanNumUnif(iru))
       end do
       averageA(1)=averageA(1)+Aone
    end do
    averageA(1)=averageA(1)/real(numberofrealizations)
      
    write(*,*) K,averageA(K)
  
 do K=2,maxK

    time=0.0_long 
  
    h=L/real(K)
    diff=D/(h*h)
    k1h=D*k1/(D*h*h*h-beta*h*h*k1)
      
    do ii=1,2 
       iru=iru+1
       if (iru>ranmax) then
          call random_number(RanNumUnif)
          iru=1
       end if
       pom1=int(RanNumUnif(iru)*K)+1
       iru=iru+1
       if (iru>ranmax) then
          call random_number(RanNumUnif)
          iru=1
       end if
       pom2=int(RanNumUnif(iru)*K)+1
       iru=iru+1
       if (iru>ranmax) then
          call random_number(RanNumUnif)
          iru=1
       end if
       pom3=int(RanNumUnif(iru)*K)+1
       if (ii.EQ.1) then
	   maxused=1
	   positions(maxused,1)=pom1
	   positions(maxused,2)=pom2
	   positions(maxused,3)=pom3
	   positions(maxused,4)=1
       else
           pombool=.TRUE.	   
	   jj=0  
	   do while (pombool.AND.(jj<maxused))
	      jj=jj+1
	      if (positions(jj,1).EQ.pom1) then
	         if (positions(jj,2).EQ.pom2) then
		    if (positions(jj,3).EQ.pom3) then
	               positions(jj,4)=positions(jj,4)+1
	               pombool=.FALSE.
		    endif
		 endif
              endif
	   enddo 
	   if (pombool) then
	      maxused=maxused+1
	      positions(maxused,1)=pom1
	      positions(maxused,2)=pom2
	      positions(maxused,3)=pom3
	      positions(maxused,4)=1
	   endif
       endif
    enddo   	   
	        		       
    do ii=3,maxocc
       positions(ii,1)=0
       positions(ii,2)=0
       positions(ii,3)=0
       positions(ii,4)=0
    end do 
    
    a0=k2
    do ii=1,maxused
       a0=a0+k1h*real(positions(ii,4)*(positions(ii,4)-1))+3.0_long*diff*real(positions(ii,4))
       if ((positions(ii,1)<K).AND.(positions(ii,1)>1)) then
          a0=a0+diff*positions(ii,4)
       endif
       if ((positions(ii,2)<K).AND.(positions(ii,2)>1)) then
          a0=a0+diff*positions(ii,4)
       endif
       if ((positions(ii,3)<K).AND.(positions(ii,3)>1)) then
          a0=a0+diff*positions(ii,4)
       endif
    enddo
          iru=iru+1
          if (iru>ranmax) then
             call random_number(RanNumUnif)
             iru=1
          end if
          time=time+(1.0_long/a0)*log(1/RanNumUnif(iru))
 
    do nr=1,numberofrealizations
       do while (time<real(nr*savetime))
          iru=iru+1
          if (iru>ranmax) then
             call random_number(RanNumUnif)
             iru=1
          end if
          ss=k2
	  if (ss>RanNumUnif(iru)*a0) then
             iru=iru+1
             if (iru>ranmax) then
                call random_number(RanNumUnif)
                iru=1
             end if
             pom1=int(RanNumUnif(iru)*K)+1
             iru=iru+1
             if (iru>ranmax) then
                call random_number(RanNumUnif)
                iru=1
             end if
             pom2=int(RanNumUnif(iru)*K)+1
             iru=iru+1
             if (iru>ranmax) then
                call random_number(RanNumUnif)
                iru=1
             end if
             pom3=int(RanNumUnif(iru)*K)+1
             if (maxused.EQ.0) then
	        maxused=1
   	        positions(maxused,1)=pom1
	        positions(maxused,2)=pom2
	        positions(maxused,3)=pom3
	        positions(maxused,4)=1
	     else
                pombool=.TRUE.	   
	        jj=0  
	        do while (pombool.AND.(jj<maxused))
	           jj=jj+1
	           if (positions(jj,1).EQ.pom1) then
	              if (positions(jj,2).EQ.pom2) then
		         if (positions(jj,3).EQ.pom3) then
  		            a0=a0-k1h*real(positions(jj,4)*(positions(jj,4)-1))
			    positions(jj,4)=positions(jj,4)+1
		            a0=a0+k1h*real(positions(jj,4)*(positions(jj,4)-1))
	                    pombool=.FALSE.
		         endif
		      endif
                   endif
	        enddo 
	        if (pombool) then
	           maxused=maxused+1
	           positions(maxused,1)=pom1
	           positions(maxused,2)=pom2
	           positions(maxused,3)=pom3
	           positions(maxused,4)=1
	        endif
             endif
	     a0=a0+3.0_long*diff
             if ((pom1<K).AND.(pom1>1)) then
                a0=a0+diff
             endif
             if ((pom2<K).AND.(pom2>1)) then
                a0=a0+diff
             endif
             if ((pom3<K).AND.(pom3>1)) then
                a0=a0+diff
             endif
          else   	   
             pombool=.TRUE.
	     ii=0   
	     do while (pombool.AND.(ii<maxused))
	        ii=ii+1
	        if (positions(ii,4)>1) then
	           ss=ss+k1h*real(positions(ii,4)*(positions(ii,4)-1))
	           if (ss>RanNumUnif(iru)*a0) then
		      a0=a0-k1h*real(positions(ii,4)*(positions(ii,4)-1))
		      positions(ii,4)=positions(ii,4)-2
		      a0=a0+k1h*real(positions(ii,4)*(positions(ii,4)-1))
		      a0=a0-6.0_long*diff
                      if ((positions(ii,1)<K).AND.(positions(ii,1)>1)) then
                         a0=a0-2.0_long*diff
                      endif
                      if ((positions(ii,2)<K).AND.(positions(ii,2)>1)) then
                         a0=a0-2.0_long*diff
                      endif
                      if ((positions(ii,3)<K).AND.(positions(ii,3)>1)) then
                         a0=a0-2.0_long*diff
                      endif
		      if (positions(ii,4).EQ.0) then
		         if (ii<maxused) then
			    do jj=ii+1,maxused
			       do ll=1,4
			          positions(jj-1,ll)=positions(jj,ll)
			       end do
			    end do   	      
			 endif
		         maxused=maxused-1
		      endif
		      pombool=.FALSE.
		   endif
                endif
	     enddo
	     ii=0    
	     do while (pombool.AND.(ii<maxused))
	        ii=ii+1
	        if (positions(ii,1)>1) then
	           ss=ss+diff*real(positions(ii,4))
	           if (ss>RanNumUnif(iru)*a0) then
		      pom1=positions(ii,1)-1
		      pom2=positions(ii,2)
		      pom3=positions(ii,3)
		      pombool2=.TRUE.	   
	              jj=0  
	              do while (pombool2.AND.(jj<maxused))
	                 jj=jj+1
	                 if (positions(jj,1).EQ.pom1) then
	                    if (positions(jj,2).EQ.pom2) then
		               if (positions(jj,3).EQ.pom3) then
  		                  a0=a0-k1h*real(positions(jj,4)*(positions(jj,4)-1))
			          positions(jj,4)=positions(jj,4)+1
		                  a0=a0+k1h*real(positions(jj,4)*(positions(jj,4)-1))
		                  a0=a0-k1h*real(positions(ii,4)*(positions(ii,4)-1))
			          positions(ii,4)=positions(ii,4)-1
		                  a0=a0+k1h*real(positions(ii,4)*(positions(ii,4)-1))
	                          pombool2=.FALSE.
		                  if (positions(ii,4).EQ.0) then
		                     if (ii<maxused) then
			                do pp=ii+1,maxused
			                   do ll=1,4
			                      positions(pp-1,ll)=positions(pp,ll)
			                   end do
			                end do   	      
			             endif
		                     maxused=maxused-1
                                  endif		         		      
   		               endif
		            endif
                         endif
	              enddo 
	              if (pombool2) then
		         if (positions(ii,4).EQ.1) then
			    positions(ii,1)=positions(ii,1)-1
			 else   
	                    maxused=maxused+1
	                    positions(maxused,1)=pom1
	                    positions(maxused,2)=pom2
	                    positions(maxused,3)=pom3
	                    positions(maxused,4)=1
		            a0=a0-k1h*real(positions(ii,4)*(positions(ii,4)-1))
			    positions(ii,4)=positions(ii,4)-1
		            a0=a0+k1h*real(positions(ii,4)*(positions(ii,4)-1))
	                 endif
		      endif
		      if (pom1.EQ.K-1) then
		         a0=a0+diff
		      end if
		      if (pom1.EQ.1) then
		         a0=a0-diff
		      endif	 		          
		   pombool=.FALSE.
		   endif
                endif
	     enddo
	     ii=0    
	     do while (pombool.AND.(ii<maxused))
	        ii=ii+1
	        if (positions(ii,1)<K) then
	           ss=ss+diff*real(positions(ii,4))
	           if (ss>RanNumUnif(iru)*a0) then
		      pom1=positions(ii,1)+1
		      pom2=positions(ii,2)
		      pom3=positions(ii,3)
		      pombool2=.TRUE.	   
	              jj=0  
	              do while (pombool2.AND.(jj<maxused))
	                 jj=jj+1
	                 if (positions(jj,1).EQ.pom1) then
	                    if (positions(jj,2).EQ.pom2) then
		               if (positions(jj,3).EQ.pom3) then
  		                  a0=a0-k1h*real(positions(jj,4)*(positions(jj,4)-1))
			          positions(jj,4)=positions(jj,4)+1
		                  a0=a0+k1h*real(positions(jj,4)*(positions(jj,4)-1))
		                  a0=a0-k1h*real(positions(ii,4)*(positions(ii,4)-1))
			          positions(ii,4)=positions(ii,4)-1
		                  a0=a0+k1h*real(positions(ii,4)*(positions(ii,4)-1))
	                          pombool2=.FALSE.
		                  if (positions(ii,4).EQ.0) then
		                     if (ii<maxused) then
			                do pp=ii+1,maxused
			                   do ll=1,4
			                      positions(pp-1,ll)=positions(pp,ll)
			                   end do
			                end do   	      
			             endif
		                     maxused=maxused-1
                                  endif		         		      
   		               endif
		            endif
                         endif
	              enddo 
	              if (pombool2) then
		         if (positions(ii,4).EQ.1) then
			    positions(ii,1)=positions(ii,1)+1
			 else   
	                    maxused=maxused+1
	                    positions(maxused,1)=pom1
	                    positions(maxused,2)=pom2
	                    positions(maxused,3)=pom3
	                    positions(maxused,4)=1
		            a0=a0-k1h*real(positions(ii,4)*(positions(ii,4)-1))
			    positions(ii,4)=positions(ii,4)-1
		            a0=a0+k1h*real(positions(ii,4)*(positions(ii,4)-1))
	                 endif
		      endif
		      if (pom1.EQ.2) then
		         a0=a0+diff
		      end if
		      if (pom1.EQ.K) then
		         a0=a0-diff
		      endif	 		          
		   pombool=.FALSE.
		   endif
                endif
	     enddo
	     ii=0    
	     do while (pombool.AND.(ii<maxused))
	        ii=ii+1
	        if (positions(ii,2)>1) then
	           ss=ss+diff*real(positions(ii,4))
	           if (ss>RanNumUnif(iru)*a0) then
		      pom1=positions(ii,1)
		      pom2=positions(ii,2)-1
		      pom3=positions(ii,3)
		      pombool2=.TRUE.	   
	              jj=0  
	              do while (pombool2.AND.(jj<maxused))
	                 jj=jj+1
	                 if (positions(jj,1).EQ.pom1) then
	                    if (positions(jj,2).EQ.pom2) then
		               if (positions(jj,3).EQ.pom3) then
  		                  a0=a0-k1h*real(positions(jj,4)*(positions(jj,4)-1))
			          positions(jj,4)=positions(jj,4)+1
		                  a0=a0+k1h*real(positions(jj,4)*(positions(jj,4)-1))
		                  a0=a0-k1h*real(positions(ii,4)*(positions(ii,4)-1))
			          positions(ii,4)=positions(ii,4)-1
		                  a0=a0+k1h*real(positions(ii,4)*(positions(ii,4)-1))
	                          pombool2=.FALSE.
		                  if (positions(ii,4).EQ.0) then
		                     if (ii<maxused) then
			                do pp=ii+1,maxused
			                   do ll=1,4
			                      positions(pp-1,ll)=positions(pp,ll)
			                   end do
			                end do   	      
			             endif
		                     maxused=maxused-1
                                  endif		         		      
   		               endif
		            endif
                         endif
	              enddo 
	              if (pombool2) then
		         if (positions(ii,4).EQ.1) then
			    positions(ii,2)=positions(ii,2)-1
			 else   
	                    maxused=maxused+1
	                    positions(maxused,1)=pom1
	                    positions(maxused,2)=pom2
	                    positions(maxused,3)=pom3
	                    positions(maxused,4)=1
		            a0=a0-k1h*real(positions(ii,4)*(positions(ii,4)-1))
			    positions(ii,4)=positions(ii,4)-1
		            a0=a0+k1h*real(positions(ii,4)*(positions(ii,4)-1))
	                 endif
		      endif
		      if (pom2.EQ.K-1) then
		         a0=a0+diff
		      end if
		      if (pom2.EQ.1) then
		         a0=a0-diff
		      endif	 		          
		   pombool=.FALSE.
		   endif
                endif
	     enddo
	     ii=0    
	     do while (pombool.AND.(ii<maxused))
	        ii=ii+1
	        if (positions(ii,2)<K) then
	           ss=ss+diff*real(positions(ii,4))
	           if (ss>RanNumUnif(iru)*a0) then
		      pom1=positions(ii,1)
		      pom2=positions(ii,2)+1
		      pom3=positions(ii,3)
		      pombool2=.TRUE.	   
	              jj=0  
	              do while (pombool2.AND.(jj<maxused))
	                 jj=jj+1
	                 if (positions(jj,1).EQ.pom1) then
	                    if (positions(jj,2).EQ.pom2) then
		               if (positions(jj,3).EQ.pom3) then
  		                  a0=a0-k1h*real(positions(jj,4)*(positions(jj,4)-1))
			          positions(jj,4)=positions(jj,4)+1
		                  a0=a0+k1h*real(positions(jj,4)*(positions(jj,4)-1))
		                  a0=a0-k1h*real(positions(ii,4)*(positions(ii,4)-1))
			          positions(ii,4)=positions(ii,4)-1
		                  a0=a0+k1h*real(positions(ii,4)*(positions(ii,4)-1))
	                          pombool2=.FALSE.
		                  if (positions(ii,4).EQ.0) then
		                     if (ii<maxused) then
			                do pp=ii+1,maxused
			                   do ll=1,4
			                      positions(pp-1,ll)=positions(pp,ll)
			                   end do
			                end do   	      
			             endif
		                     maxused=maxused-1
                                  endif		         		      
   		               endif
		            endif
                         endif
	              enddo 
	              if (pombool2) then
		         if (positions(ii,4).EQ.1) then
			    positions(ii,2)=positions(ii,2)+1
			 else   
	                    maxused=maxused+1
	                    positions(maxused,1)=pom1
	                    positions(maxused,2)=pom2
	                    positions(maxused,3)=pom3
	                    positions(maxused,4)=1
		            a0=a0-k1h*real(positions(ii,4)*(positions(ii,4)-1))
			    positions(ii,4)=positions(ii,4)-1
		            a0=a0+k1h*real(positions(ii,4)*(positions(ii,4)-1))
	                 endif
		      endif
		      if (pom2.EQ.2) then
		         a0=a0+diff
		      end if
		      if (pom2.EQ.K) then
		         a0=a0-diff
		      endif	 		          
		   pombool=.FALSE.
		   endif
                endif
	     enddo
	     ii=0    
	     do while (pombool.AND.(ii<maxused))
	        ii=ii+1
	        if (positions(ii,3)>1) then
	           ss=ss+diff*real(positions(ii,4))
	           if (ss>RanNumUnif(iru)*a0) then
		      pom1=positions(ii,1)
		      pom2=positions(ii,2)
		      pom3=positions(ii,3)-1
		      pombool2=.TRUE.	   
	              jj=0  
	              do while (pombool2.AND.(jj<maxused))
	                 jj=jj+1
	                 if (positions(jj,1).EQ.pom1) then
	                    if (positions(jj,2).EQ.pom2) then
		               if (positions(jj,3).EQ.pom3) then
  		                  a0=a0-k1h*real(positions(jj,4)*(positions(jj,4)-1))
			          positions(jj,4)=positions(jj,4)+1
		                  a0=a0+k1h*real(positions(jj,4)*(positions(jj,4)-1))
		                  a0=a0-k1h*real(positions(ii,4)*(positions(ii,4)-1))
			          positions(ii,4)=positions(ii,4)-1
		                  a0=a0+k1h*real(positions(ii,4)*(positions(ii,4)-1))
	                          pombool2=.FALSE.
		                  if (positions(ii,4).EQ.0) then
		                     if (ii<maxused) then
			                do pp=ii+1,maxused
			                   do ll=1,4
			                      positions(pp-1,ll)=positions(pp,ll)
			                   end do
			                end do   	      
			             endif
		                     maxused=maxused-1
                                  endif		         		      
   		               endif
		            endif
                         endif
	              enddo 
	              if (pombool2) then
		         if (positions(ii,4).EQ.1) then
			    positions(ii,3)=positions(ii,3)-1
			 else   
	                    maxused=maxused+1
	                    positions(maxused,1)=pom1
	                    positions(maxused,2)=pom2
	                    positions(maxused,3)=pom3
	                    positions(maxused,4)=1
		            a0=a0-k1h*real(positions(ii,4)*(positions(ii,4)-1))
			    positions(ii,4)=positions(ii,4)-1
		            a0=a0+k1h*real(positions(ii,4)*(positions(ii,4)-1))
	                 endif
		      endif
		      if (pom3.EQ.K-1) then
		         a0=a0+diff
		      end if
		      if (pom3.EQ.1) then
		         a0=a0-diff
		      endif	 		          
		   pombool=.FALSE.
		   endif
                endif
	     enddo
	     ii=0    
	     do while (pombool.AND.(ii<maxused))
	        ii=ii+1
	        if (positions(ii,3)<K) then
	           ss=ss+diff*real(positions(ii,4))
	           if (ss>RanNumUnif(iru)*a0) then
		      pom1=positions(ii,1)
		      pom2=positions(ii,2)
		      pom3=positions(ii,3)+1
		      pombool2=.TRUE.	   
	              jj=0  
	              do while (pombool2.AND.(jj<maxused))
	                 jj=jj+1
	                 if (positions(jj,1).EQ.pom1) then
	                    if (positions(jj,2).EQ.pom2) then
		               if (positions(jj,3).EQ.pom3) then
  		                  a0=a0-k1h*real(positions(jj,4)*(positions(jj,4)-1))
			          positions(jj,4)=positions(jj,4)+1
		                  a0=a0+k1h*real(positions(jj,4)*(positions(jj,4)-1))
		                  a0=a0-k1h*real(positions(ii,4)*(positions(ii,4)-1))
			          positions(ii,4)=positions(ii,4)-1
		                  a0=a0+k1h*real(positions(ii,4)*(positions(ii,4)-1))
	                          pombool2=.FALSE.
		                  if (positions(ii,4).EQ.0) then
		                     if (ii<maxused) then
			                do pp=ii+1,maxused
			                   do ll=1,4
			                      positions(pp-1,ll)=positions(pp,ll)
			                   end do
			                end do   	      
			             endif
		                     maxused=maxused-1
                                  endif		         		      
   		               endif
		            endif
                         endif
	              enddo 
	              if (pombool2) then
		         if (positions(ii,4).EQ.1) then
			    positions(ii,3)=positions(ii,3)+1
			 else   
	                    maxused=maxused+1
	                    positions(maxused,1)=pom1
	                    positions(maxused,2)=pom2
	                    positions(maxused,3)=pom3
	                    positions(maxused,4)=1
		            a0=a0-k1h*real(positions(ii,4)*(positions(ii,4)-1))
			    positions(ii,4)=positions(ii,4)-1
		            a0=a0+k1h*real(positions(ii,4)*(positions(ii,4)-1))
	                 endif
		      endif
		      if (pom3.EQ.2) then
		         a0=a0+diff
		      end if
		      if (pom3.EQ.K) then
		         a0=a0-diff
		      endif	 		          
		   pombool=.FALSE.
		   endif
                endif
	     enddo
	     if (pombool) then
	        write(*,*) 'Error. This should not happen.',ss,a0
		STOP
             endif
	  endif   		
          iru=iru+1
          if (iru>ranmax) then
             call random_number(RanNumUnif)
             iru=1
          end if
          time=time+(1.0_long/a0)*log(1/RanNumUnif(iru))
       enddo
       total=0
       do ii=1,maxused
          total=total+positions(ii,4)
       enddo
       averageA(K)=averageA(K)+real(total) 
    end do
    averageA(K)=averageA(K)/real(numberofrealizations)
    write(*,*) K,averageA(K)

    open(unit=1, file='data_Figure6_4_b.dat', status='replace')

    100 format(i6,X,e12.5)
 
    do ii=1,K
          write(1,100) ii,averageA(ii)
    end do

    close(unit=1)

end do

end program Fortran_Figure6_4_b