Cleaned up unused variables, type conversion inconsistencies (except …

…for unit conversions using quad precision) and other things flagged by gfortran warnings

src/base/base_module.f90

@@ -2078,7 +2078,7 @@ subroutine base_util_unique_DP(input_array, output_array, index_map)
             n = n + 1
             lo = minval(input_array(:), mask=input_array(:) > lo)
             unique_array(n) = lo
-            where(input_array(:) == lo) index_map(:) = n
+            where(abs(input_array(:) - lo) < epsilon(1.0_DP) * lo) index_map(:) = n
             if (lo >= hi) exit
          enddo
          allocate(output_array(n), source=unique_array(1:n)) 

src/collision/collision_check.f90

@@ -77,7 +77,8 @@ module subroutine collision_check_plpl(self, nbody_system, param, t, dt, irec, l
       ! Internals
       logical, dimension(:), allocatable        :: lcollision, lmask
       real(DP), dimension(NDIM)                 :: xr, vr
-      integer(I4B)                              :: i, j, k, nenc
+      integer(I4B)                              :: i, j
+      integer(I8B)                              :: k, nenc
       real(DP)                                  :: rlim, Gmtot
       logical                                   :: lany_closest
 
@@ -102,9 +103,9 @@ module subroutine collision_check_plpl(self, nbody_system, param, t, dt, irec, l
             self%lclosest(:) = .false.
 
 #ifdef DOCONLOC
-            do concurrent(k = 1:nenc, lmask(k)) shared(self,pl,lmask, dt, lcollision) local(i,j,xr,vr,rlim,Gmtot)
+            do concurrent(k = 1_I8B:nenc, lmask(k)) shared(self,pl,lmask, dt, lcollision) local(i,j,xr,vr,rlim,Gmtot)
 #else
-            do concurrent(k = 1:nenc, lmask(k))
+            do concurrent(k = 1_I8B:nenc, lmask(k))
 #endif
                i = self%index1(k)
                j = self%index2(k)
@@ -120,7 +121,7 @@ module subroutine collision_check_plpl(self, nbody_system, param, t, dt, irec, l
 
             if (lany_collision .or. lany_closest) then
                call pl%rh2rb(nbody_system%cb) ! Update the central body barycenteric position vector to get us out of DH and into bary
-               do k = 1, nenc
+               do k = 1_I8B, nenc
                   if (.not.lcollision(k) .and. .not. self%lclosest(k)) cycle
                   i = self%index1(k)
                   j = self%index2(k)
@@ -178,7 +179,8 @@ module subroutine collision_check_pltp(self, nbody_system, param, t, dt, irec, l
       ! Internals
       logical, dimension(:), allocatable        :: lcollision, lmask
       real(DP), dimension(NDIM)                 :: xr, vr
-      integer(I4B)                              :: i, j, k, nenc
+      integer(I4B)                              :: i, j
+      integer(I8B)                              :: k, nenc
       logical                                   :: lany_closest
       character(len=STRMAX)                     :: timestr, idstri, idstrj, message
       class(collision_list_pltp), allocatable       :: tmp       
@@ -194,13 +196,13 @@ module subroutine collision_check_pltp(self, nbody_system, param, t, dt, irec, l
 
             nenc = self%nenc
             allocate(lmask(nenc))
-             lmask(:) = (self%status(1:nenc) == ACTIVE) 
+            lmask(:) = (self%status(1:nenc) == ACTIVE) 
             select type(pl)
             class is (symba_pl)
-            select type(tp)
-            class is (symba_tp)
-               lmask(:) = lmask(:) .and. (tp%levelg(self%index2(1:nenc)) >= irec)
-            end select
+               select type(tp)
+               class is (symba_tp)
+                  lmask(:) = lmask(:) .and. (tp%levelg(self%index2(1:nenc)) >= irec)
+               end select
             end select
             if (.not.any(lmask(:))) return
 
@@ -209,9 +211,9 @@ module subroutine collision_check_pltp(self, nbody_system, param, t, dt, irec, l
             self%lclosest(:) = .false.
 
 #ifdef DOCONLOC
-            do concurrent(k = 1:nenc, lmask(k)) shared(self,pl,tp,lmask, dt, lcollision) local(i,j,xr,vr)
+            do concurrent(k = 1_I8B:nenc, lmask(k)) shared(self,pl,tp,lmask, dt, lcollision) local(i,j,xr,vr)
 #else
-            do concurrent(k = 1:nenc, lmask(k))
+            do concurrent(k = 1_I8B:nenc, lmask(k))
 #endif
                i = self%index1(k)
                j = self%index2(k)

src/collision/collision_generate.f90

@@ -170,7 +170,8 @@ module subroutine collision_generate_merge(self, nbody_system, param, t)
       class(base_parameters),   intent(inout) :: param        !! Current run configuration parameters 
       real(DP),                 intent(in)    :: t            !! The time of the collision
       ! Internals
-      integer(I4B)                              :: i, j, k, ibiggest
+      integer(I4B)                              :: i, j, ibiggest
+      integer(I8B)                              :: k
       real(DP), dimension(NDIM)                 :: L_spin_new, L_residual
       real(DP)                                  :: volume
       character(len=STRMAX) :: message
@@ -231,7 +232,7 @@ module subroutine collision_generate_merge(self, nbody_system, param, t)
                   call collision_util_velocity_torque(-L_residual(:), fragments%mass(1), fragments%rb(:,1), fragments%vb(:,1))
 
                   ! Update any encounter lists that have the removed bodies in them so that they instead point to the new body
-                  do k = 1, nbody_system%plpl_encounter%nenc
+                  do k = 1_I8B, nbody_system%plpl_encounter%nenc
                      do j = 1, impactors%ncoll
                         i = impactors%id(j)
                         if (i == ibiggest) cycle

src/collision/collision_io.f90

@@ -365,7 +365,7 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param)
       class(encounter_storage), intent(inout) :: history !! Collision history object
       class(base_parameters),      intent(inout) :: param   !! Current run configuration parameters
       ! Internals
-      integer(I4B)           :: i, idslot, old_mode, npl, stage
+      integer(I4B)           :: i, idslot, old_mode, npl, stage, tmp
       character(len=NAMELEN) :: charstring
       class(swiftest_pl), allocatable :: pl
 
@@ -427,7 +427,7 @@ module subroutine collision_io_netcdf_write_frame_snapshot(self, history, param)
                call netcdf_io_check( nf90_put_var(nc%id, nc%L_spin_varid,   collider%L_spin(:,:),  start=[1, 1, eslot], count=[NDIM, 2, 1]), "collision_io_netcdf_write_frame_snapshot nf90_put_var L_spin_varid before" )
             end if
 
-            call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode) )
+            call netcdf_io_check( nf90_set_fill(nc%id, old_mode, tmp) )
          end associate
       end select
       return

src/collision/collision_module.f90

@@ -126,7 +126,7 @@ module collision
       real(DP),                  dimension(:),   allocatable :: rmag         !! Array of radial distance magnitudes of individual fragments in the collisional coordinate frame 
       real(DP),                  dimension(:),   allocatable :: vmag         !! Array of radial distance magnitudes of individual fragments in the collisional coordinate frame 
       real(DP),                  dimension(:),   allocatable :: rotmag       !! Array of rotation magnitudes of individual fragments 
-      integer(I1B),              dimension(:),   allocatable :: origin_body  !! Array of indices indicating which impactor body (1 or 2) the fragment originates from
+      integer(I4B),              dimension(:),   allocatable :: origin_body  !! Array of indices indicating which impactor body (1 or 2) the fragment originates from
       real(DP),                  dimension(NDIM)             :: L_orbit_tot  !! Orbital angular momentum vector of all fragments
       real(DP),                  dimension(NDIM)             :: L_spin_tot   !! Spin angular momentum vector of all fragments
       real(DP),                  dimension(:,:), allocatable :: L_orbit      !! Orbital angular momentum vector of each individual fragment

src/collision/collision_resolve.f90

@@ -191,8 +191,9 @@ module subroutine collision_resolve_extract_plpl(self, nbody_system, param)
       logical,      dimension(:), allocatable :: lplpl_collision
       logical,      dimension(:), allocatable :: lplpl_unique_parent
       integer(I4B), dimension(:), pointer     :: plparent
-      integer(I4B), dimension(:), allocatable :: collision_idx, unique_parent_idx
-      integer(I4B)                            :: i, index_coll, ncollisions, nunique_parent, nplplenc
+      integer(I4B)                            :: nunique_parent
+      integer(I8B)                            :: ncollisions, index_coll, k, nplplenc
+      integer(I8B), dimension(:), allocatable :: unique_parent_idx, collision_idx
 
       select type(nbody_system)
       class is (swiftest_nbody_system)
@@ -201,14 +202,14 @@ module subroutine collision_resolve_extract_plpl(self, nbody_system, param)
          associate(idx1 => self%index1, idx2 => self%index2, plparent => pl%kin%parent)
             nplplenc = self%nenc
             allocate(lplpl_collision(nplplenc))
-            lplpl_collision(:) = self%status(1:nplplenc) == COLLIDED
+            lplpl_collision(:) = self%status(1_I8B:nplplenc) == COLLIDED
             if (.not.any(lplpl_collision)) return 
             ! Collisions have been detected in this step. So we need to determine which of them are between unique bodies.
 
             ! Get the subset of pl-pl encounters that lead to a collision
             ncollisions = count(lplpl_collision(:))
             allocate(collision_idx(ncollisions))
-            collision_idx = pack([(i, i=1, nplplenc)], lplpl_collision)
+            collision_idx = pack([(k, k=1_I8B, nplplenc)], lplpl_collision)
 
             ! Get the subset of collisions that involve a unique pair of parents
             allocate(lplpl_unique_parent(ncollisions))
@@ -223,7 +224,7 @@ module subroutine collision_resolve_extract_plpl(self, nbody_system, param)
             ! due to restructuring of parent/child relationships when there are large numbers of multi-body collisions in a single
             ! step
             lplpl_unique_parent(:) = .true.
-            do index_coll = 1, ncollisions
+            do index_coll = 1_I8B, ncollisions
                associate(ip1 => plparent(idx1(collision_idx(index_coll))), ip2 => plparent(idx2(collision_idx(index_coll))))
                   lplpl_unique_parent(:) = .not. ( any(plparent(idx1(unique_parent_idx(:))) == ip1) .or. &
                                                    any(plparent(idx2(unique_parent_idx(:))) == ip1) .or. &
@@ -537,7 +538,8 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec)
       character(len=STRMAX) :: timestr, idstr
       integer(I4B), dimension(2) :: idx_parent       !! Index of the two bodies considered the "parents" of the collision
       logical  :: lgoodcollision
-      integer(I4B) :: i, j, nnew, loop, ncollisions
+      integer(I4B) :: i, j, nnew, loop
+      integer(I8B) :: k, ncollisions
       integer(I4B), dimension(:), allocatable :: idnew
       integer(I4B), parameter :: MAXCASCADE = 1000
 
@@ -576,9 +578,9 @@ module subroutine collision_resolve_plpl(self, nbody_system, param, t, dt, irec)
                   call swiftest_io_log_one_message(COLLISION_LOG_OUT, "***********************************************************" // &
                                                             "***********************************************************")
 
-                  do i = 1, ncollisions
-                     idx_parent(1) = pl%kin(idx1(i))%parent
-                     idx_parent(2) = pl%kin(idx2(i))%parent
+                  do k = 1_I8B, ncollisions
+                     idx_parent(1) = pl%kin(idx1(k))%parent
+                     idx_parent(2) = pl%kin(idx2(k))%parent
                      call impactors%consolidate(nbody_system, param, idx_parent, lgoodcollision)
                      if ((.not. lgoodcollision) .or. any(pl%status(idx_parent(:)) /= COLLIDED)) cycle
 

src/encounter/encounter_check.f90

@@ -636,7 +636,7 @@ module subroutine encounter_check_collapse_ragged_list(ragged_list, n1, nenc, in
       ! Internals
       integer(I4B) :: i
       integer(I8B) :: j1, j0, nenci
-      integer(I4B), dimension(n1) :: ibeg
+      integer(I8B), dimension(n1) :: ibeg
 
       associate(nenc_arr => ragged_list(:)%nenc)
          nenc = sum(nenc_arr(:))
@@ -865,7 +865,7 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r
          if (loverlap(i)) then
             ibeg =  self%aabb%ibeg(i) + 1_I8B
             iend =  self%aabb%iend(i) - 1_I8B
-            nbox = iend - ibeg + 1
+            nbox = int(iend - ibeg, kind=I4B) + 1
             call encounter_check_all_sweep_one(i, nbox, r1(1,i), r1(2,i), r1(3,i), v1(1,i), v1(2,i), v1(3,i), &
                                                          xind(ibeg:iend), yind(ibeg:iend), zind(ibeg:iend),&
                                                          vxind(ibeg:iend), vyind(ibeg:iend), vzind(ibeg:iend), &
@@ -881,7 +881,7 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, r1, v1, r
          if (loverlap(i)) then
             ibeg =  self%aabb%ibeg(i) + 1_I8B
             iend =  self%aabb%iend(i) - 1_I8B
-            nbox = iend - ibeg + 1
+            nbox = int(iend - ibeg, kind=I4B) + 1
             ii = i - n1
             call encounter_check_all_sweep_one(ii, nbox, r2(1,ii), r2(2,ii), r2(3,ii), v2(1,ii), v2(2,ii), v2(3,ii), &
                                                           xind(ibeg:iend), yind(ibeg:iend), zind(ibeg:iend),&
@@ -958,7 +958,7 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, r, v, renc, dt
          if (loverlap(i)) then
             ibeg =  self%aabb%ibeg(i) + 1_I8B
             iend =  self%aabb%iend(i) - 1_I8B
-            nbox = int(iend - ibeg + 1, kind=I4B)
+            nbox = int(iend - ibeg, kind=I4B) + 1
             lencounteri(ibeg:iend) = .true.
             call encounter_check_all_sweep_one(i, nbox, r(1,i), r(2,i), r(3,i), v(1,i), v(2,i), v(3,i), &
                                                       xind(ibeg:iend), yind(ibeg:iend), zind(ibeg:iend),&

src/encounter/encounter_io.f90

@@ -231,7 +231,7 @@ module subroutine encounter_io_netcdf_write_frame_snapshot(self, history, param)
       class(base_parameters),      intent(inout) :: param             !! Current run configuration parameters
 
       ! Internals
-      integer(I4B)           :: i, idslot, old_mode, npl, ntp
+      integer(I4B)           :: i, idslot, old_mode, npl, ntp, tmp
       character(len=STRMAX) :: charstring
 
       select type(param)
@@ -284,7 +284,7 @@ module subroutine encounter_io_netcdf_write_frame_snapshot(self, history, param)
                call netcdf_io_check( nf90_put_var(nc%id, nc%ptype_varid, charstring, start=[1, idslot], count=[NAMELEN, 1]), "encounter_io_netcdf_write_frame_snapshot nf90_put_var tp particle_type_varid"  )
             end do
 
-            call netcdf_io_check( nf90_set_fill(nc%id, old_mode, old_mode) )
+            call netcdf_io_check( nf90_set_fill(nc%id, old_mode, tmp) )
          end associate
       end select
       end select

src/fraggle/fraggle_generate.f90

@@ -361,8 +361,6 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu
       integer(I4B) :: i, j, loop, istart, nfrag, npl, ntp
       logical :: lsupercat, lhitandrun
       integer(I4B), parameter :: MAXLOOP = 10000
-      real(DP), parameter :: cloud_size_scale_factor = 3.0_DP ! Scale factor to apply to the size of the cloud relative to the distance from the impact point. 
-                                                              ! A larger value puts more space between fragments initially
       real(DP), parameter :: rbuffer = 1.01_DP ! Body radii are inflated by this scale factor to prevent secondary collisions 
       real(DP), parameter :: pack_density = 0.5236_DP ! packing density of loose spheres
 

src/fraggle/fraggle_util.f90

@@ -108,7 +108,6 @@ module subroutine fraggle_util_set_mass_dist(self, param)
       integer(I4B), parameter :: iMrem = 3
       integer(I4B), parameter :: NFRAGMIN = iMrem + 2 !! Minimum number of fragments that can be generated 
       integer(I4B), dimension(:), allocatable :: ind
-      integer(I4B), parameter :: MAXLOOP = 20
       logical :: flipper
       logical, dimension(size(IEEE_ALL))      :: fpe_halting_modes
 
@@ -174,7 +173,7 @@ module subroutine fraggle_util_set_mass_dist(self, param)
             mass(2) = impactors%mass_dist(iMslr)
 
             ! Recompute the slope parameter beta so that we span the complete size range
-            if (Mslr == min_mfrag) Mslr = Mslr + impactors%mass_dist(iMrem) / nfrag
+            if (abs(Mslr - min_mfrag) < epsilon(min_mfrag) * min_mfrag) Mslr = Mslr + impactors%mass_dist(iMrem) / nfrag
             mremaining = impactors%mass_dist(iMrem)
 
             ! The mass will be distributed in a power law where N>M=(M/Mslr)**(-beta/3)
@@ -225,8 +224,8 @@ module subroutine fraggle_util_set_mass_dist(self, param)
 
             ! For catastrophic impacts, we will assign each of the n>2 fragments to one of the two original bodies so that the fragment cloud occupies 
             ! roughly the same space as both original bodies. For all other disruption cases, we use body 2 as the center of the cloud.
-            fragments%origin_body(1) = 1
-            fragments%origin_body(2) = 2
+            fragments%origin_body(1) = 1_I1B
+            fragments%origin_body(2) = 2_I1B
             if (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) then
                mcumul = fragments%mass(1)
                flipper = .true.

src/rmvs/rmvs_step.f90

@@ -541,9 +541,8 @@ subroutine rmvs_peri_tp(tp, pl, t, dt, lfirst, inner_index, ipleP, param)
       integer(I4B),               intent(in)    :: ipleP     !!  index of RMVS planet being closely encountered
       class(swiftest_parameters), intent(in)    :: param    !! Current run configuration parameters
       ! Internals
-      integer(I4B)              :: i, id1, id2
-      real(DP)                  :: r2, mu, rhill2, vdotr, a, peri, capm, tperi, rpl
-      real(DP), dimension(NDIM) :: rh1, rh2, vh1, vh2
+      integer(I4B)              :: i
+      real(DP)                  :: r2, mu, rhill2, vdotr, a, peri, capm, tperi
 
       rhill2 = pl%rhill(ipleP)**2
       mu = pl%Gmass(ipleP)

src/swiftest/swiftest_module.f90

@@ -671,11 +671,13 @@ module subroutine swiftest_io_netcdf_get_t0_values_system(self, nc, param)
          class(swiftest_parameters),        intent(inout) :: param !! Current run configuration parameters 
       end subroutine swiftest_io_netcdf_get_t0_values_system
 
-      module subroutine swiftest_io_netcdf_get_valid_masks(self, plmask, tpmask)
+      module subroutine swiftest_io_netcdf_get_valid_masks(self, plmask, tpmask, plmmask, Gmtiny)
          implicit none
-         class(swiftest_netcdf_parameters),  intent(inout) :: self   !! Parameters used to identify a particular NetCDF dataset
-         logical, dimension(:), allocatable, intent(out)   :: plmask !! Logical mask indicating which bodies are massive bodies
-         logical, dimension(:), allocatable, intent(out)   :: tpmask !! Logical mask indicating which bodies are test particles
+         class(swiftest_netcdf_parameters),  intent(inout)          :: self    !! Parameters used to identify a particular NetCDF dataset
+         logical, dimension(:), allocatable, intent(out)            :: plmask  !! Logical mask indicating which bodies are massive bodies
+         logical, dimension(:), allocatable, intent(out)            :: tpmask  !! Logical mask indicating which bodies are test particles
+         logical, dimension(:), allocatable, intent(out), optional  :: plmmask !! Logical mask indicating which bodies are fully interacting massive bodies
+         real(DP),                           intent(in),  optional  :: Gmtiny  !! The cutoff G*mass between semi-interacting and fully interacting massive bodies
       end subroutine swiftest_io_netcdf_get_valid_masks
 
       module subroutine swiftest_io_netcdf_initialize_output(self, param)