program Fortran_Figure3_4b

implicit none

save

 integer, parameter :: long = selected_real_kind(15,50)
 integer, parameter :: ranpom=1000000
 integer, parameter :: noreal=10000
 integer, parameter :: nomesh=1951

 real(kind=long), parameter :: k1=0.1_long
 real(kind=long), parameter :: k2=1.0_long
 real(kind=long), dimension(ranpom) :: RanNumUnif 
 real(kind=long), dimension(19) :: ys,exittime  
 real(kind=long) :: A,time,a0
 real(kind=long), dimension(nomesh) :: y,tau,ps,intps 
 real(kind=long) :: dx
 integer :: ii, iru,et
 
 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(RanNumUnif)
 
 iru=0
 
 do et=1,19  
    exittime(et)=0.0_long
    ys(et)=real(et-1)    
    write(*,*) et,' out of 19'
    
    do ii=1,noreal
       time=0.0_long
       
       A=ys(et)
     
       do while (A<18.5_long)
          
          iru=iru+1
          if (iru>ranpom) then
             call random_number(RanNumUnif)
             iru=1
          end if 
          a0=k1*A+k2
          time=time+(1/a0)*log(1/RanNumUnif(iru))
          iru=iru+1
          if (iru>ranpom) then
             call random_number(RanNumUnif)
             iru=1
          end if 
          if (RanNumUnif(iru)*a0<k1*A) then
              A=A-1.0_long
          else
              A=A+1.0_long
          end if
       end do	     
       exittime(et)=exittime(et)+time
    end do
    exittime(et)=exittime(et)/real(noreal)
 end do
 
 open(unit=1, file='data_Figure3_4b_1.dat', status='replace')

100 format(i6,4X,e12.6)

 do et=1,19
    write(1,100) nint(ys(et)),exittime(et)
 end do

 close(unit=1)
 
 dx=19.5_long/real(nomesh-1)
 
 do ii=1,nomesh
    y(ii)=-0.5_long+dx*real(ii-1)
    ps(ii)=exp(-2_long*y(ii)+(4_long*k2/k1-1_long)*log(k1*y(ii)+k2))
 end do
 intps(1)=0.0_long
 do ii=2,nomesh
    intps(ii)=intps(ii-1)+dx*ps(ii-1)/2.0_long+dx*ps(ii)/2.0_long
 end do 
 tau(nomesh)=0.0_long
 do ii=1,nomesh-1
    tau(nomesh-ii)=tau(nomesh-ii+1)+dx*intps(nomesh-ii+1)/ps(nomesh-ii+1)/(k1*y(nomesh-ii+1)+k2) 
    tau(nomesh-ii)=tau(nomesh-ii)+dx*intps(nomesh-ii)/ps(nomesh-ii)/(k1*y(nomesh-ii)+k2) 
 end do   

 open(unit=2, file='data_Figure3_4b_2.dat', status='replace')

200 format(e12.6,X,e12.6)

 do ii=1,nomesh
    if (abs(10.0_long*y(ii)-nint(10.0_long*y(ii)))<dx*5.0_long) then
        write(2,200) y(ii),tau(ii)
    end if
 end do

 close(unit=2)
  
end program Fortran_Figure3_4b