external_force_with_gravity Subroutine

    subroutine external_force_with_gravity(self, particle, pairs, acc, n)
        class(external_force_type), intent(inout) :: self
        type(particle_type), intent(in) :: particle     !! 区域
        type(nnps_pairs_type), intent(in) :: pairs      !! 粒子对
        real(rk), intent(inout) :: acc(:, :)            !! 加速度 \( F_{external} = -g + F_{boundary} \)
        integer, intent(in) :: n                        !! 实粒子数 \( N_{real} \)
        integer :: i

        ! -------------------------------- 计算重力 --------------------------------- !
        acc(2, :n) = acc(2, :n) - g   ! 仅对实体粒子起作用，虚粒子不需要求解加速度而运动

        ! ----------------------------- 计算边界排斥力 ------------------------------- !
        do i = 1, pairs%len
            associate (ip => pairs%items(2*i - 1), jp => pairs%items(2*i), r => pairs%rdx(3*i - 2))
                select case (pairs%contact_type(i))
                case (3:4)
                    if (r < self%r0) then
                        associate (f => ((self%r0/r)**p1 - (self%r0/r)**p2)/(r*r), &
                                   dx => pairs%rdx(3*i - 1:3*i))
                            acc(:, ip) = acc(:, ip) + self%dd*dx(:)*f
                            acc(:, jp) = acc(:, jp) - self%dd*dx(:)*f
                        end associate
                    end if
                end select
            end associate
        end do

    end subroutine external_force_with_gravity