program Fortran_Figure3_5b

implicit none

save

 integer, parameter :: long = selected_real_kind(15,50)
 integer, parameter :: ranpom=1000000
 integer, parameter :: noreal=10000
 
 real(kind=long), parameter :: k1=0.00025_long
 real(kind=long), parameter :: k2=0.18_long
 real(kind=long), parameter :: k3=37.5_long
 real(kind=long), parameter :: k4=2200.0_long
 real(kind=long), dimension(2,ranpom) :: RanNumUnif 
 real(kind=long), dimension(235) :: exittimes ! from [0,235] for X=0,1,...,234
 integer, dimension(235) :: initialx
 real(kind=long) :: X,ss,a0,time,tau
 integer :: ii,iru,jj
 
 integer size, seed(2)
 data seed /100, 100/   !(you can choose different values here)
      size = 2

 call random_seed(SIZE=size)
 call random_seed(PUT=seed(1:size))
 call random_number(RanNumUnif)
 
 iru=0
   
 do ii=1,235
    
    initialx(ii)=ii-1
    exittimes(ii)=0.0_long
    write(*,*) ii,' out of 235'
        
    do jj=1,noreal
 
       time=0.0_long
       X=real(initialx(ii))
   
       do while (X<234.5_long)
          
          a0=k1*X*(X-1)*(X-2)+k2*X*(X-1)+k3*X+k4
          iru=iru+1
          if (iru>ranpom) then
             call random_number(RanNumUnif)
	     iru=1
          endif
          tau=(1/a0)*log(1/RanNumUnif(1,iru))
          time=time+tau
          ss=k2*X*(X-1)
          if (ss>RanNumUnif(2,iru)*a0) then
             X=X+1
          else
             ss=ss+k1*X*(X-1)*(X-2)
             if (ss>RanNumUnif(2,iru)*a0) then
                X=X-1
             else
                ss=ss+k3*X     
                if (ss>RanNumUnif(2,iru)*a0) then
                   X=X-1
                else
                   ss=ss+k4
                   if (ss>RanNumUnif(2,iru)*a0) then
                      X=X+1
                   else
                      STOP
                   endif
                endif    
             endif       
          endif
       end do 
       exittimes(ii)=exittimes(ii)+time
    end do 
    exittimes(ii)=exittimes(ii)/real(noreal)
 end do  
 
 open(unit=1, file='data_Figure3_5b.dat', status='replace')

100 format(i6,3x,e14.8)
 
 do ii=1,235
    write(1,100) initialx(ii),exittimes(ii)
 end do

 close(unit=1)
 
end program Fortran_Figure3_5b