diff --git a/src/fragmentation/fragmentation.f90 b/src/fragmentation/fragmentation.f90 index ad28edf6c..e1b69dd44 100644 --- a/src/fragmentation/fragmentation.f90 +++ b/src/fragmentation/fragmentation.f90 @@ -525,21 +525,21 @@ subroutine set_fragment_tan_vel(lerr) call calculate_system_energy(linclude_fragments=.true.) ke_frag_budget = -dEtot - Qloss - !write(*,*) '***************************************************' - !write(*,*) 'Original dis : ',norm2(x(:,2) - x(:,1)) - !write(*,*) 'r_max : ',r_max - !write(*,*) 'f_spin : ',f_spin - !write(*,*) '***************************************************' - !write(*,*) 'Energy balance so far: ' - !write(*,*) 'ke_frag_budget : ',ke_frag_budget - !write(*,*) 'ke_orbit_before: ',ke_orbit_before - !write(*,*) 'ke_orbit_after : ',ke_orbit_after - !write(*,*) 'ke_spin_before : ',ke_spin_before - !write(*,*) 'ke_spin_after : ',ke_spin_after - !write(*,*) 'pe_before : ',pe_before - !write(*,*) 'pe_after : ',pe_after - !write(*,*) 'Qloss : ',Qloss - !write(*,*) '***************************************************' + write(*,*) '***************************************************' + write(*,*) 'Original dis : ',norm2(x(:,2) - x(:,1)) + write(*,*) 'r_max : ',r_max + write(*,*) 'f_spin : ',f_spin + write(*,*) '***************************************************' + write(*,*) 'Energy balance so far: ' + write(*,*) 'ke_frag_budget : ',ke_frag_budget + write(*,*) 'ke_orbit_before: ',ke_orbit_before + write(*,*) 'ke_orbit_after : ',ke_orbit_after + write(*,*) 'ke_spin_before : ',ke_spin_before + write(*,*) 'ke_spin_after : ',ke_spin_after + write(*,*) 'pe_before : ',pe_before + write(*,*) 'pe_after : ',pe_after + write(*,*) 'Qloss : ',Qloss + write(*,*) '***************************************************' if (ke_frag_budget < 0.0_DP) then write(*,*) 'Negative ke_frag_budget: ',ke_frag_budget r_max_start = r_max_start / 2 @@ -604,11 +604,11 @@ subroutine set_fragment_tan_vel(lerr) ! If we are over the energy budget, flag this as a failure so we can try again lerr = (ke_radial < 0.0_DP) - !write(*,*) 'Tangential' - !write(*,*) 'ke_frag_budget: ',ke_frag_budget - !write(*,*) 'ke_frag_orbit : ',ke_frag_orbit - !write(*,*) 'ke_frag_spin : ',ke_frag_spin - !write(*,*) 'ke_radial : ',ke_radial + write(*,*) 'Tangential' + write(*,*) 'ke_frag_budget: ',ke_frag_budget + write(*,*) 'ke_frag_orbit : ',ke_frag_orbit + write(*,*) 'ke_frag_spin : ',ke_frag_spin + write(*,*) 'ke_radial : ',ke_radial return diff --git a/src/util/util_get_energy_momentum.f90 b/src/util/util_get_energy_momentum.f90 index fd331bcc3..4ec9c6c36 100644 --- a/src/util/util_get_energy_momentum.f90 +++ b/src/util/util_get_energy_momentum.f90 @@ -76,7 +76,7 @@ module subroutine util_get_energy_momentum_system(self, param) !$omp simd do i = 1, npl associate(px => pl%xh(1,i), py => pl%xh(2,i), pz => pl%xh(3,i)) - pecb(i) = -cb%mass * pl%mass(i) / sqrt(px**2 + py**2 + pz**2) + pecb(i) = -param%GU * cb%mass * pl%mass(i) / sqrt(px**2 + py**2 + pz**2) end associate end do @@ -100,7 +100,7 @@ module subroutine util_get_energy_momentum_system(self, param) irh(i) = 1.0_DP / norm2(pl%xh(:,i)) end do call obl_pot(npl, cb%mass, pl%mass, cb%j2rp2, cb%j4rp4, pl%xh, irh, oblpot) - system%pe = system%pe + oblpot + system%pe = system%pe + param%GU * oblpot end if system%Lorbit(1) = sum(Lplorbitx(1:npl), lstatus(1:npl))