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.