program Fortran_Figure9_1_b

implicit none

save

 integer, parameter :: long = selected_real_kind(15,50)
 integer, parameter :: molnum=10  ! total number of molecules
 integer, parameter :: ranpom=1000000
 integer, parameter :: ranpom2=500000
 integer, parameter :: numbercomp=10
 integer, parameter :: numberhist=20
 integer, parameter :: timesteps=10000000
 
 real(kind=long), parameter :: DiffConst=1.0_long ! diffusion constant
 real(kind=long), parameter :: dt=0.0001_long ! timestep
 
 real(kind=long), dimension(2,ranpom2) :: RanNum 
 real(kind=long), dimension(ranpom) :: RanNumNorm 
 real(kind=long), dimension(numbercomp) :: A,deltaA 
 real(kind=long), dimension(numberhist) :: statdist,xvalue
 real(kind=long), dimension(ranpom) :: RanNumUnif 
 real(kind=long), dimension(molnum) :: pos
 logical, dimension(molnum) :: inuse
 integer, dimension(molnum) :: indexempty
 
 integer :: jj,zz,numbermol,iru,irn,emptyuse,lastmol,ntime,qq,ii
 real(kind=long) :: pom1,pom2,pom3,time,alpha,beta,h,initialposition,Phi,xi1,xi2,xi3,oldpos
 
 integer size, seed(2)
 data seed /100, 100/ 
      size = 2

 call random_seed(SIZE=size)
 call random_seed(PUT=seed(1:size))
 call random_number(RanNum)
 call random_number(RanNumUnif)
 
 do zz=1,ranpom2
    pom1=sqrt(2.0_long*log(1.0_long/RanNum(1,zz)))
    pom2=2.0_long*3.14159265358979_long*RanNum(2,zz)
    RanNumNorm(2*zz-1)=pom1*cos(pom2)
    RanNumNorm(2*zz)=pom1*sin(pom2)
 end do   
 
 irn=0
 iru=0
 
 h=1.0_long/real(numbercomp)
 alpha=sqrt(2.0_long*DiffConst*dt)
 beta=dt*DiffConst/(h*h)
 Phi=1.0_long
   
 numbermol=5
 lastmol=numbermol      
  
 do jj=1,numbermol
    pos(jj)=0.1_long+0.2_long*real(jj-1)
    inuse(jj)=.TRUE.
 end do
  
 do jj=numbermol+1,molnum
    pos(jj)=0.5_long
    inuse(jj)=.FALSE.
 end do
 
 emptyuse=0 
 do jj=1,molnum
    indexempty(jj)=0
 end do
 
 do jj=1,numbercomp
    A(jj)=0.0_long
 end do
                  
 do jj=1,(molnum-numbermol)
    iru=iru+1
    if (iru>ranpom) then
       call random_number(RanNumUnif)
       iru=1
    end if    
    qq=1+int(RanNumUnif(iru)*numbercomp)
    A(qq)=A(qq)+1.0_long
 end do
    
 time=0.0_long
     
! main loop over timesteps time steps
  
 do ntime=1,timesteps
 
       time=time+dt

       ! initialize deltaA, step (a22)
       
       do jj=1,numbercomp
          deltaA(jj)=0.0_long
       end do

       ! evolve compartments, jumps to the left, step (b22)
       
       do jj=1,numbercomp-1
          iru=iru+1
          if (iru>ranpom) then
              call random_number(RanNumUnif)
              iru=1
          end if    
	  if (RanNumUnif(iru)<beta*A(jj)) then
	     deltaA(jj)=deltaA(jj)-1.0_long
	     deltaA(jj+1)=deltaA(jj+1)+1.0_long
	  end if
       end do
              
       ! evolve compartments, jumps to the left, step (c22)
       
       do jj=2,numbercomp
          iru=iru+1
          if (iru>ranpom) then
              call random_number(RanNumUnif)
              iru=1
          end if    
	  if (RanNumUnif(iru)<beta*A(jj)) then
	     deltaA(jj)=deltaA(jj)-1.0_long
	     deltaA(jj-1)=deltaA(jj-1)+1.0_long
	  end if
       end do
       
       do jj=1,molnum 
          
	  if (inuse(jj)) then             
       
	      ! evolve positions of molecules over one time step, step (d22)
	      
	      irn=irn+1
              if (irn>ranpom) then
                 call random_number(RanNum)
	         do zz=1,ranpom2
                    pom1=sqrt(2.0_long*log(1.0_long/RanNum(1,zz)))
                    pom2=2.0_long*3.14159265358979_long*RanNum(2,zz)
	            RanNumNorm(2*zz-1)=pom1*cos(pom2)
	            RanNumNorm(2*zz)=pom1*sin(pom2)
	         end do   
                 irn=1
              end if
              oldpos=pos(jj)
              pos(jj)=pos(jj)+alpha*RanNumNorm(irn)
	      
	      ! boundary condition at x=1 - no flux, step (e22)
	      	      	      
              if (pos(jj)>1.0_long) then
	          pos(jj)=2.0_long-pos(jj)
              end if
	      	      
	      ! boundary condition at x=0, steps (f22) and (g22)
	      
              if (pos(jj)<0.0_long) then
	          inuse(jj)=.FALSE.
		  emptyuse=emptyuse+1
		  numbermol=numbermol-1
		  indexempty(emptyuse)=jj
		  deltaA(1)=deltaA(1)+1.0_long
	      else
                  iru=iru+1
                  if (iru>ranpom) then
                     call random_number(RanNumUnif)
                     iru=1
                  end if    
                  if (RanNumUnif(iru)<exp(-pos(jj)*oldpos/(DiffConst*dt))) then
	             inuse(jj)=.FALSE.
		     emptyuse=emptyuse+1
		     numbermol=numbermol-1
		     indexempty(emptyuse)=jj
                     deltaA(1)=deltaA(1)+1.0_long			
		  end if       
              end if
	      
	  end if    
       end do
       
                    
       ! transfer from compartments to BD, step (h22)
       
       iru=iru+1
       if (iru>ranpom) then
          call random_number(RanNumUnif)
          iru=1
       end if    
       if (RanNumUnif(iru)<beta*Phi*A(1)) then
          deltaA(1)=deltaA(1)-1.0_long            
          iru=iru+1
          if (iru>ranpom) then
              call random_number(RanNumUnif)
              iru=1
          end if
          initialposition=h*RanNumUnif(iru)
	  if (emptyuse>0) then
	      pos(indexempty(emptyuse))=initialposition
              inuse(indexempty(emptyuse))=.TRUE.
              numbermol=numbermol+1
	      emptyuse=emptyuse-1
          else
	      lastmol=lastmol+1
	      pos(lastmol)=initialposition
              inuse(lastmol)=.TRUE.
              numbermol=numbermol+1
              if (lastmol>molnum) then
	         write(*,*) lastmol,molnum,time
		 pom1=0.0_long
		 do zz=1,numbercomp
		    pom1=pom1+A(zz)
		 end do   
		 write(*,*) pom1
	         write(*,*) 'Error! This should not happen.'
                 STOP
	      end if	 
	  end if
       end if	  	   
       
       ! update number of molecules in all compartments, step (i22)

       do jj=1,numbercomp
          A(jj)=A(jj)+deltaA(jj)
       end do
                  
       do ii=1,numbercomp
          statdist(numbercomp+1-ii)=statdist(numbercomp+1-ii)+A(ii)
       end do	  
          
       do jj=1,molnum
	  if (inuse(jj)) then
	     qq=numbercomp+1+int(pos(jj)*real(numberhist-numbercomp))            
             if (qq>numberhist) then
	         qq=numberhist  
             end if
             statdist(qq)=statdist(qq)+1.0_long
	  end if
       end do	  
              	      
       if (((ntime-1)/100000)<(ntime/100000)) then
           write(*,*) 'Time computed:',ntime,'out of 10000000'
       end if
  
 end do

 do ii=1,numberhist
    statdist(ii)=statdist(ii)/(real(timesteps)*h)
    xvalue(ii)=-1.0_long+(0.5_long+real(ii-1))*h
 end do	  
               
 open(unit=1, file='data_Figure9_1_b.dat', status='replace')

100 format(f10.6,x,f10.6)
 
do ii=1,numberhist
   write(1,100) xvalue(ii),statdist(ii)
end do
 
 close(unit=1)

end program Fortran_Figure9_1_b