program Fortran_Figure3_3b

implicit none

save

 integer, parameter :: long = selected_real_kind(15,50)
 integer, parameter :: ranpom=1000000
 integer, parameter :: ranpom2=500000
 integer, parameter :: noreal=1000
 integer, parameter :: nomesh=10001

 real(kind=long), parameter :: k1=0.001_long
 real(kind=long), parameter :: k2=0.75_long
 real(kind=long), parameter :: k3=165.0_long
 real(kind=long), parameter :: k4=10000.0_long
 real(kind=long), parameter :: k5=200.0_long
 real(kind=long), dimension(2,ranpom2) :: RanNum 
 real(kind=long), dimension(ranpom) :: RanNumNorm 
 real(kind=long), dimension(ranpom) :: RanNumUnif 
 real(kind=long), dimension(249) :: exittime 
 real(kind=long) :: X,time,pom1,pom2,dt,oldX,dx
 real(kind=long), dimension(nomesh) :: y,tau,ps,intps
 integer :: ii,zz,irn,iru,et
 logical :: insidedomain
 
 integer size, seed(2)
 data seed /10, 10/ 
      size = 2

 call random_seed(SIZE=size)
 call random_seed(PUT=seed(1:size))
 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=0

 call random_seed(SIZE=size)
 call random_seed(PUT=seed(1:size))
 call random_number(RanNumUnif)
 
 iru=0
 dt=0.001_long
 
 do et=1,249  
    exittime(et)=0.0_long    
    write(*,*) et,' out of 249'
    
    do ii=1,noreal
       time=0.0_long
       X=real(et)
       insidedomain=.TRUE.
    
       do while (insidedomain)
          
          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 
       
          oldX=X
          X=X+dt*(-k1*X*X*X+k2*X*X-k3*X+k4)+k5*sqrt(dt)*RanNumNorm(irn)
       
          if (X>250.0_long) then
             insidedomain=.FALSE.
          else
             iru=iru+1
             if (iru>ranpom) then
                call random_number(RanNumUnif)
                iru=1
             end if 
             if (RanNumUnif(iru)<exp(-2*(250.0_long-X)*(250.0_long-oldX)/(k5*k5*dt))) then	  
                insidedomain=.FALSE.
	     end if    
          end if
          time=time+dt
       end do
       exittime(et)=exittime(et)+time
    end do      
    exittime(et)=exittime(et)/real(noreal) 
 end do
 
 open(unit=1, file='data_Figure3_3b_1.dat', status='replace')

100 format(i6,X,e12.6)
 
 do et=1,249
    write(1,100) et,exittime(et)
 end do

 close(unit=1)

 dx=250.0_long/real(nomesh-1)
 do ii=1,nomesh
    y(ii)=250.0_long*real(ii-1)/real(nomesh-1)
    ps(ii)=exp((-3*k1*y(ii)*y(ii)*y(ii)*y(ii)+4*k2*y(ii)*y(ii)*y(ii)-6*k3*y(ii)*y(ii)+12*k4*y(ii))/(6*k5*k5))
 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)/(k5*k5) 
    tau(nomesh-ii)=tau(nomesh-ii)+dx*intps(nomesh-ii)/ps(nomesh-ii)/(k5*k5) 
 end do   

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

200 format(e12.6,X,e12.6)
 
 do ii=1,nomesh
    if (abs(y(ii)-nint(y(ii)))<dx/2.0_long) then
        write(2,200) y(ii),tau(ii)
    end if
 end do

 close(unit=2)
  
  
end program Fortran_Figure3_3b