program Fortran_Figure7_2_a

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,k2S,kjump 
 real(kind=long),  dimension(numbermesh) :: kplus,kminus 
 integer, dimension(numbermesh) :: 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.3125_long ! creation rate
 
 real(kind=long) :: time,h,ss,a0,finaltime,tau
 integer :: ii,ir
 logical :: savingdata
 integer size, seed(2)
 data seed /10,10/   
      size = 2
  
 h=0.025_long ! h=1/real(numbermesh)
 time=0_long
 ir=0
 a0=0.0_long
 finaltime=15.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=1,numbermesh/2
    kplus(ii)=DiffConst/(h*h)+0.001_long/(2.0_long*h)
    kminus(ii)=DiffConst/(h*h)-0.001_long/(2.0_long*h)
 end do
    
 do ii=numbermesh/2+1,numbermesh
    kplus(ii)=DiffConst/(h*h)-0.001_long/(2.0_long*h)
    kminus(ii)=DiffConst/(h*h)+0.001_long/(2.0_long*h)
 end do
    
 do ii=2,numbermesh-1
    kjump(ii)=kplus(ii)+kminus(ii)
 end do
 ii=1
    kjump(ii)=kplus(ii)
 ii=numbermesh
    kjump(ii)=kminus(ii)
    
 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+kjump(ii)*A(ii)+k1*A(ii)+k2S(ii)
 end do
    
 ir=ir+1
 if (ir>nrn) then
    call random_number(RanNum)
    ir=1
 end if
 tau=(1/a0)*log(1/RanNum(ir))
 time=time+tau
 
 do while (time<finaltime)
    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+k1*A(ii)
    end do
    if (ss>RanNum(ir)*a0) then
	A(ii)=A(ii)-1_long
        a0=a0-k1-kjump(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+kjump(ii)
	else
           ii=0
           do while ((ss<RanNum(ir)*a0).AND.(ii<(numbermesh-1)))
              ii=ii+1
	      ss=ss+kplus(ii)*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+kjump(ii+1)-kjump(ii)
	   else
               ii=1
               do while ((ss<RanNum(ir)*a0).AND.(ii<(numbermesh)))
                  ii=ii+1
	          ss=ss+kminus(ii)*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+kjump(ii-1)-kjump(ii)
	       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
    tau=(1/a0)*log(1/RanNum(ir))
    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 900 sec)'
    end if
 end do
 
 do ii=1,numbermesh
    Asave(ii)=nint(A(ii))
 end do
	          
 open(unit=1, file='data_Figure7_2_a.dat', status='replace')

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

 close(unit=1)
 

end program Fortran_Figure7_2_a