Refactored rscale to dscale for consistency

src/fragmentation/fragmentation.f90

@@ -22,7 +22,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin,
       real(DP),                                                  intent(inout) :: Qloss !! Energy lost during the collision
       logical,                                                   intent(out)   :: lfailure !! Answers the question: Should this have been a merger instead?
       ! Internals
-      real(DP)                                :: mscale, rscale, vscale, tscale, Lscale, Escale ! Scale factors that reduce quantities to O(~1) in the collisional system
+      real(DP)                                :: mscale, dscale, vscale, tscale, Lscale, Escale ! Scale factors that reduce quantities to O(~1) in the collisional system
       real(DP)                                :: mtot 
       real(DP), dimension(NDIM)               :: xcom, vcom
       integer(I4B)                            :: ii
@@ -46,6 +46,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin,
       integer(I4B), parameter                 :: MAXTRY = 3000
       integer(I4B), parameter                 :: TANTRY = 3
       logical, dimension(size(IEEE_ALL))      :: fpe_halting_modes, fpe_quiet_modes
+      class(swiftest_parameters), allocatable :: tmpparam
 
       if (nfrag < NFRAG_MIN) then
          write(*,*) "symba_frag_pos needs at least ",NFRAG_MIN," fragments, but only ",nfrag," were given."
@@ -59,7 +60,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin,
 
       f_spin = 0.05_DP
       mscale = 1.0_DP
-      rscale = 1.0_DP
+      dscale = 1.0_DP
       vscale = 1.0_DP
       tscale = 1.0_DP
       Lscale = 1.0_DP
@@ -166,27 +167,27 @@ subroutine set_scale_factors()
 
          ! Set scale factors
          Escale = 0.5_DP * (mass(1) * dot_product(v(:,1), v(:,1)) + mass(2) * dot_product(v(:,2), v(:,2)))
-         rscale = sum(radius(:))
-         mscale = sqrt(Escale * rscale / param%GU) 
+         dscale = sum(radius(:))
+         mscale = sqrt(Escale * dscale / param%GU) 
          vscale = sqrt(Escale / mscale) 
-         tscale = rscale / vscale 
-         Lscale = mscale * rscale * vscale
+         tscale = dscale / vscale 
+         Lscale = mscale * dscale * vscale
 
-         xcom(:) = xcom(:) / rscale
+         xcom(:) = xcom(:) / dscale
          vcom(:) = vcom(:) / vscale
 
          mtot = mtot / mscale
          mass = mass / mscale
-         radius = radius / rscale
-         x = x / rscale
+         radius = radius / dscale
+         x = x / dscale
          v = v / vscale
          L_spin = L_spin / Lscale
          do i = 1, 2
             rot(:,i) = L_spin(:,i) / (mass(i) * radius(i)**2 * Ip(3, i))
          end do
 
          m_frag = m_frag / mscale
-         rad_frag = rad_frag / rscale
+         rad_frag = rad_frag / dscale
          Qloss = Qloss / Escale
 
          return
@@ -201,23 +202,23 @@ subroutine restore_scale_factors()
 
          call ieee_set_halting_mode(IEEE_ALL,.false.)
          ! Restore scale factors
-         xcom(:) = xcom(:) * rscale
+         xcom(:) = xcom(:) * dscale
          vcom(:) = vcom(:) * vscale
 
          mtot = mtot * mscale
          mass = mass * mscale
-         radius = radius * rscale
-         x = x * rscale
+         radius = radius * dscale
+         x = x * dscale
          v = v * vscale
          L_spin = L_spin * Lscale
          do i = 1, 2
             rot(:,i) = L_spin(:,i) * (mass(i) * radius(i)**2 * Ip(3, i))
          end do
 
          m_frag = m_frag * mscale
-         rad_frag = rad_frag * rscale
+         rad_frag = rad_frag * dscale
          rot_frag = rot_frag / tscale
-         x_frag = x_frag * rscale
+         x_frag = x_frag * dscale
          v_frag = v_frag * vscale
          Qloss = Qloss * Escale
 
@@ -240,7 +241,7 @@ subroutine restore_scale_factors()
          Lmag_after = Lmag_after * Lscale 
 
          mscale = 1.0_DP
-         rscale = 1.0_DP
+         dscale = 1.0_DP
          vscale = 1.0_DP
          tscale = 1.0_DP
          Lscale = 1.0_DP
@@ -314,7 +315,6 @@ subroutine calculate_system_energy(linclude_fragments)
          logical, dimension(:), allocatable :: ltmp
          logical :: lk_plpl
          class(swiftest_nbody_system), allocatable :: tmpsys
-         class(swiftest_parameters), allocatable :: tmpparam
 
          ! Because we're making a copy of symba_pl with the excludes/fragments appended, we need to deallocate the
          ! big k_plpl array and recreate it when we're done, otherwise we run the risk of blowing up the memory by
@@ -344,12 +344,12 @@ subroutine calculate_system_energy(linclude_fragments)
             ltmp(1:npl) = .true.
             call tmpsys%pl%setup(npl_new, param)
             call tmpsys%pl%fill(pl, ltmp)
-            call tmpsys%rescale(tmpparam, mscale, rscale, tscale)
+            call tmpsys%rescale(tmpparam, mscale, dscale, tscale)
 
             if (linclude_fragments) then ! Append the fragments if they are included
                ! Energy calculation requires the fragments to be in the system barcyentric frame, s
                tmpsys%pl%mass(npl+1:npl_new) = m_frag(1:nfrag)
-               tmpsys%pl%Gmass(npl+1:npl_new) = m_frag(1:nfrag) * param%GU * rscale**3 * mscale / tscale**2 
+               tmpsys%pl%Gmass(npl+1:npl_new) = m_frag(1:nfrag) * param%GU * dscale**3 * mscale / tscale**2 
                tmpsys%pl%radius(npl+1:npl_new) = rad_frag(1:nfrag)
                tmpsys%pl%xb(:,npl+1:npl_new) =  xb_frag(:,1:nfrag)
                tmpsys%pl%vb(:,npl+1:npl_new) =  vb_frag(:,1:nfrag)