Fixed several issues involved in calculating oblateness acceleration …

…terms for test particles undergoing encounters

src/obl/obl.f90

@@ -28,7 +28,6 @@
             aobl(:, i) = fac1 * xh(:, i)
             aobl(3, i) = fac2 * xh(3, i) + aobl(3, i)
          end do
-         aoblcb(:) = 0.0_DP
          select type(self)
          class is (swiftest_pl)
             associate(Mpl => self%Gmass)

src/rmvs/rmvs_setup.f90

@@ -30,11 +30,16 @@ module subroutine rmvs_setup_pl(self,n)
             do i = 0, NTENC
                allocate(pl%outer(i)%x(NDIM, n))
                allocate(pl%outer(i)%v(NDIM, n))
+               pl%outer(i)%x(:,:) = 0.0_DP
+               pl%outer(i)%v(:,:) = 0.0_DP
             end do
             do i = 0, NTPHENC
                allocate(pl%inner(i)%x(NDIM, n))
                allocate(pl%inner(i)%v(NDIM, n))
                allocate(pl%inner(i)%aobl(NDIM, n))
+               pl%inner(i)%x(:,:) = 0.0_DP
+               pl%inner(i)%v(:,:) = 0.0_DP
+               pl%inner(i)%aobl(:,:) = 0.0_DP
             end do
          end if
       end associate

src/rmvs/rmvs_step.f90

@@ -278,8 +278,14 @@ subroutine rmvs_interp_in(pl, cb, dt, outer_index, config)
                pl%inner(inner_index)%aobl(:, :) = pl%aobl(:, :) 
             end if
          end do
-         ! Put the planet positions back into place
-         if (config%loblatecb) call move_alloc(xh_original, pl%xh)
+         if (config%loblatecb) then
+            ! Calculate the final value of oblateness accelerations at the final inner substep
+            pl%xh(:,:) = pl%inner(NTPHENC)%x(:, :)
+            call pl%obl_acc(cb)
+            pl%inner(NTPHENC)%aobl(:, :) = pl%aobl(:, :) 
+            ! Put the planet positions back into place
+            call move_alloc(xh_original, pl%xh)
+         end if
       end associate
       return