Merge branch 'debug'

src/io/io.f90

@@ -172,6 +172,7 @@ module subroutine io_dump_system(self, param)
       real(DP)                      :: deltawall, wallperstep, tfrac
       integer(I8B)                  :: clock_count, count_rate, count_max
       character(*),     parameter   :: statusfmt   = '("Time = ", ES12.5, "; fraction done = ", F6.3, "; Number of active pl, tp = ", I5, ", ", I5)'
+      character(*),     parameter   :: symbastatfmt   = '("Time = ", ES12.5, "; fraction done = ", F6.3, "; Number of active plm, pl, tp = ", I5, ", ", I5, ", ", I5)'
       character(len=*), parameter   :: walltimefmt = '("      Wall time (s): ", es12.5, "; Wall time/step in this interval (s):  ", es12.5)'
       logical, save                 :: lfirst = .true.
       real(DP), save                :: start, finish
@@ -209,7 +210,12 @@ module subroutine io_dump_system(self, param)
          deltawall = clock_count / (count_rate * 1.0_DP) - finish
          wallperstep = deltawall / param%istep_dump
          finish = clock_count / (count_rate * 1.0_DP)
-         write(*, statusfmt) param%t, tfrac, self%pl%nbody, self%tp%nbody
+         select type(pl => self%pl)
+         class is (symba_pl)
+            write(*, symbastatfmt) param%t, tfrac, pl%nplm, pl%nbody, self%tp%nbody
+         class default
+            write(*, statusfmt) param%t, tfrac, pl%nbody, self%tp%nbody
+         end select
          write(*, walltimefmt) finish - start, wallperstep
          if (param%lenergy) call self%conservation_report(param, lterminal=.true.)
       end if

src/symba/symba_collision.f90

@@ -220,7 +220,7 @@ module function symba_collision_casemerge(system, param, family, x, v, mass, rad
       ! Result
       integer(I4B)                                   :: status           !! Status flag assigned to this outcome
       ! Internals
-      integer(I4B)                              :: i, j, ibiggest, nfamily
+      integer(I4B)                              :: i, j, k, ibiggest, nfamily
       real(DP)                                  :: volume_new, pe
       real(DP), dimension(NDIM)                 :: xc, vc, xcrossv
       real(DP), dimension(2)                    :: vol
@@ -283,8 +283,27 @@ module function symba_collision_casemerge(system, param, family, x, v, mass, rad
          param%Ecollisions = param%Ecollisions + pe 
          param%Euntracked  = param%Euntracked - pe 
 
+         ! Update any encounter lists that have the removedbodies in them so that they instead point to the new 
+         do k = 1, system%plplenc_list%nenc
+            do j = 1, nfamily
+               i = family(j)
+               if (i == ibiggest) cycle
+               if (system%plplenc_list%id1(k) == pl%id(i)) then
+                  system%plplenc_list%id1(k) = pl%id(ibiggest)
+                  system%plplenc_list%index1(k) = i
+               end if
+               if (system%plplenc_list%id2(k) == pl%id(i)) then
+                  system%plplenc_list%id2(k) = pl%id(ibiggest)
+                  system%plplenc_list%index2(k) = i
+               end if
+               if (system%plplenc_list%id1(k) == system%plplenc_list%id2(k)) system%plplenc_list%status(k) = INACTIVE
+            end do
+         end do
+
+
          status = MERGED
          call symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, m_frag, rad_frag, xb_frag, vb_frag, rot_frag, status)
+
 
       end select
 

src/symba/symba_util.f90

@@ -468,6 +468,7 @@ module subroutine symba_util_rearray_pl(self, system, param)
       logical, dimension(:), allocatable :: lmask
       class(symba_plplenc), allocatable :: plplenc_old
       logical :: lencounter
+      integer(I4B), dimension(:), allocatable :: levelg_orig_pl, levelm_orig_pl, levelg_orig_tp, levelm_orig_tp, nplenc_orig_pl, nplenc_orig_tp, ntpenc_orig_pl
 
       associate(pl => self, pl_adds => system%pl_adds, nadd => system%pl_adds%nbody)
 
@@ -522,12 +523,27 @@ module subroutine symba_util_rearray_pl(self, system, param)
          pl%kin(1:npl)%nchild = 0
          pl%kin(1:npl)%parent = [(i, i=1, npl)]
 
-         ! Re-build the encounter list
+         ! Re-build the zero-level encounter list, being sure to save the original level information for all bodies
+         allocate(levelg_orig_pl, source=pl%levelg)
+         allocate(levelm_orig_pl, source=pl%levelm)
+         allocate(nplenc_orig_pl, source=pl%nplenc)
+         allocate(ntpenc_orig_pl, source=pl%ntpenc)
          lencounter = pl%encounter_check(system, param%dt, 0) 
-         select type(tp => system%tp)
-         class is (symba_tp)
-            lencounter = tp%encounter_check(system, param%dt, 0)
-         end select
+         if (system%tp%nbody > 0) then
+            select type(tp => system%tp)
+            class is (symba_tp)
+               allocate(levelg_orig_tp, source=tp%levelg)
+               allocate(levelm_orig_tp, source=tp%levelm)
+               allocate(nplenc_orig_tp, source=tp%nplenc)
+               lencounter = tp%encounter_check(system, param%dt, 0)
+               call move_alloc(levelg_orig_tp, tp%levelg)
+               call move_alloc(levelm_orig_tp, tp%levelm)
+               call move_alloc(nplenc_orig_tp, tp%nplenc)
+            end select
+         end if
+         call move_alloc(levelg_orig_pl, pl%levelg)
+         call move_alloc(levelm_orig_pl, pl%levelm)
+         call move_alloc(nplenc_orig_pl, pl%nplenc)
 
          associate(idnew1 => system%plplenc_list%id1, idnew2 => system%plplenc_list%id2, idold1 => plplenc_old%id1, idold2 => plplenc_old%id2)
             do k = 1, system%plplenc_list%nenc