program FortranFigure4_1

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) :: A,k2S,DhS 
 integer, dimension(numbermesh,2) :: Asave 
 real(kind=long), parameter :: DiffConst=0.0001_long ! diffusion constant
 real(kind=long), parameter :: k1=0.001_long ! degradation rate
 real(kind=long), parameter :: k2=0.3_long ! creation rate
 
 real(kind=long) :: time,Dh,h,ss,a0,savingtime,finaltime,ate,tau
 integer :: ii,ir,ns
 logical :: savingdata
 integer size, seed(2)
 data seed /10,10/   
      size = 2
  
 h=0.025_long ! h=1/real(numbermesh)
 Dh=DiffConst/(h*h)
 time=0_long
 ir=0
 a0=0.0_long
 savingtime=10.0_long*60.0_long
 finaltime=30.0_long*60.0_long
 savingdata=.TRUE.
 
 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/5   
    k2S(ii)=k2
 end do
 do ii=numbermesh/5+1,numbermesh   
    k2S(ii)=0.0_long
 end do
 
 do ii=1,numbermesh    
    a0=a0+DhS(ii)*A(ii)+k1*A(ii)+k2S(ii)
 end do
 
 do while (time<finaltime)
    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+k1*A(ii)
    end do
    if (ss>RanNum(ir)*a0) then
	A(ii)=A(ii)-1_long
        a0=a0-k1-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
           A(ii)=A(ii)+1_long
           a0=a0+k1+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
	       A(ii)=A(ii)-1_long
	       A(ii+1)=A(ii+1)+1_long
	       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
	           A(ii)=A(ii)-1_long
	           A(ii-1)=A(ii-1)+1_long   
	           a0=a0+DhS(ii-1)-DhS(ii)
	       else	   
	           write(*,*) 'Error.'
		   STOP
	       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
        ate=0_long
        do ii=1,numbermesh    
           ate=ate+DhS(ii)*A(ii)+k1*A(ii)+k2S(ii)
        end do
        write(*,*) 'Time computed:',nint(time),'seconds (final time 1800 sec)'
    end if
    if ((time>savingtime).AND.savingdata) then
        ns=1
	savingdata=.FALSE.
        do ii=1,numbermesh
           Asave(ii,ns)=nint(A(ii))
	end do
    end if	
 end do
 
 ns=2
 do ii=1,numbermesh
    Asave(ii,ns)=nint(A(ii))
 end do
	          
 open(unit=1, file='dataFigure4_1.dat', status='replace')

100 format(i6,X,i6)
 
    do ii=1,numbermesh 
        write(1,100) Asave(ii,1),Asave(ii,2)
    end do

 close(unit=1)
 

end program FortranFigure4_1
