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
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
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
Note that as each
particle_ad array is now
allocatable, you will need to initialise each array before you can call
CLICK HERE to find out more related problems solutions.