Added explicit array range for a bunch of whole array operations. Als…

…o fixed a lingering bug in RMVS that was causing the planet accelerations to be miscalculated after an encounter was done if oblateness was turned on

src/discard/discard.f90

@@ -53,7 +53,7 @@ module subroutine discard_pl(self, system, param)
       class(swiftest_parameters),   intent(inout) :: param  !! Current run configuration parameter
 
       if (self%nbody == 0) return
-      self%ldiscard(:) = .false.
+      self%ldiscard(1:self%nbody) = .false.
 
       return
    end subroutine discard_pl
@@ -86,9 +86,9 @@ module subroutine discard_tp(self, system, param)
          if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or.  (param%rmaxu >= 0.0_DP)) call discard_cb_tp(tp, system, param)
          if (param%qmin >= 0.0_DP) call discard_peri_tp(tp, system, param)
          if (param%lclose) call discard_pl_tp(tp, system, param)
-         if (any(tp%ldiscard)) then
+         if (any(tp%ldiscard(1:ntp))) then
             allocate(ldiscard, source=tp%ldiscard)
-            call tp%spill(system%tp_discards, ldiscard, ldestructive=.true.)
+            call tp%spill(system%tp_discards, ldiscard(1:ntp), ldestructive=.true.)
          end if
       end associate
 

src/helio/helio_kick.f90

@@ -65,9 +65,9 @@ module subroutine helio_kick_getacch_tp(self, system, param, t, lbeg)
       associate(tp => self, cb => system%cb, pl => system%pl, npl => system%pl%nbody)
          system%lbeg = lbeg
          if (system%lbeg) then
-            call tp%accel_int(pl%Gmass(:), pl%xbeg(:,:), npl)
+            call tp%accel_int(pl%Gmass(1:npl), pl%xbeg(:,1:npl), npl)
          else
-            call tp%accel_int(pl%Gmass(:), pl%xend(:,:), npl)
+            call tp%accel_int(pl%Gmass(1:npl), pl%xend(:,1:npl), npl)
          end if
          if (param%loblatecb) call tp%accel_obl(system)
          if (param%lextra_force) call tp%accel_user(system, param, t, lbeg)
@@ -99,7 +99,7 @@ module subroutine helio_kick_vb_pl(self, system, param, t, dt, lbeg)
       if (self%nbody == 0) return
 
       associate(pl => self, npl => self%nbody)
-         pl%ah(:,:) = 0.0_DP
+         pl%ah(:, 1:npl) = 0.0_DP
          call pl%accel(system, param, t, lbeg)
          if (lbeg) then
             call pl%set_beg_end(xbeg = pl%xh)
@@ -136,7 +136,7 @@ module subroutine helio_kick_vb_tp(self, system, param, t, dt, lbeg)
       if (self%nbody == 0) return
 
       associate(tp => self, ntp => self%nbody)
-         tp%ah(:,:) = 0.0_DP
+         tp%ah(:, 1:ntp) = 0.0_DP
          call tp%accel(system, param, t, lbeg)
          do concurrent(i = 1:ntp, tp%lmask(i)) 
             tp%vb(:, i) = tp%vb(:, i) + tp%ah(:, i) * dt

src/kick/kick.f90

@@ -64,6 +64,8 @@ module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl)
       real(DP)                  :: rji2, irij3, fac, r2
       real(DP), dimension(NDIM) :: dx
 
+      if ((self%nbody == 0) .or. (npl == 0)) return
+
       associate(tp => self, ntp => self%nbody)
          do concurrent(i = 1:ntp, tp%lmask(i))
             do j = 1, npl

src/rmvs/rmvs_encounter_check.f90

@@ -23,13 +23,14 @@ module function rmvs_encounter_check_tp(self, system, dt) result(lencounter)
       real(DP), dimension(system%pl%nbody)    :: r2crit
       logical                                 :: lflag
 
+      lencounter = .false.
       if (self%nbody == 0) return
 
       select type(pl => system%pl)
       class is (rmvs_pl)
          associate(tp => self, ntp => self%nbody, npl => pl%nbody, rts => system%rts)
-            r2crit(:) = (rts * pl%rhill(:))**2
-            tp%plencP(:) = 0
+            r2crit(1:npl) = (rts * pl%rhill(1:npl))**2
+            tp%plencP(1:ntp) = 0
             do j = 1, npl
                do i = 1, ntp
                   if ((.not.tp%lmask(i)).or.(tp%plencP(i) /= 0)) cycle
@@ -41,9 +42,9 @@ module function rmvs_encounter_check_tp(self, system, dt) result(lencounter)
                   lflag = rmvs_chk_ind(r2, v2, vdotr, dt, r2crit(j))
                   if (lflag) tp%plencP(i) = j
                end do
-               pl%nenc(j) = count(tp%plencP(:) == j)
+               pl%nenc(j) = count(tp%plencP(1:ntp) == j)
             end do
-            lencounter = any(pl%nenc(:) > 0)
+            lencounter = any(pl%nenc(1:npl) > 0)
          end associate
       end select
       return

src/rmvs/rmvs_kick.f90

@@ -61,7 +61,9 @@ module subroutine rmvs_kick_getacch_tp(self, system, param, t, lbeg)
                         end if
 
                         ! Swap the planetocentric and heliocentric position vectors and central body masses
-                        tp%xh(:,:) = tp%xheliocentric(:,:)
+                        do concurrent(i = 1:ntp, tp%lmask(i))
+                           tp%xh(:, i) = tp%xheliocentric(:, i)
+                        end do
                         GMcb_original = cb%Gmass
                         cb%Gmass = tp%cb_heliocentric%Gmass
 
@@ -71,7 +73,7 @@ module subroutine rmvs_kick_getacch_tp(self, system, param, t, lbeg)
                         if (param%lgr) call tp%accel_gr(param)
 
                         ! Put everything back the way we found it
-                        tp%xh(:,:) = xh_original(:,:)
+                        call move_alloc(xh_original, tp%xh)
                         cb%Gmass = GMcb_original
 
                      end associate