defined all values of c_lm less than epsilon as 0.0_DP to reduce extr…

…a calculations

src/swiftest/swiftest_io.f90

@@ -3353,6 +3353,13 @@ module subroutine swiftest_io_read_in_system(self, nc, param)
 
       param%lshgrav = allocated(self%cb%c_lm) !! .and. (size(self%cb%c_lm) /= 0) 
 
+      if(param%lshgrav) then
+         ! Replace elements of c_lm smaller than epsilon with 0.0
+         WHERE (abs(self%cb%c_lm) < EPSILON(0.0_DP))
+            self%cb%c_lm = 0.0_DP
+         END WHERE
+      end if
+
       param%loblatecb = ((self%cb%j2rp2 /= 0.0_DP) .or. (self%cb%j4rp4 /= 0.0_DP)) .and. (.not. param%lshgrav)
       if (.not.param%loblatecb .and. .not.param%lshgrav) then
          if (allocated(self%pl%aobl)) deallocate(self%pl%aobl)

src/swiftest/swiftest_sph.f90

@@ -185,7 +185,6 @@ module subroutine swiftest_sph_g_acc_pl_all(self, nbody_system)
         class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
         ! Internals
         integer(I4B)    :: i = 1
-        real(DP)        :: theta, phi              !! zenith angle, and azimuthal angle
         real(DP), dimension(NDIM)  :: g_sph        !! Gravitational terms from Spherical Harmonics
 
         associate(pl => self, npl => self%nbody, cb => nbody_system%cb, rh => self%rh)
@@ -216,8 +215,6 @@ module subroutine swiftest_sph_g_acc_tp_all(self, nbody_system)
         class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object
         ! Internals
         integer(I4B)    :: i = 1
-        real(DP)        :: theta, phi              !! zenith angle, and azimuthal angle
-        real(DP), dimension(NDIM)  :: rh           !! Position vector of the test particle
         real(DP), dimension(NDIM)  :: g_sph        !! Gravitational terms from Spherical Harmonics
         real(DP), dimension(NDIM)  :: aoblcb       !! Temporary variable for central body oblateness acceleration