program  Fortran_Figure6_6_b

implicit none

save

 integer, parameter :: long = selected_real_kind(15,50)
 integer, parameter :: totalmol=10000  ! 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 :: sigma=0.002_long ! rate of adsorbtion
 real(kind=long), parameter :: pdes=0.005_long ! rate of desorption
 
 real(kind=long), dimension(2,ranpom2) :: RanNum 
 real(kind=long), dimension(ranpom) :: RanNumNorm 
 real(kind=long), dimension(ranpom) :: RanNumUnif 
 
 integer, dimension(noboxes) :: histsave 
 real(kind=long), dimension(totalmol) :: zpos
 logical, dimension(totalmol) :: inuse
 
 integer :: ww,jj,zz,ss,irn,iru,ll,noadd
 real(kind=long) :: dt,pom1,pom2,time,alpha,oldposZ,prob1,c0,probaddnext,pads
 
 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
 
 alpha=sqrt(2.0_long*DiffConst)
 
 do jj=1,totalmol
    iru=iru+1
    if (iru>ranpom) then
        call random_number(RanNumUnif)
        iru=1
    end if    
    zpos(jj)=RanNumUnif(iru)*0.5_long
    inuse(jj)=.TRUE.
 end do
  
 do ss=1,noboxes
    histsave(ss)=0
 end do
      
 time=0.0_long
 dt=0.1_long
 pads=sigma*sqrt(3.14159265358979_long*dt/DiffConst)/2.0_long
      
    do while (time<600.0_long)

       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 
	      !              - adsorbing with probability pads
	      !              - reflecting otherwise
	      	      
	      
              if (zpos(jj)<0.0_long) then
	          iru=iru+1
                  if (iru>ranpom) then
                      call random_number(RanNumUnif)
                      iru=1
                  end if    
                  if (RanNumUnif(iru)<pads) then
	              inuse(jj)=.FALSE.
		  else
	              zpos(jj)=-zpos(jj)
		  end if    
             else
	          iru=iru+1
                  if (iru>ranpom) then
                      call random_number(RanNumUnif)
                      iru=1
                  end if
		  prob1=exp(-2*zpos(jj)*oldposZ/(2.0_long*DiffConst*dt))    
                  if (RanNumUnif(iru)<prob1) then
	              iru=iru+1
                      if (iru>ranpom) then
                          call random_number(RanNumUnif)
                          iru=1
                      end if    
                      if (RanNumUnif(iru)<pads) then
	                  inuse(jj)=.FALSE.
		      end if    
                  end if
              end if
	      
	  else
	  
	  ! molecules is adsorbed, there is a finite probability of desorption
	  
	      iru=iru+1
              if (iru>ranpom) then
                 call random_number(RanNumUnif)
                 iru=1
              end if    
              if (RanNumUnif(iru)<pdes*dt) then
	         inuse(jj)=.TRUE.
	         zpos(jj)=0.0_long
              end if    
	  
	  
	  end if    
       end do
	                    	      
    end do
    
    do jj=1,totalmol
       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)=histsave(ss)+1
       end if	   
    end do   	         
     
  
 open(unit=1, file='data_Figure6_6_b.dat', status='replace')

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

 close(unit=1)


 end program Fortran_Figure6_6_b