Made somewhat functional SyMBA version of the discard method for pl

src/io/io.f90

@@ -1074,7 +1074,7 @@ module subroutine io_write_discard(self, param)
 
          write(LUN, HDRFMT) param%t, nsp, param%lbig_discard
          do i = 1, nsp
-            write(LUN, NAMEFMT) sub, tp_discards%id(i), tp_discards%status(i)
+            write(LUN, NAMEFMT) SUB, tp_discards%id(i), tp_discards%status(i)
             write(LUN, VECFMT) tp_discards%xh(1, i), tp_discards%xh(2, i), tp_discards%xh(3, i)
             write(LUN, VECFMT) tp_discards%vh(1, i), tp_discards%vh(2, i), tp_discards%vh(3, i)
          end do

src/modules/symba_classes.f90

@@ -164,7 +164,6 @@ module symba_classes
       class(symba_pl),      allocatable :: mergesub_list !! List of subtracted bodies in mergers or collisions
       class(symba_pltpenc), allocatable :: pltpenc_list  !! List of massive body-test particle encounters in a single step 
       class(symba_plplenc), allocatable :: plplenc_list  !! List of massive body-massive body encounters in a single step
-      class(symba_pl),      allocatable :: pl_discards   !! Discarded test particle data structure
       integer(I4B)                      :: irec          !! System recursion level
    contains
       procedure :: write_discard    => symba_io_write_discard             !! Write out information about discarded and merged planets and test particles in SyMBA

src/setup/setup.f90

@@ -53,7 +53,6 @@ module subroutine setup_construct_system(system, param)
             allocate(symba_cb :: system%cb)
             allocate(symba_pl :: system%pl)
             allocate(symba_tp :: system%tp)
-            allocate(symba_pl :: system%pl_discards)
             allocate(symba_tp :: system%tp_discards)
             allocate(symba_pl :: system%mergeadd_list)
             allocate(symba_pl :: system%mergesub_list)

src/symba/symba_collision.f90

@@ -67,11 +67,11 @@ module subroutine symba_collision_check_plplenc(self, system, param, t, dt, irec
                   lmask(:) = .false.
                   lmask(ind1(k)) = .not.pl%lcollision(ind1(k))
                   lmask(ind2(k)) = .not.pl%lcollision(ind2(k))
-                  call system%mergesub_list%append(pl, lmask)
-
                   ! Set the collision flag for these to bodies to true in case they become involved in another collision later in the step
                   pl%lcollision([ind1(k), ind2(k)]) = .true.
-                  pl%ldiscard([ind1(k), ind2(k)]) = .true.
+                  pl%lcollision([ind1(k), ind2(k)]) = .true.
+                  pl%status([ind1(k), ind2(k)]) = COLLISION
+                  call system%mergesub_list%append(pl, lmask)
                end do
 
             end if

src/symba/symba_drift.f90

@@ -17,9 +17,9 @@ module subroutine symba_drift_pl(self, system, param, dt)
 
       select type(system)
       class is (symba_nbody_system)
-         self%lmask(:) = self%status(:) == ACTIVE .and. self%levelg(:) == system%irec
+         self%lmask(:) = self%status(:) /= INACTIVE .and. self%levelg(:) == system%irec
          call helio_drift_body(self, system, param, dt)
-         self%lmask(:) = self%status(:) == ACTIVE 
+         self%lmask(:) = self%status(:) /= INACTIVE 
       end select
 
       return
@@ -41,9 +41,9 @@ module subroutine symba_drift_tp(self, system, param, dt)
 
       select type(system)
       class is (symba_nbody_system)
-         self%lmask(:) = self%status(:) == ACTIVE .and. self%levelg(:) == system%irec
+         self%lmask(:) = self%status(:) /= INACTIVE .and. self%levelg(:) == system%irec
          call helio_drift_body(self, system, param, dt)
-         self%lmask(:) = self%status(:) == ACTIVE 
+         self%lmask(:) = self%status(:) /= INACTIVE 
       end select
 
       return

src/symba/symba_io.f90

@@ -225,10 +225,10 @@ module subroutine symba_io_write_discard(self, param)
       character(*), parameter :: PLNAMEFMT = '(I8, 2(1X, E23.16))'
       class(swiftest_body), allocatable :: pltemp
 
-      associate(pl_discards => self%pl_discards, nsppl => self%pl_discards%nbody, pl => self%pl, npl => self%pl%nbody)
+      associate(pl => self%pl, npl => self%pl%nbody, mergesub_list => self%mergesub_list, mergeadd_list => self%mergeadd_list)
          if (self%tp_discards%nbody > 0) call io_write_discard(self, param)
 
-         if (nsppl == 0) return
+         if (mergesub_list%nbody == 0) return
          select case(param%out_stat)
          case('APPEND')
             open(unit = LUN, file = param%discard_out, status = 'OLD', position = 'APPEND', form = 'FORMATTED', iostat = ierr)
@@ -239,14 +239,19 @@ module subroutine symba_io_write_discard(self, param)
             call util_exit(FAILURE)
          end select
          lfirst = .false.
-         if (param%lgr) call pl_discards%pv2v(param) 
-
-      !    write(LUN, HDRFMT) param%t, nsp, param%lbig_discard
-      !    do i = 1, nsp
-      !       write(LUN, NAMEFMT) sub, tp_discards%id(i), tp_discards%status(i)
-      !       write(LUN, VECFMT) tp_discards%xh(1, i), tp_discards%xh(2, i), tp_discards%xh(3, i)
-      !       write(LUN, VECFMT) tp_discards%vh(1, i), tp_discards%vh(2, i), tp_discards%vh(3, i)
-      !    end do
+         if (param%lgr) then
+            call mergesub_list%pv2v(param) 
+            call mergeadd_list%pv2v(param) 
+         end if
+
+         write(LUN, HDRFMT) param%t, mergesub_list%nbody, param%lbig_discard
+         do i = 1, mergesub_list%nbody
+             write(LUN, NAMEFMT) SUB, mergesub_list%id(i), mergesub_list%status(i)
+             write(LUN, VECFMT) mergesub_list%xh(1, i), mergesub_list%xh(2, i), mergesub_list%xh(3, i)
+             write(LUN, VECFMT) mergesub_list%vh(1, i), mergesub_list%vh(2, i), mergesub_list%vh(3, i)
+         end do
+
+         ! This is incomplete until the mergeadd_list methods are completed
       !    if (param%lbig_discard) then
       !          if (param%lgr) then
       !             allocate(pltemp, source = pl)
@@ -265,7 +270,7 @@ module subroutine symba_io_write_discard(self, param)
       !          end do
       !          deallocate(vh)
       !    end if
-      !    close(LUN)
+          close(LUN)
        end associate
 
       return

src/symba/symba_kick.f90

@@ -124,8 +124,8 @@ module subroutine symba_kick_pltpenc(self, system, dt, irec, sgn)
          select type(tp => system%tp)
          class is (symba_tp)
             associate(ind1 => self%index1, ind2 => self%index2)
-               if (pl%nbody > 0) pl%lmask(:) = pl%status(:) == ACTIVE
-               if (tp%nbody > 0) tp%lmask(:) = tp%status(:) == ACTIVE
+               if (pl%nbody > 0) pl%lmask(:) = pl%status(:) /= INACTIVE
+               if (tp%nbody > 0) tp%lmask(:) = tp%status(:) /= INACTIVE
 
                irm1 = irec - 1
                if (sgn < 0) then
@@ -145,7 +145,7 @@ module subroutine symba_kick_pltpenc(self, system, dt, irec, sgn)
                   else
                      lgoodlevel = (pl%levelg(ind1(k)) >= irm1) .and. (tp%levelg(ind2(k)) >= irm1)
                   end if
-                  if ((self%status(k) == ACTIVE) .and. lgoodlevel) then
+                  if ((self%status(k) /= INACTIVE) .and. lgoodlevel) then
                      if (isplpl) then
                         ri = ((pl%rhill(ind1(k))  + pl%rhill(ind2(k)))**2) * (RHSCALE**2) * (RSHELL**(2*irecl))
                         rim1 = ri * (RSHELL**2)

src/util/util_append.f90

@@ -232,7 +232,7 @@ module subroutine util_append_body(self, source, lsource_mask)
       call util_append(self%omega, source%omega, lsource_mask)
       call util_append(self%capm, source%capm, lsource_mask)
 
-      self%nbody = count(self%status(:) == ACTIVE)
+      self%nbody = count(self%status(:) /= INACTIVE)
 
       return
    end subroutine util_append_body

src/util/util_coord.f90

@@ -54,7 +54,7 @@ module subroutine util_coord_h2b_tp(self, cb)
       associate(ntp => self%nbody, xbcb => cb%xb, vbcb => cb%vb, status => self%status, &
                xb => self%xb, xh => self%xh, vb => self%vb, vh => self%vh)
 
-         where(status(1:ntp) == ACTIVE)
+         where(status(1:ntp) /= INACTIVE)
             xb(1, 1:ntp) = xh(1, 1:ntp) + xbcb(1)
             xb(2, 1:ntp) = xh(2, 1:ntp) + xbcb(2)
             xb(3, 1:ntp) = xh(3, 1:ntp) + xbcb(3)
@@ -109,7 +109,7 @@ module subroutine util_coord_b2h_tp(self, cb)
 
       associate(ntp => self%nbody, xbcb => cb%xb, vbcb => cb%vb, xb => self%xb, xh => self%xh, &
                vb => self%vb, vh => self%vh, status => self%status)
-         where(status(1:ntp) == ACTIVE)
+         where(status(1:ntp) /= INACTIVE)
             xh(1, 1:ntp) = xb(1, 1:ntp) - xbcb(1)
             xh(2, 1:ntp) = xb(2, 1:ntp) - xbcb(2)
             xh(3, 1:ntp) = xb(3, 1:ntp) - xbcb(3)

src/util/util_resize.f90

@@ -201,7 +201,7 @@ module subroutine util_resize_body(self, nnew)
       call util_resize(self%capom, nnew)
       call util_resize(self%omega, nnew)
       call util_resize(self%capm, nnew)
-      self%nbody = count(self%status(1:nnew) == ACTIVE)
+      self%nbody = count(self%status(1:nnew) /= INACTIVE)
 
       return
    end subroutine util_resize_body