program FortranFigure5_4

implicit none

save

 integer, parameter :: long = selected_real_kind(15,50)
 integer, parameter :: numbermesh=40
 integer, parameter :: nahbon=100000 
 real(kind=long), dimension(nahbon) :: RanNum  
 real(kind=long), dimension(numbermesh) :: x,cubic,Ur,Vr,U,V 
 real(kind=long) :: Du,Dv,time,tau,h,Dpu,Dpv,ss,a0,timeend
 real(kind=long) :: pom,k1,k2,k3,k4
 integer :: ii,ir
 integer size, seed(2)
 data seed /20, 20/ 
      size = 2
  
 Du=0.00001_long                      ! Du = 10^{-5} mm^2 / sec
 Dv=0.001_long                        ! Dv = 10^{-3} mm^2 / sec  
 h=0.025_long ! h=1/real(numbermesh)  ! L = 1 mm, h=0.025 mm
 
 Dpu=2_long*Du/(h*h)
 Dpv=2_long*Dv/(h*h)
 
 k1=1_long                             
 k2=3_long
 k3=0.02_long
 k4=0.000001_long
 
 time=0_long
 timeend=30*60_long
 ir=0
 a0=0_long
 
 call random_seed(SIZE=size)
 call random_seed(PUT=seed(1:size))
 call random_number(RanNum)
 
 do ii=1,numbermesh
    x(ii)=h*(real(ii)-0.5_long)
    U(ii)=200_long
    V(ii)=75_long
    cubic(ii)=k4*U(ii)*(U(ii)-1_long)*V(ii)
    a0=a0+cubic(ii)+k3*U(ii)+k1+k2+U(ii)*Dpu+V(ii)*Dpv
 end do
 
 do while (time<timeend)
    ir=ir+1
    if (ir>nahbon) then
        call random_number(RanNum)
        ir=1
    end if
    tau=(1/a0)*log(1/RanNum(ir))
    ir=ir+1
    if (ir>nahbon) 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+cubic(ii)
    end do
    if (ss>RanNum(ir)*a0) then
        a0=a0-cubic(ii)
        U(ii)=U(ii)+1_long
	V(ii)=V(ii)-1_long
	cubic(ii)=k4*U(ii)*(U(ii)-1_long)*V(ii)
        a0=a0+cubic(ii)+k3+Dpu-Dpv
    else
        ii=0
        do while ((ss<RanNum(ir)*a0).AND.(ii<(numbermesh)))
           ii=ii+1
	   ss=ss+k3*U(ii)
	end do
        if (ss>RanNum(ir)*a0) then
           a0=a0-cubic(ii)
           U(ii)=U(ii)-1_long
	   cubic(ii)=k4*U(ii)*(U(ii)-1_long)*V(ii)
           a0=a0+cubic(ii)-k3-Dpu
	else
           ii=0
           do while ((ss<RanNum(ir)*a0).AND.(ii<(numbermesh)))
              ii=ii+1
	      ss=ss+k1
	   end do
           if (ss>RanNum(ir)*a0) then
               a0=a0-cubic(ii)
               U(ii)=U(ii)+1_long
	       cubic(ii)=k4*U(ii)*(U(ii)-1_long)*V(ii)
               a0=a0+cubic(ii)+k3+Dpu
	   else
               ii=0
               do while ((ss<RanNum(ir)*a0).AND.(ii<(numbermesh)))
                  ii=ii+1
	          ss=ss+k2
	       end do
               if (ss>RanNum(ir)*a0) then
                   a0=a0-cubic(ii)
                   V(ii)=V(ii)+1_long
	           cubic(ii)=k4*U(ii)*(U(ii)-1_long)*V(ii)
                   a0=a0+cubic(ii)+Dpv
	       else
                   ii=0
                   do while ((ss<RanNum(ir)*a0).AND.(ii<(numbermesh)))
                      ii=ii+1
	              ss=ss+U(ii)*Dpu
		   end do
                   if (ss>RanNum(ir)*a0) then
	               ir=ir+1
                       if (ir>nahbon) then
                           call random_number(RanNum)
                           ir=1
                       end if
		       if (RanNum(ir)<0.5_long) then
		          if (ii>1) then
                              a0=a0-cubic(ii)
                              a0=a0-cubic(ii-1)
		              U(ii)=U(ii)-1_long
			      U(ii-1)=U(ii-1)+1_long   
    	                      cubic(ii)=k4*U(ii)*(U(ii)-1_long)*V(ii)
    	                      cubic(ii-1)=k4*U(ii-1)*(U(ii-1)-1_long)*V(ii-1)
                              a0=a0+cubic(ii)
                              a0=a0+cubic(ii-1)
		          end if
		       else   
		          if (ii<numbermesh) then
                              a0=a0-cubic(ii)
                              a0=a0-cubic(ii+1)
		              U(ii)=U(ii)-1_long
			      U(ii+1)=U(ii+1)+1_long   
    	                      cubic(ii)=k4*U(ii)*(U(ii)-1_long)*V(ii)
    	                      cubic(ii+1)=k4*U(ii+1)*(U(ii+1)-1_long)*V(ii+1)
                              a0=a0+cubic(ii)
                              a0=a0+cubic(ii+1)
		          end if
		       end if   
		   else
                       ii=0
                       do while ((ss<RanNum(ir)*a0).AND.(ii<(numbermesh)))
                          ii=ii+1
	                  ss=ss+V(ii)*Dpv
		       end do
		       if (ss>RanNum(ir)*a0) then
                           ir=ir+1
                           if (ir>nahbon) then
                               call random_number(RanNum)
                               ir=1
                           end if
		           if (RanNum(ir)<0.5_long) then
		               if (ii>1) then
				   a0=a0-cubic(ii)
                                   a0=a0-cubic(ii-1)
		                   V(ii)=V(ii)-1_long
			           V(ii-1)=V(ii-1)+1_long   
    	                           cubic(ii)=k4*U(ii)*(U(ii)-1_long)*V(ii)
    	                           cubic(ii-1)=k4*U(ii-1)*(U(ii-1)-1_long)*V(ii-1)
                                   a0=a0+cubic(ii)
                                   a0=a0+cubic(ii-1)
		               end if
		           else   
		               if (ii<numbermesh) then
                                   a0=a0-cubic(ii)
                                   a0=a0-cubic(ii+1)
		                   V(ii)=V(ii)-1_long
			           V(ii+1)=V(ii+1)+1_long   
    	                           cubic(ii)=k4*U(ii)*(U(ii)-1_long)*V(ii)
    	                           cubic(ii+1)=k4*U(ii+1)*(U(ii+1)-1_long)*V(ii+1)
                                   a0=a0+cubic(ii)
                                   a0=a0+cubic(ii+1)
		               end if
                           end if        
		       else	   
			   write(*,*) 'Error.'
			   STOP
		       end if	   	   
	           end if
	       end if	   		
	   end if		
        end if
    end if			
    time=time+tau
    if (int(time*0.01_long)-int(time*0.01_long-tau*0.01_long)>0.5_long) then
        write(*,*) 'Time computed:',nint(time),'seconds (final time 1800 sec)'
    end if
 end do				   
	          
 open(unit=1, file='dataFigure5_4.dat', status='replace')

100 format(f10.6,X,f10.6,X,f10.6)
 
    do ii=1,numbermesh 
        Ur(ii)=U(ii)
        Vr(ii)=V(ii)
        write(1,100) x(ii),Ur(ii),Vr(ii)
    end do

 close(unit=1)
 

end program FortranFigure5_4