Merge branch 'debug'

src/fragmentation/fragmentation.f90

@@ -42,7 +42,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin,
       integer(I4B), parameter                 :: NFRAG_MIN = 7 !! The minimum allowable number of fragments (set to 6 because that's how many unknowns are needed in the tangential velocity calculation)
       real(DP)                                :: r_max_start, r_max_start_old, r_max, f_spin 
       real(DP), parameter                     :: Ltol = 10 * epsilon(1.0_DP)
-      real(DP), parameter                     :: Etol = 1e-10_DP
+      real(DP), parameter                     :: Etol = 1e-8_DP
       integer(I4B), parameter                 :: MAXTRY = 3000
       logical, dimension(size(IEEE_ALL))      :: fpe_halting_modes, fpe_quiet_modes
       class(swiftest_nbody_system), allocatable :: tmpsys
@@ -126,17 +126,17 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin,
          if (.not.lfailure) call set_fragment_tan_vel(lfailure)
 
          if (lfailure) then
-            !write(*,*) 'Failed to find tangential velocities'
+            ! write(*,*) 'Failed to find tangential velocities'
          else 
             call set_fragment_radial_velocities(lfailure)
-            !if (lfailure) write(*,*) 'Failed to find radial velocities'
+            ! if (lfailure) write(*,*) 'Failed to find radial velocities'
             if (.not.lfailure) then
                call calculate_system_energy(linclude_fragments=.true.)
                if ((abs(dEtot + Qloss) > Etol) .or. (dEtot > 0.0_DP)) then
-                  !write(*,*) 'Failed due to high energy error: ',dEtot, abs(dEtot + Qloss) / Etol
+                  ! write(*,*) 'Failed due to high energy error: ',dEtot, abs(dEtot + Qloss) / Etol
                   lfailure = .true.
                else if (abs(dLmag) / Lmag_before > Ltol) then
-                  !write(*,*) 'Failed due to high angular momentum error: ', dLmag / Lmag_before
+                  ! write(*,*) 'Failed due to high angular momentum error: ', dLmag / Lmag_before
                   lfailure = .true.
                end if
             end if
@@ -153,11 +153,11 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin,
       ! write(*,        "('  Final diagnostic')")
       ! write(*,        "(' -------------------------------------------------------------------------------------')")
       ! call calculate_system_energy(linclude_fragments=.true.)
-      ! if (lfailure) then
-      !    write(*,*) "symba_frag_pos failed after: ",try," tries"
-      !    do ii = 1, nfrag
-      !       vb_frag(:, ii) = vcom(:)
-      !    end do
+      if (lfailure) then
+         ! write(*,*) "symba_frag_pos failed after: ",try," tries"
+         do ii = 1, nfrag
+            vb_frag(:, ii) = vcom(:)
+         end do
       ! else
       !    write(*,*) "symba_frag_pos succeeded after: ",try," tries"
       !    write(*,        "(' dL_tot should be very small' )")
@@ -166,7 +166,7 @@ module subroutine fragmentation_initialize(system, param, family, x, v, L_spin,
       !    write(*,fmtlabel) ' dE_tot      |', dEtot / abs(Etot_before)
       !    write(*,fmtlabel) ' Qloss       |', -Qloss / abs(Etot_before)
       !    write(*,fmtlabel) ' dE - Qloss  |', (Etot_after - Etot_before + Qloss) / abs(Etot_before)
-      ! end if
+      end if
       ! write(*,        "(' -------------------------------------------------------------------------------------')")
 
       call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes)  ! Save the current halting modes so we can turn them off temporarily

src/io/io.f90

@@ -1660,8 +1660,6 @@ module subroutine io_write_frame_system(self, iu, param)
       call pl%write_frame(iu, param)
       call tp%write_frame(iu, param)
 
-      deallocate(cb, pl, tp)
-
       close(iu, err = 667, iomsg = errmsg)
 
       return

src/orbel/orbel.f90

@@ -16,7 +16,7 @@ module subroutine orbel_el2xv_vec(self, cb)
       if (self%nbody == 0) return
 
       call self%set_mu(cb)
-      do i = 1, self%nbody
+      do concurrent (i = 1:self%nbody)
          call orbel_el2xv(self%mu(i), self%a(i), self%e(i), self%inc(i), self%capom(i), &
                            self%omega(i), self%capm(i), self%xh(:, i), self%vh(:, i))
       end do
@@ -288,7 +288,7 @@ real(DP) pure function orbel_flon(e,icapn)
       end if
 
       return
-   end function  orbel_flon
+   end function orbel_flon
 
 
    !**********************************************************************
@@ -358,7 +358,7 @@ real(DP) pure function orbel_fget(e,capn)
 
       !write(*,*) 'fget : returning without complete convergence'
       return
-   end function  orbel_fget
+   end function orbel_fget
 
 
    !**********************************************************************
@@ -878,18 +878,19 @@ module subroutine orbel_xv2el_vec(self, cb)
       if (self%nbody == 0) return
 
       call self%set_mu(cb)
-      if (.not.allocated(self%a))     allocate(self%a(self%nbody))
-      if (.not.allocated(self%e))     allocate(self%e(self%nbody))
-      if (.not.allocated(self%inc))   allocate(self%inc(self%nbody))
-      if (.not.allocated(self%capom)) allocate(self%capom(self%nbody))
-      if (.not.allocated(self%omega)) allocate(self%omega(self%nbody))
-      if (.not.allocated(self%capm))  allocate(self%capm(self%nbody))
-      do i = 1, self%nbody
+      if (allocated(self%a)) deallocate(self%a); allocate(self%a(self%nbody))
+      if (allocated(self%e)) deallocate(self%e); allocate(self%e(self%nbody))
+      if (allocated(self%inc)) deallocate(self%inc); allocate(self%inc(self%nbody))
+      if (allocated(self%capom)) deallocate(self%capom); allocate(self%capom(self%nbody))
+      if (allocated(self%omega)) deallocate(self%omega); allocate(self%omega(self%nbody))
+      if (allocated(self%capm)) deallocate(self%capm);  allocate(self%capm(self%nbody))
+      do concurrent (i = 1:self%nbody)
          call orbel_xv2el(self%mu(i), self%xh(:, i), self%vh(:, i), self%a(i), self%e(i), self%inc(i),  &
                           self%capom(i), self%omega(i), self%capm(i))
       end do
    end subroutine orbel_xv2el_vec 
 
+
    module pure subroutine orbel_xv2el(mu, x, v, a, e, inc, capom, omega, capm)
       !! author: David A. Minton
       !!