Fixed calcululation of density that was incorrect

src/collision/collision_resolve.f90

@@ -321,13 +321,14 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status)
       class(swiftest_pl), allocatable           :: plnew, plsub
       character(*), parameter :: FRAGFMT = '("Newbody",I0.7)'
       character(len=NAMELEN) :: newname, origin_type
+      real(DP) :: volume
 
       select type(nbody_system)
       class is (swiftest_nbody_system)
       select type(param)
       class is (swiftest_parameters)
          associate(pl => nbody_system%pl, pl_discards => nbody_system%pl_discards, info => nbody_system%pl%info, pl_adds => nbody_system%pl_adds, cb => nbody_system%cb, npl => pl%nbody, &
-            collision_basic => nbody_system%collider, impactors => nbody_system%collider%impactors,fragments => nbody_system%collider%fragments)
+            collider => nbody_system%collider, impactors => nbody_system%collider%impactors,fragments => nbody_system%collider%fragments)
 
             ! Add the impactors%id bodies to the subtraction list
             nimpactors = impactors%ncoll
@@ -351,7 +352,10 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status)
             plnew%mass(1:nfrag) = fragments%mass(1:nfrag)
             plnew%Gmass(1:nfrag) = param%GU * fragments%mass(1:nfrag)
             plnew%radius(1:nfrag) = fragments%radius(1:nfrag)
-            plnew%density(1:nfrag) = fragments%mass(1:nfrag) / fragments%radius(1:nfrag)
+            do concurrent(i = 1:nfrag)
+               volume = 4.0_DP/3.0_DP * PI * plnew%radius(i)**3
+               plnew%density(i) = fragments%mass(i) / volume
+            end do
             call plnew%set_rhill(cb)
 
             select case(status)
@@ -437,7 +441,7 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status)
             end where
 
             ! Log the properties of the new bodies
-            select type(after => collision_basic%after)
+            select type(after => collider%after)
             class is (swiftest_nbody_system)
                allocate(after%pl, source=plnew)
             end select

src/swiftest/swiftest_util.f90

@@ -1660,7 +1660,7 @@ module subroutine swiftest_util_rearray_pl(self, nbody_system, param)
       ! Internals
       class(swiftest_pl), allocatable :: tmp !! The discarded body list.
       integer(I4B) :: i, k, npl, nadd, nencmin, nenc_old, idnew1, idnew2, idold1, idold2
-      logical, dimension(:), allocatable :: lmask, ldump_mask
+      logical, dimension(:), allocatable :: lmask
       class(encounter_list), allocatable :: plplenc_old
       logical :: lencounter
 
@@ -1696,14 +1696,7 @@ module subroutine swiftest_util_rearray_pl(self, nbody_system, param)
          if (nadd > 0) then
             ! Append the adds to the main pl object
             call pl%append(pl_adds, lsource_mask=[(.true., i=1, nadd)])
-
-            allocate(ldump_mask(npl+nadd))  ! This mask is used only to append the original Fortran binary particle.dat file with new bodies. This is ignored for NetCDF output
-            ldump_mask(1:npl) = .false.
-            ldump_mask(npl+1:npl+nadd) = pl%status(npl+1:npl+nadd) == NEW_PARTICLE
             npl = pl%nbody
-         else
-            allocate(ldump_mask(npl))
-            ldump_mask(:) = .false.
          end if
 
          ! Reset all of the status flags for this body
@@ -1727,7 +1720,6 @@ module subroutine swiftest_util_rearray_pl(self, nbody_system, param)
          call nc%open(param)
          call pl%write_info(nc, param)
          call nc%close()
-         deallocate(ldump_mask)
 
          ! Reindex the new list of bodies 
          call pl%sort("mass", ascending=.false.)