Hey friends, everyday i show up here again to search for help because i’m still beggining to learn the language 
I don’t know if anyone of you is familiar to DSMC for gas dynamics but this is my code for now, but the thing is that it is stoping before the time it should, in the first ‘do’ and i don’t know why, if someone find what’s wrong or how to improve something please feel free to share. I think it is some parameter that is stopping the code because of some bug in the loop for changing the vellocities. N_mc is ny suspect but i’m not sure
program collisions
implicit none
integer :: t, c, i, j, k, N_col=0
integer, parameter :: N=200, dim = 3 !number of particles
integer, parameter :: Ncell=10 !number of cells
integer :: PPC(Ncell**3), cell_number, Nm !number of particles per cell
integer, parameter :: Nmft=20, Nt = 10000
real(8), parameter :: pi = 4.D0*DATAN(1.D0)
real(8), parameter :: radius = 0.06
real(8), parameter :: Tw = 1.0 !temperature of the wall
real(8), parameter :: sigmat = pi*radius**2 !total cross section
real(8), parameter :: Lx = 1.0 , Ly = 1.0 , Lz = 1.0, VOLUME = Lx*Ly*Lz !dimensions of the box
real(8), parameter :: n0 = N/VOLUME
real(8), parameter :: mfp = 1/(sqrt(2.0)*n0*sigmat**2) !mean free path
real(8), parameter :: Kn = mfp / Lx !Knudssen number
real(8), parameter :: v = 100.0, v_mean = sqrt(8*Tw/pi), v_rel_max = 200
real(8), parameter :: tau = mfp / v_mean, dt = 0.00001*tau
real(8), parameter :: dx = Lx/Ncell, dy=Ly/Ncell, dz = Lz/Ncell
real(8), parameter :: vol = dx*dy*dz, Nr = n0*VOLUME/N! !volume of each cell and number of real particles per cell
real(8) :: N_mc, Fn !Maximum number of model particle collisions and number of model particle per cell
real(8), dimension(N) :: x, y, z
real(8), dimension(N) :: vx, vy, vz
real(8), dimension(:), allocatable :: x_c, y_c, z_c, vx_c, vy_c, vz_c
integer, dimension(N) :: cell
integer, allocatable, dimension(:) :: particles_in_cell
real(8) :: v_rel, Rf1, Rf2, Rf3, r1, r2, vx_cm, vy_cm, vz_cm, g1, g2, g3, e, cosx, sinx
integer :: i_prop, j_prop
!setting initial positions
call random_number(x)
call random_number(y)
call random_number(z)
x = x * Lx
y = y * Ly
z = z * Lz
!setting initial velocities
where(x .le. Lx/2.0)
vx = v
else where
vx = -v
end where
vy = 0
vz = 0
do t=1, Nt
!drifting
x = x + vx*dt
y = y + vy*dt
z = z + vz*dt
!interacting with bounderies
where(x .le. 0)
vx = abs(vx)
endwhere
where(x .ge. Lx)
vx = -abs(vx)
endwhere
where(y .le. 0)
vy = abs(vy)
endwhere
where(y .ge. Ly)
vy = -abs(vy)
endwhere
where(z .le. 0)
vz = abs(vz)
endwhere
where(y .ge. Lz)
vz = -abs(vz)
endwhere
!indexing particles in cell
PPC = 0.0
do c=1, N
i = int(x(c) / dx)
j = int(y(c) / dy)
k = int(z(c) / dz)
cell_number = i + j * Ncell + k * Ncell**2
cell(c) = cell_number
PPC(cell_number) = PPC(cell_number) + 1
!print*, c, x(c), cell_number
enddo
!passing through all cells
do c=1, Ncell**3, 1
x_c = pack(x, cell == c)
y_c = pack(y, cell == c)
z_c = pack(z, cell == c)
!print*, c, PPC(c), x_c
vx_c = pack(vx, cell == c)
vy_c = pack(vy, cell == c)
vz_c = pack(vz, cell == c)
Nm = PPC(c) !number of number particles in the cell
Fn = Nr / Nm
N_mc = (Nm**2*Fn *sigmat * v_rel_max * dt ) / ( 2 * vol) !(Nm **2 *Fn * sigmat * v_rel_max * dt )
!print*, N_mc
do k=0, int(N_mc)
call random_number(Rf1)
!i and j prop must be integers in range of the number of particles in the cell Nm
call random_number(r1)
call random_number(r2)
i_prop = 1+floor(Nm*r1) !random integer in the interval [1, Nm]
j_prop = 1+floor(Nm*r2)
v_rel = sqrt((vx_c(i_prop) - vx_c(j_prop))**2 + (vy_c(i_prop) - vy_c(j_prop))**2 + (vz_c(i_prop) - vz_c(j_prop))**2)
!print*, v_rel
!accept-reject
if((v_rel/v_rel_max) > Rf1) then
vx_cm = 0.5 * (vx_c(i_prop) + vy_c(j_prop))
vy_cm = 0.5 * (vy_c(i_prop) + vy_c(j_prop))
vz_cm = 0.5 * (vz_c(i_prop) + vz_c(j_prop))
call random_number(Rf2)
call random_number(Rf3)
cosx = 2.0 * Rf2 - 1.0
sinx = sqrt(1 - cosx**2)
e = 2.0 * pi * Rf3
g1 = v_rel * cosx
g2 = v_rel * sinx * sin(e)
g3 = v_rel * sinx * cos(e)
vx_c(i_prop) = vx_cm + 0.5*g1
vy_c(i_prop) = vy_cm + 0.5*g2
vz_c(i_prop) = vz_cm + 0.5*g3
vx_c(j_prop) = vx_cm - 0.5*g1
vy_c(j_prop) = vy_cm - 0.5*g2
vz_c(j_prop) = vz_cm - 0.5*g3
N_col = N_col + 1
continue
endif
x = unpack(x_c, cell==c, x)
y = unpack(y_c, cell==c, y)
z = unpack(z_c, cell==c, z)
vx = unpack(vx_c, cell==c, vx)
vy = unpack(vy_c, cell==c, vy)
vz = unpack(vz_c, cell==c, vz)
deallocate(vx_c)
deallocate(vy_c)
deallocate(vz_c)
enddo
enddo
write(*,*)'T=',t, 'Number of collisions = ', N_col
enddo
end program collisions