program Fortran_Figure6_2

implicit none

save

 integer, parameter :: long = selected_real_kind(15,50)
 integer, parameter :: totalmol=100000  ! number of molecules
 integer, parameter :: ranpom=1000000
 integer, parameter :: ranpom2=500000
 integer, parameter :: noboxes=40
 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), dimension(2,ranpom2) :: RanNum 
 real(kind=long), dimension(ranpom) :: RanNumNorm 
 real(kind=long), dimension(ranpom) :: RanNumUnif 
 
 integer, dimension(noboxes,2) :: histsave 
 real(kind=long), dimension(totalmol) :: zpos
 logical, dimension(totalmol) :: inuse
 integer, dimension(totalmol) :: indexempty
 
 integer :: ww,jj,zz,ss,irn,iru,ll,noadd,lastmol,numbermol,emptyuse,ns
 real(kind=long) :: dt,pom1,pom2,time,alpha,oldposZ,prob1,c0,probaddnext
 real(kind=long) :: creationrate,savingtime,finaltime
 logical :: savingdata
 
 integer size, seed(2)
 data seed /2, 5/ 
      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
 
 alpha=sqrt(2.0_long*DiffConst)
 
 lastmol=1
 numbermol=1
 
 do jj=1,lastmol
    zpos(jj)=1.5_long
    inuse(jj)=.TRUE.
 end do
 do jj=lastmol+1,totalmol
    zpos(jj)=0.0_long
    inuse(jj)=.FALSE.
 end do
  
 do ns=1,2
    do ss=1,noboxes
       histsave(ss,ns)=0
    end do
 end do
      
 time=0.0_long
 dt=0.01_long
 emptyuse=0       
 creationrate=k2*dt*(0.2_long/0.025_long)
 savingtime=10.0_long*60.0_long
 finaltime=30.0_long*60.0_long
 savingdata=.TRUE.
     
 do while (time<finaltime)

       time=time+dt
       
       ! evolve every molecule over one time step
       
       do jj=1,totalmol 
          
	  if (inuse(jj)) then             
       
	      ! evolve positions of molecules over one time step
	      
	      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 
	      oldposZ=zpos(jj)   
              zpos(jj)=zpos(jj)+alpha*RanNumNorm(irn)*sqrt(dt)
	      
	      ! boundary condition at z=1 - no flux
	      
              if (zpos(jj)>1.0_long) then
	          zpos(jj)=2.0_long-zpos(jj)
              end if
	      
	      ! boundary condition at z=0 - no flux
	      	      	      
              if (zpos(jj)<0.0_long) then
	          zpos(jj)=-zpos(jj)
              end if
	      
	      ! degradation with rate k1
	  
	      iru=iru+1
              if (iru>ranpom) then
                  call random_number(RanNumUnif)
                  iru=1
              end if    
              if (RanNumUnif(iru)<k1*dt) then
	          inuse(jj)=.FALSE.
		  emptyuse=emptyuse+1
		  numbermol=numbermol-1
		  indexempty(emptyuse)=jj
              end if
	      
	  end if    
       end do
              	      
       ! creation of molecules in x=[0,0.2]
       
       iru=iru+1
       if (iru>ranpom) then
           call random_number(RanNumUnif)
           iru=1
       end if    
       if (RanNumUnif(iru)<creationrate) then
 	   if (emptyuse>0) then
               iru=iru+1
               if (iru>ranpom) then
                   call random_number(RanNumUnif)
                   iru=1
               end if    
	       zpos(indexempty(emptyuse))=0.2_long*RanNumUnif(iru)
               inuse(indexempty(emptyuse))=.TRUE.
               numbermol=numbermol+1
	       emptyuse=emptyuse-1
           else
	       lastmol=lastmol+1
               iru=iru+1
               if (iru>ranpom) then
                   call random_number(RanNumUnif)
                   iru=1
               end if    
	       zpos(lastmol)=0.2_long*RanNumUnif(iru)
               inuse(lastmol)=.TRUE.
               numbermol=numbermol+1
	       if (lastmol>totalmol) then
	           write(*,*) 'Not enough memory. Increase totalmol and compute again!'
		   STOP
	       end if	   
	   end if	  
       end if      
	                    	      
    if (int(time*0.01_long)-int(time*0.01_long-dt*0.01_long)>0.5_long) then
        write(*,*) 'Time computed:',nint(time),'seconds (final time 1800 sec)'
    end if
    
    if ((time>savingtime).AND.savingdata) then
        ns=1
	savingdata=.FALSE.
        do jj=1,lastmol
           if (inuse(jj)) then
	       ss=nint(real(noboxes)*zpos(jj)+0.5000_long)
	       if (ss>noboxes) then
	           ss=noboxes
	       end if	   
	       if (ss<1) then
	           ss=1
	       end if	   
	       histsave(ss,ns)=histsave(ss,ns)+1
           end if	   
        end do   	         
    end if
    
 end do
    
 ns=2
 do jj=1,lastmol
    if (inuse(jj)) then
	ss=nint(real(noboxes)*zpos(jj)+0.5000_long)
	if (ss>noboxes) then
	    ss=noboxes
	end if	   
	if (ss<1) then
	    ss=1
	end if	   
	histsave(ss,ns)=histsave(ss,ns)+1
    end if	   
 end do   	         
     
  
 open(unit=1, file='data_Figure6_2.dat', status='replace')

100 format(i6,X,i6)
 
       do ss=1,noboxes
          write(1,100) histsave(ss,1),histsave(ss,2)
       end do   	    

 close(unit=1)

end program Fortran_Figure6_2
