diff --git a/src/io/io_conservation_report.f90 b/src/io/io_conservation_report.f90 index 94020a34b..f08f6e64d 100644 --- a/src/io/io_conservation_report.f90 +++ b/src/io/io_conservation_report.f90 @@ -15,10 +15,12 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param, type(user_input_parameters), intent(in) :: param !! Input colleciton of user-defined parameters logical, intent(in) :: lterminal !! Indicates whether to output information to the terminal screen ! Internals - real(DP), dimension(NDIM), save :: Ltot_orig, Ltot_last - real(DP), save :: Eorbit_orig, Mtot_orig, Lmag_orig, ke_orb_last, ke_spin_last, pe_last, Eorbit_last - real(DP) :: ke_orbit, ke_spin, pe, Eorbit - real(DP), dimension(NDIM) :: Ltot_now + real(DP), dimension(NDIM) :: Ltot_now, Lorbit_now, Lspin_now + real(DP), dimension(NDIM), save :: Ltot_orig, Lorbit_orig, Lspin_orig + real(DP), dimension(NDIM), save :: Ltot_last, Lorbit_last, Lspin_last + real(DP), save :: Eorbit_orig, Mtot_orig, Lmag_orig + real(DP), save :: ke_orbit_last, ke_spin_last, pe_last, Eorbit_last + real(DP) :: ke_orbit_now, ke_spin_now, pe_now, Eorbit_now real(DP) :: Eorbit_error, Etotal_error, Ecoll_error real(DP) :: Mtot_now, Merror real(DP) :: Lmag_now, Lerror @@ -42,39 +44,53 @@ module subroutine io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param, write(egyiu,egyheader) end if end if - call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe, Eorbit, Ltot_now) + call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_now, ke_spin_now, pe_now, Lorbit_now, Lspin_now) + Eorbit_now = ke_orbit_now + ke_spin_now + pe_now + Ltot_now(:) = Lorbit_now(:) + Lspin_now(:) + Lescape(:) Mtot_now = dMcb + sum(mass(2:npl)) + Mcb_initial + Mescape Ltot_now(:) = Lescape(:) + Ltot_now(:) if (lfirst) then - Eorbit_orig = Eorbit + Eorbit_orig = Eorbit_now Mtot_orig = Mtot_now + Lorbit_orig(:) = Lorbit_now(:) + Lspin_orig(:) = Lspin_now(:) Ltot_orig(:) = Ltot_now(:) Lmag_orig = norm2(Ltot_orig(:)) lfirst = .false. end if - write(egyiu,egyfmt) t, Eorbit, Ecollisions, Ltot_now, Mtot_now + write(egyiu,egyfmt) t, Eorbit_now, Ecollisions, Ltot_now, Mtot_now flush(egyiu) if (.not.lfirst .and. lterminal) then Lmag_now = norm2(Ltot_now) Lerror = norm2(Ltot_now - Ltot_orig) / Lmag_orig - Eorbit_error = (Eorbit - Eorbit_orig) / abs(Eorbit_orig) + Eorbit_error = (Eorbit_now - Eorbit_orig) / abs(Eorbit_orig) Ecoll_error = Ecollisions / abs(Eorbit_orig) - Etotal_error = (Eorbit - Ecollisions - Eorbit_orig - Euntracked) / abs(Eorbit_orig) + Etotal_error = (Eorbit_now - Ecollisions - Eorbit_orig - Euntracked) / abs(Eorbit_orig) Merror = (Mtot_now - Mtot_orig) / Mtot_orig write(*, egytermfmt) Lerror, Ecoll_error, Etotal_error, Merror if (Ecoll_error > 0.0_DP) then write(*,*) 'Something has gone wrong! Collisional energy should not be positive!' - write(*,*) 'dke_orbit: ',(ke_orbit - ke_orb_last) / abs(Eorbit_orig) - write(*,*) 'dke_spin : ',(ke_spin - ke_spin_last) / abs(Eorbit_orig) - write(*,*) 'dpe : ',(pe - pe_last) / abs(Eorbit_orig) + write(*,*) 'dke_orbit: ',(ke_orbit_now - ke_orbit_last) / abs(Eorbit_orig) + write(*,*) 'dke_spin : ',(ke_spin_now - ke_spin_last) / abs(Eorbit_orig) + write(*,*) 'dpe : ',(pe_now - pe_last) / abs(Eorbit_orig) write(*,*) end if end if - ke_orb_last = ke_orbit - ke_spin_last = ke_spin - pe_last = pe - Eorbit_last = Eorbit + write(*,*) 'Angular momentum changes since last time' + write(*,*) ' dLorbit : ',(Lorbit_now(:) - Lorbit_last(:)) / Lmag_orig + write(*,*) '|dLorbit|: ',norm2(Lorbit_now(:) - Lorbit_last(:)) / Lmag_orig + write(*,*) ' dLspin : ',(Lspin_now(:) - Lspin_last(:)) / Lmag_orig + write(*,*) '|dLspin| : ',norm2(Lspin_now(:) - Lspin_last(:)) / Lmag_orig + write(*,*) ' dLtot : ',(Ltot_now(:) - Ltot_last(:)) / Lmag_orig + write(*,*) '|dLtot| : ',norm2(Ltot_now(:) - Ltot_last(:)) / Lmag_orig + ke_orbit_last = ke_orbit_now + ke_spin_last = ke_spin_now + pe_last = pe_now + Eorbit_last = Eorbit_now + Lorbit_last(:) = Lorbit_now(:) + Lspin_last(:) = Lspin_now(:) + Ltot_last(:) = Ltot_now(:) end associate return diff --git a/src/main/swiftest_symba.f90 b/src/main/swiftest_symba.f90 index b39890096..580346b1c 100644 --- a/src/main/swiftest_symba.f90 +++ b/src/main/swiftest_symba.f90 @@ -21,9 +21,9 @@ program swiftest_symba logical :: lfrag_add, ldiscard_pl, ldiscard_tp integer(I4B) :: nplm, ntp, ntp0, nsppl, nsptp, iout, idump, iloop, i integer(I4B) :: nplplenc, npltpenc, nmergeadd, nmergesub - real(DP) :: t, tfrac, tbase, msys, Eorbit_orig + real(DP) :: t, tfrac, tbase, msys real(DP) :: Ecollision, Eorbit_before, Eorbit_after, ke_orbit_before, ke_spin_before, ke_orbit_after, ke_spin_after, pe_before, pe_after - real(DP), dimension(NDIM) :: Ltot + real(DP), dimension(NDIM) :: Lorbit, Lspin character(STRMAX) :: inparfile type(symba_pl) :: symba_plA type(symba_tp) :: symba_tpA @@ -153,7 +153,8 @@ program swiftest_symba if (param%lenergy) then call coord_h2b(npl, symba_plA%helio%swiftest, msys) call io_conservation_report(t, symba_plA, npl, j2rp2, j4rp4, param, lterminal=.false.) - call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_before, ke_spin_before, pe_before, Eorbit_orig, Ltot) + call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_before, ke_spin_before, pe_before, Lorbit, Lspin) + end if write(*, *) " *************** Main Loop *************** " @@ -177,7 +178,10 @@ program swiftest_symba ldiscard_pl = ldiscard_pl .or. lfrag_add if (ldiscard_pl .or. ldiscard_tp) then - if (param%lenergy) call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_before, ke_spin_before, pe_before, Eorbit_before, Ltot) + if (param%lenergy) then + call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_before, ke_spin_before, pe_before, Lorbit, Lspin) + Eorbit_before = ke_orbit_before + ke_spin_before + pe_before + end if call symba_rearray(t, npl, nplm, ntp, nsppl, nsptp, symba_plA, symba_tpA, nmergeadd, mergeadd_list, discard_plA, & discard_tpA, ldiscard_pl, ldiscard_tp, mtiny, param, discard_l_pl, discard_stat_list) call io_discard_write_symba(t, mtiny, npl, nsppl, nsptp, nmergesub, symba_plA, & @@ -189,7 +193,8 @@ program swiftest_symba nplm = count(symba_plA%helio%swiftest%mass(1:npl) > mtiny) if (param%lenergy) then - call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_after, ke_spin_after, pe_after, Eorbit_after, Ltot) + call symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit_after, ke_spin_after, pe_after, Lorbit, Lspin) + Eorbit_after = ke_orbit_after + ke_spin_after + pe_after Ecollision = Eorbit_after - Eorbit_before ! Energy change resulting in this collisional event Total running energy offset from collision in this step symba_plA%helio%swiftest%Ecollisions = symba_plA%helio%swiftest%Ecollisions + Ecollision end if diff --git a/src/modules/module_interfaces.f90 b/src/modules/module_interfaces.f90 index 1481104cb..b0d36bd2d 100644 --- a/src/modules/module_interfaces.f90 +++ b/src/modules/module_interfaces.f90 @@ -1260,7 +1260,7 @@ END SUBROUTINE symba_user_getacch_tp END INTERFACE INTERFACE - SUBROUTINE symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe, te, Ltot) + SUBROUTINE symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe, Lorbit, Lspin) USE swiftest_globals USE swiftest_data_structures use module_symba @@ -1268,8 +1268,8 @@ SUBROUTINE symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe integer(I4B), intent(in) :: npl type(symba_pl), intent(inout) :: symba_plA real(DP), intent(in) :: j2rp2, j4rp4 - real(DP), intent(out) :: ke_orbit, ke_spin, pe, te - real(DP), dimension(:), intent(out) :: Ltot + real(DP), intent(out) :: ke_orbit, ke_spin, pe + real(DP), dimension(:), intent(out) :: Lorbit, Lspin END SUBROUTINE symba_energy_eucl END INTERFACE diff --git a/src/symba/symba_energy_eucl.f90 b/src/symba/symba_energy_eucl.f90 index 3dc691295..2d40016e5 100644 --- a/src/symba/symba_energy_eucl.f90 +++ b/src/symba/symba_energy_eucl.f90 @@ -1,4 +1,4 @@ -subroutine symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe, te, Ltot) +subroutine symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe, Lorbit, Lspin) !! author: David A. Minton !! !! Compute total system angular momentum vector and kinetic, potential and total system energy @@ -14,8 +14,8 @@ subroutine symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe ! arguments integer(I4B), intent(in) :: npl real(DP), intent(in) :: j2rp2, j4rp4 - real(DP), intent(out) :: ke_orbit, ke_spin, pe, te - real(DP), dimension(:), intent(out) :: Ltot + real(DP), intent(out) :: ke_orbit, ke_spin, pe + real(DP), dimension(:), intent(out) :: Lorbit, Lspin type(symba_pl), intent(inout) :: symba_plA ! internals @@ -23,22 +23,27 @@ subroutine symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe integer(I8B) :: k real(DP) :: rmag, v2, rot2, oblpot, hx, hy, hz, hsx, hsy, hsz real(DP), dimension(npl) :: irh, kepl, kespinpl, pecb - real(DP), dimension(npl) :: Lplx, Lply, Lplz + real(DP), dimension(npl) :: Lplorbitx, Lplorbity, Lplorbitz + real(DP), dimension(npl) :: Lplspinx, Lplspiny, Lplspinz real(DP), dimension(symba_plA%helio%swiftest%num_plpl_comparisons) :: pepl logical, dimension(symba_plA%helio%swiftest%num_plpl_comparisons) :: lstatpl logical, dimension(npl) :: lstatus ! executable code - Ltot = 0.0_DP + Lorbit(:) = 0.0_DP + Lspin(:) = 0.0_DP ke_orbit = 0.0_DP ke_spin = 0.0_DP associate(xb => symba_plA%helio%swiftest%xb, vb => symba_plA%helio%swiftest%vb, mass => symba_plA%helio%swiftest%mass, radius => symba_plA%helio%swiftest%radius, & Ip => symba_plA%helio%swiftest%Ip, rot => symba_plA%helio%swiftest%rot, xh => symba_plA%helio%swiftest%xh, status => symba_plA%helio%swiftest%status) kepl(:) = 0.0_DP - Lplx(:) = 0.0_DP - Lply(:) = 0.0_DP - Lplz(:) = 0.0_DP + Lplorbitx(:) = 0.0_DP + Lplorbity(:) = 0.0_DP + Lplorbitz(:) = 0.0_DP + Lplspinx(:) = 0.0_DP + Lplspiny(:) = 0.0_DP + Lplspinz(:) = 0.0_DP lstatus(1:npl) = status(1:npl) /= INACTIVE !!$omp simd private(v2, rot2, hx, hy, hz) do i = 1, npl @@ -54,26 +59,24 @@ subroutine symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe hsz = Ip(3,i) * radius(i)**2 * rot(3,i) ! Angular momentum from orbit and spin - Lplx(i) = mass(i) * (hsx + hx) - Lply(i) = mass(i) * (hsy + hy) - Lplz(i) = mass(i) * (hsz + hz) + Lplorbitx(i) = mass(i) * hx + Lplorbity(i) = mass(i) * hy + Lplorbitz(i) = mass(i) * hz + + Lplspinx(i) = mass(i) * hsx + Lplspiny(i) = mass(i) * hsy + Lplspinz(i) = mass(i) * hsz ! Kinetic energy from orbit and spin kepl(i) = mass(i) * v2 kespinpl(i) = mass(i) * Ip(3, i) * radius(i)**2 * rot2 end do - ke_orbit = 0.5_DP * sum(kepl(1:npl), lstatus(:)) - ke_spin = 0.5_DP * sum(kespinpl(1:npl), lstatus(:)) - Ltot(1) = sum(Lplx(1:npl), lstatus(1:npl)) - Ltot(2) = sum(Lply(1:npl), lstatus(1:npl)) - Ltot(3) = sum(Lplz(1:npl), lstatus(1:npl)) - ! Do the central body potential energy component first !$omp simd do i = 2, npl associate(px => xh(1,i), py => xh(2,i), pz => xh(3,i)) - pecb(i) = - mass(1) * mass(i) / sqrt(px**2 + py**2 + pz**2) + pecb(i) = -mass(1) * mass(i) / sqrt(px**2 + py**2 + pz**2) end associate end do @@ -85,6 +88,9 @@ subroutine symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe end associate end do + ke_orbit = 0.5_DP * sum(kepl(1:npl), lstatus(:)) + ke_spin = 0.5_DP * sum(kespinpl(1:npl), lstatus(:)) + pe = sum(pepl(:), lstatpl(:)) + sum(pecb(2:npl), lstatus(2:npl)) ! Potential energy from the oblateness term @@ -97,7 +103,14 @@ subroutine symba_energy_eucl(npl, symba_plA, j2rp2, j4rp4, ke_orbit, ke_spin, pe pe = pe + oblpot end if - te = ke_orbit + ke_spin + pe + Lorbit(1) = sum(Lplorbitx(1:npl), lstatus(1:npl)) + Lorbit(2) = sum(Lplorbity(1:npl), lstatus(1:npl)) + Lorbit(3) = sum(Lplorbitz(1:npl), lstatus(1:npl)) + + Lspin(1) = sum(Lplspinx(1:npl), lstatus(1:npl)) + Lspin(2) = sum(Lplspiny(1:npl), lstatus(1:npl)) + Lspin(3) = sum(Lplspinz(1:npl), lstatus(1:npl)) + end associate return diff --git a/src/symba/symba_frag_pos.f90 b/src/symba/symba_frag_pos.f90 index c61e61312..01c729c93 100644 --- a/src/symba/symba_frag_pos.f90 +++ b/src/symba/symba_frag_pos.f90 @@ -124,8 +124,8 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi real(DP), dimension(:), allocatable :: rmag, v_r_mag, v_t_mag real(DP), dimension(NDIM) :: Ltot_before real(DP), dimension(NDIM) :: Ltot_after - real(DP) :: Etot_before, ke_orb_before, ke_spin_before, pe_before, Lmag_before - real(DP) :: Etot_after, ke_orb_after, ke_spin_after, pe_after, Lmag_after, dEtot, dLmag + real(DP) :: Etot_before, ke_orbit_before, ke_spin_before, pe_before, Lmag_before + real(DP) :: Etot_after, ke_orbit_after, ke_spin_after, pe_after, Lmag_after, dEtot, dLmag real(DP), dimension(NDIM) :: L_frag_spin, L_frag_tot, L_frag_orb, L_offset real(DP) :: ke_family, ke_target, ke_frag, ke_offset real(DP), dimension(NDIM) :: x_col_unit, y_col_unit, z_col_unit @@ -194,9 +194,9 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi !write(*, "(' -------------------------------------------------------------------------------------')") !write(*, "(' | T_orb T_spin T pe Etot Ltot')") !write(*, "(' -------------------------------------------------------------------------------------')") - !write(*,fmtlabel) ' change |',(ke_orb_after - ke_orb_before) / abs(Etot_before), & + !write(*,fmtlabel) ' change |',(ke_orbit_after - ke_orbit_before) / abs(Etot_before), & ! (ke_spin_after - ke_spin_before)/ abs(Etot_before), & - ! (ke_orb_after + ke_spin_after - ke_orb_before - ke_spin_before)/ abs(Etot_before), & + ! (ke_orbit_after + ke_spin_after - ke_orbit_before - ke_spin_before)/ abs(Etot_before), & ! (pe_after - pe_before) / abs(Etot_before), & ! dEtot, dLmag ! @@ -230,9 +230,9 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi !write(*, "(' -------------------------------------------------------------------------------------')") !write(*, "(' | T_orb T_spin T pe Etot Ltot')") !write(*, "(' -------------------------------------------------------------------------------------')") - !write(*,fmtlabel) ' change |',(ke_orb_after - ke_orb_before) / abs(Etot_before), & + !write(*,fmtlabel) ' change |',(ke_orbit_after - ke_orbit_before) / abs(Etot_before), & ! (ke_spin_after - ke_spin_before)/ abs(Etot_before), & - ! (ke_orb_after + ke_spin_after - ke_orb_before - ke_spin_before)/ abs(Etot_before), & + ! (ke_orbit_after + ke_spin_after - ke_orbit_before - ke_spin_before)/ abs(Etot_before), & ! (pe_after - pe_before) / abs(Etot_before), & ! dEtot, dLmag !else @@ -243,21 +243,21 @@ subroutine symba_frag_pos(param, symba_plA, family, x, v, L_spin, Ip, mass, radi call restructure_failed_fragments() try = try + 1 end do - !write(*, "(' -------------------------------------------------------------------------------------')") - !write(*, "(' Final diagnostic')") - !write(*, "(' -------------------------------------------------------------------------------------')") + write(*, "(' -------------------------------------------------------------------------------------')") + write(*, "(' Final diagnostic')") + write(*, "(' -------------------------------------------------------------------------------------')") if (lmerge) then write(*,*) "symba_frag_pos failed after: ",try," tries" else write(*,*) "symba_frag_pos succeeded after: ",try," tries" - ! write(*, "(' dL_tot should be very small' )") - ! write(*,fmtlabel) ' dL_tot |', dLmag - ! write(*, "(' dE_tot should be negative and equal to Qloss' )") - ! write(*,fmtlabel) ' dE_tot |', dEtot - ! write(*,fmtlabel) ' Qloss |', -Qloss / abs(Etot_before) - ! write(*,fmtlabel) ' dE - Qloss |', (Etot_after - Etot_before + Qloss) / abs(Etot_before) + write(*, "(' dL_tot should be very small' )") + write(*,fmtlabel) ' dL_tot |', dLmag + write(*, "(' dE_tot should be negative and equal to Qloss' )") + write(*,fmtlabel) ' dE_tot |', dEtot + write(*,fmtlabel) ' Qloss |', -Qloss / abs(Etot_before) + write(*,fmtlabel) ' dE - Qloss |', (Etot_after - Etot_before + Qloss) / abs(Etot_before) end if - !write(*, "(' -------------------------------------------------------------------------------------')") + write(*, "(' -------------------------------------------------------------------------------------')") call restore_scale_factors() @@ -416,12 +416,11 @@ subroutine calculate_system_energy(linclude_fragments) ! Arguments logical, intent(in) :: linclude_fragments ! Internals - real(DP) :: ke_orb, ke_spin, pe - real(DP), dimension(NDIM) :: Ltot + real(DP) :: ke_orbit, ke_spin, pe, te + real(DP), dimension(NDIM) :: Ltot, Lorbit, Lspin integer(I4B) :: i, npl_new, nplm logical, dimension(:), allocatable :: ltmp logical :: lk_plpl - real(DP) :: te type(symba_pl) :: symba_plwksp ! Because we're making a copy of symba_pl with the excludes/fragments appended, we need to deallocate the @@ -477,7 +476,7 @@ subroutine calculate_system_energy(linclude_fragments) nplm = count(symba_plwksp%helio%swiftest%mass > param%mtiny / mscale) call util_dist_index_plpl(npl_new, nplm, symba_plwksp) - call symba_energy_eucl(npl_new, symba_plwksp, param%j2rp2, param%j4rp4, ke_orb, ke_spin, pe, te, Ltot) + call symba_energy_eucl(npl_new, symba_plwksp, param%j2rp2, param%j4rp4, ke_orbit, ke_spin, pe, Lorbit, Lspin) ! Restore the big array deallocate(symba_plwksp%helio%swiftest%k_plpl) @@ -486,12 +485,12 @@ subroutine calculate_system_energy(linclude_fragments) ! Calculate the current fragment energy and momentum balances if (linclude_fragments) then - Ltot_after(:) = Ltot(:) - Lmag_after = norm2(Ltot(:)) - ke_orb_after = ke_orb + Ltot_after(:) = Lorbit(:) + Lspin(:) + Lmag_after = norm2(Ltot_after(:)) + ke_orbit_after = ke_orbit ke_spin_after = ke_spin pe_after = pe - Etot_after = ke_orb_after + ke_spin_after + pe_after + Etot_after = ke_orbit_after + ke_spin_after + pe_after dEtot = (Etot_after - Etot_before) / abs(Etot_before) dLmag = norm2(Ltot_after(:) - Ltot_before(:)) / Lmag_before ke_frag = 0._DP @@ -502,12 +501,12 @@ subroutine calculate_system_energy(linclude_fragments) ke_offset = ke_frag - ke_target L_offset(:) = Ltot_before(:) - Ltot(:) else - Ltot_before(:) = Ltot(:) - Lmag_before = norm2(Ltot(:)) - ke_orb_before = ke_orb + Ltot_before(:) = Lorbit(:) + Lspin(:) + Lmag_before = norm2(Ltot_before(:)) + ke_orbit_before = ke_orbit ke_spin_before = ke_spin pe_before = pe - Etot_before = ke_orb_before + ke_spin_before + pe_before + Etot_before = ke_orbit_before + ke_spin_before + pe_before end if end associate return @@ -673,7 +672,7 @@ subroutine set_fragment_radial_velocities(lmerge) real(DP), dimension(:), allocatable :: v_r_initial, v_r_sigma real(DP), dimension(:,:), allocatable :: v_r - ! Set the "target" ke_orb_after (the value of the orbital kinetic energy that the fragments ought to have) + ! Set the "target" ke_orbit_after (the value of the orbital kinetic energy that the fragments ought to have) if ((ke_target < 0.0_DP) .or. (ke_offset > 0.0_DP)) then !if (ke_target < 0.0_DP) write(*,*) 'Negative ke_target: ',ke_target