Your `particle_ad`

arrays only have space for `10*L`

particles, but you appear to be trying to store up to `N=2*L**2`

particles in them (depending on how the random numbers fall). On average, each will have enough space, but your code will fail if too many particles fall into a single block, and (if I’m remembering the probabilities right) the chances of this happening increase as `L`

increases.

You could solve this problem by replacing `Integer :: particle_ad(10*L)`

with `Integer :: particle_ad(N)`

, but this will waste an awful lot of space.

A better solution would be to re-size the `particle_ad`

arrays on the fly, making them bigger every time they get full. For convenience, you could wrap this behaviour in a method of the `block_p`

class.

For example,

```
type block_p
Integer :: partical_N
Integer, allocatable :: particle_ad(:)
contains
subroutine add_particle(this, index)
class(block_p), intent(inout) :: this
integer, intent(in) :: index
integer, allocatable :: temp
! Update partical_N.
this%partical_N = this%partical_N + 1
! Resize particle_ad if it is full.
if (size(this%particle_ad)<this%partical_N) then
temp = this%particle_ad
deallocate(this%particle_ad)
allocate(this%particle_ad(2*this%partical_N))
this%particle_ad(:size(temp)) = temp
endif
! Add the new index to particle_ad.
this%particle_ad(this%partical_N) = index
end subroutine
end type
```

You could then replace the lines

```
C(I_b,J_b)%partical_N = C(I_b,J_b)%partical_N + 1
C(I_b,J_b)%particle_ad( C(I_b,J_b)%partical_N ) = i
```

with

```
call C(I_b,J_b)%add_particle(i)
```

Note that as each `particle_ad`

array is now `allocatable`

, you will need to initialise each array before you can call `add_particle`

.

CLICK HERE to find out more related problems solutions.