From 893407ad61626d03cdedd3bcd7511c8cb79a66d1 Mon Sep 17 00:00:00 2001 From: David Minton Date: Thu, 1 Apr 2021 11:14:31 -0400 Subject: [PATCH] Reformatted code to pull use statements out of individual subroutines in the submodules due to complains from ifort. Also cleaned up some cruft --- src/discard/discard.f90 | 78 ++++-- src/driftkick/drift.f90 | 12 +- src/driftkick/kick.f90 | 50 ++-- src/eucl/eucl.f90 | 4 +- src/gr/gr.f90 | 3 +- src/helio/helio_coord.f90 | 6 +- src/helio/helio_drift.f90 | 104 ++++--- src/helio/helio_getacch.f90 | 243 ++++++++-------- src/helio/helio_setup.f90 | 3 +- src/helio/helio_step.f90 | 4 +- src/io/io.f90 | 2 +- src/modules/rmvs_classes.f90 | 4 +- src/modules/swiftest_classes.f90 | 16 -- src/obl/obl.f90 | 2 +- src/operators/operator_cross.f90 | 3 +- src/orbel/orbel_el2xv.f90 | 4 +- src/orbel/orbel_scget.f90 | 2 +- src/orbel/orbel_xv2aeq.f90 | 2 +- src/orbel/orbel_xv2aqt.f90 | 2 +- src/orbel/orbel_xv2el.f90 | 3 +- src/rmvs/rmvs_discard.f90 | 9 +- src/rmvs/rmvs_encounter_check.f90 | 3 +- src/rmvs/rmvs_getacch.f90 | 3 +- src/rmvs/rmvs_interp.f90 | 3 +- src/rmvs/rmvs_obl.f90 | 4 +- src/rmvs/rmvs_peri.f90 | 2 +- src/rmvs/rmvs_setup.f90 | 5 +- src/rmvs/rmvs_spill_and_fill.f90 | 177 ++++++------ src/rmvs/rmvs_step.f90 | 451 +++++++++++++++--------------- src/setup/setup.f90 | 10 +- src/step/step.f90 | 3 +- src/user/user_getacch.f90 | 26 +- src/util/util_coord.f90 | 7 +- src/util/util_copy.f90 | 9 +- src/util/util_exit.f90 | 2 +- src/util/util_hills.f90 | 2 +- src/util/util_index.f90 | 2 +- src/util/util_peri.f90 | 2 +- src/util/util_reverse_status.f90 | 3 +- src/util/util_sort_dp.f90 | 2 +- src/util/util_sort_i4b.f90 | 2 +- src/util/util_sort_sp.f90 | 2 +- src/util/util_spill_and_fill.f90 | 62 ++-- src/util/util_toupper.f90 | 2 +- src/util/util_valid.f90 | 2 +- src/util/util_version.f90 | 86 +++--- src/whm/whm_coord.f90 | 4 +- src/whm/whm_drift.f90 | 3 +- src/whm/whm_getacch.f90 | 78 +++--- src/whm/whm_gr.f90 | 17 +- src/whm/whm_setup.f90 | 8 +- src/whm/whm_spill_and_fill.f90 | 3 +- src/whm/whm_step.f90 | 4 +- 53 files changed, 719 insertions(+), 826 deletions(-) diff --git a/src/discard/discard.f90 b/src/discard/discard.f90 index dcb3099d4..f46f802d0 100644 --- a/src/discard/discard.f90 +++ b/src/discard/discard.f90 @@ -1,23 +1,23 @@ submodule (swiftest_classes) s_discard + use swiftest contains - module procedure discard_system + module subroutine discard_system(self, config) !! author: David A. Minton !! !! Check to see if particles should be discarded based on their positions relative to the massive bodies !! !! Adapted from David E. Kaufmann's Swifter routine: discard.f90 !! Adapted from Hal Levison's Swift routine discard.f - use swiftest implicit none - logical, dimension(:), allocatable :: lspill_list + class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object + class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters - real(DP) :: msys if (self%tp%nbody == 0) return select type(self) class is (whm_nbody_system) associate(cb => self%cb, pl => self%pl, tp => self%tp, t => config%t, dt => config%dt, & msys => self%msys, discards => self%tp_discards, & - ntp => self%tp%nbody) + ntp => self%tp%nbody, ldiscard => self%tp%ldiscard) if ((config%rmin >= 0.0_DP) .or. (config%rmax >= 0.0_DP) .or. & (config%rmaxu >= 0.0_DP) .or. ((config%qmin >= 0.0_DP) .and. (config%qmin_coord == "BARY"))) then call pl%h2b(cb) @@ -27,21 +27,19 @@ if (ntp > 0) call tp%discard_sun(cb, config, t, msys) end if if (config%qmin >= 0.0_DP .and. ntp > 0) call tp%discard_peri(cb, pl, config, t, msys) - if (config%lclose .and. ntp > 0) call tp%discard_pl(cb, pl, config, t, dt) + if (config%lclose .and. ntp > 0) call tp%discard_pl(pl, t, dt) if (any(tp%ldiscard(1:ntp))) then ! Spill the discards to the spill list - allocate(lspill_list, source = tp%ldiscard) - call tp%spill(discards, lspill_list) + call tp%spill(discards, ldiscard) call self%write_discard(config, discards) - deallocate(lspill_list) end if end associate end select return - end procedure discard_system + end subroutine discard_system - module procedure discard_sun_tp + module subroutine discard_sun_tp(self, cb, config, t, msys) !! author: David A. Minton !! !! Check to see if test particles should be discarded based on their positions relative to the Sun @@ -49,7 +47,14 @@ !! !! Adapted from David E. Kaufmann's Swifter routine: discard_sun.f90 !! Adapted from Hal Levison's Swift routine discard_sun.f - use swiftest + implicit none + ! Arguments + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object + class(swiftest_configuration), intent(in) :: config !! configuration parameters + real(DP), intent(in) :: t !! Current simulation tim + real(DP), intent(in) :: msys !! Total system mass + ! Internals integer(I4B) :: i real(DP) :: energy, vb2, rb2, rh2, rmin2, rmax2, rmaxu2 @@ -83,23 +88,29 @@ end associate return - - end procedure discard_sun_tp + end subroutine discard_sun_tp - module procedure discard_peri_tp + module subroutine discard_peri_tp(self, cb, pl, config, t, msys) !! author: David A. Minton !! !! Check to see if a test particle should be discarded because its perihelion distance becomes too small !! !! Adapted from David E. Kaufmann's Swifter routine: discard_peri.f90 !! Adapted from Hal Levison's Swift routine discard_peri.f - use swiftest + implicit none + ! Arguments + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object + class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object + class(swiftest_configuration), intent(in) :: config !! configuration parameters + real(DP), intent(in) :: t !! Current simulation tim + real(DP), intent(in) :: msys !! Total system mass + ! Internals logical, save :: lfirst = .true. - integer(I4B) :: i, j, ih, ntp, npl + integer(I4B) :: i, j, ih real(DP) :: r2 real(DP), dimension(NDIM) :: dx - associate(tp => self, ntp => self%nbody, npl => pl%nbody, qmin_coord => config%qmin_coord) if (lfirst) then call util_hills(npl, pl) @@ -132,19 +143,23 @@ end associate return - end procedure discard_peri_tp + end subroutine discard_peri_tp - module procedure discard_pl_tp + module subroutine discard_pl_tp(self, pl, t, dt) !! author: David A. Minton !! !! Check to see if test particles should be discarded based on their positions relative to the massive bodies !! !! Adapted from David E. Kaufmann's Swifter routine: discard_pl.f90 !! Adapted from Hal Levison's Swift routine discard_pl.f - use swiftest implicit none - - integer(I4B) :: i, j, isp, ntp, npl + ! Arguments + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + class(swiftest_pl), intent(inout) :: pl !! Swiftest massive body object + real(DP), intent(in) :: t !! Current simulation tim + real(DP), intent(in) :: dt !! Stepsize + ! Internals + integer(I4B) :: i, j, isp real(DP) :: r2min, radius real(DP), dimension(NDIM) :: dx, dv @@ -170,9 +185,9 @@ end associate return - end procedure discard_pl_tp + end subroutine discard_pl_tp - module procedure discard_pl_close + module subroutine discard_pl_close(dx, dv, dt, r2crit, iflag, r2min) !! author: David A. Minton !! !! Check to see if a test particle and massive body are having, or will have within the next time step, an encounter such @@ -180,9 +195,15 @@ !! !! Adapted from David E. Kaufmann's Swifter routine: discard_pl_close.f90 !! Adapted from Hal Levison's Swift routine discard_pl_close.f - use swiftest + implicit none + ! Arguments + real(DP), dimension(:), intent(in) :: dx, dv + real(DP), intent(in) :: dt, r2crit + integer(I4B), intent(out) :: iflag + real(DP), intent(out) :: r2min + ! Internals real(DP) :: r2, v2, vdotr, tmin - + r2 = dot_product(dx(:), dx(:)) if (r2 <= r2crit) then iflag = 1 @@ -208,7 +229,6 @@ end if return - - end procedure discard_pl_close + end subroutine discard_pl_close end submodule s_discard diff --git a/src/driftkick/drift.f90 b/src/driftkick/drift.f90 index c70f3bb53..6a175c6fb 100644 --- a/src/driftkick/drift.f90 +++ b/src/driftkick/drift.f90 @@ -1,5 +1,5 @@ submodule (swiftest_classes) drift_implementation - + use swiftest !> Integration control parameters: real(DP), parameter :: E2MAX = 0.36_DP real(DP), parameter :: DM2MAX = 0.16_DP @@ -16,7 +16,6 @@ module pure subroutine drift_one(mu, x, v, dt, iflag) !! !! Adapted from David E. Kaufmann's Swifter routine routine drift_one.f90 !! Adapted from Hal Levison and Martin Duncan's Swift routine drift_one.f - use swiftest_globals implicit none ! Arguments real(DP), intent(in) :: mu !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body to drift @@ -46,7 +45,6 @@ pure subroutine drift_dan(mu, x0, v0, dt0, iflag) !! !! Adapted from David E. Kaufmann's Swifter routine: drift_dan.f90 !! Adapted from Hal Levison and Martin Duncan's Swift routine drift_dan.f - use swiftest implicit none integer(I4B), intent(out) :: iflag real(DP), intent(in) :: mu, dt0 @@ -118,7 +116,6 @@ pure subroutine drift_kepmd(dm, es, ec, x, s, c) !! !! Adapted from David E. Kaufmann's Swifter routine: drift_kepmd.f90 !! Adapted from Martin Duncan's Swift routine drift_kepmd.f - use swiftest implicit none real(DP), intent(in) :: dm, es, ec real(DP), intent(out) :: x, s, c @@ -155,7 +152,6 @@ pure subroutine drift_kepu(dt,r0,mu,alpha,u,fp,c1,c2,c3,iflag) !! !! Adapted from David E. Kaufmann's Swifter routine: drift_kepu.f90 !! Adapted from Hal Levison's Swift routine drift_kepu.f - use swiftest implicit none integer(I4B), intent(out) :: iflag real(DP), intent(in) :: dt, r0, mu, alpha, u @@ -183,7 +179,6 @@ pure subroutine drift_kepu_fchk(dt, r0, mu, alpha, u, s, f) !! !! Adapted from David E. Kaufmann's Swifter routine: drift_kepu_fchk.f90 !! Adapted from Martin Duncan's Swift routine drift_kepu_fchk.f - use swiftest implicit none real(DP), intent(in) :: dt, r0, mu, alpha, u, s real(DP), intent(out) :: f @@ -206,7 +201,6 @@ pure subroutine drift_kepu_guess(dt, r0, mu, alpha, u, s) !! !! Adapted from David E. Kaufmann's Swifter routine: drift_kepu_guess.f90 !! Adapted from Hal Levison and Martin Duncan's Swift routine drift_kepu_guess.f - use swiftest implicit none real(DP), intent(in) :: dt, r0, mu, alpha, u real(DP), intent(out) :: s @@ -245,7 +239,6 @@ pure subroutine drift_kepu_lag(s, dt, r0, mu, alpha, u, fp, c1, c2, c3, iflag) !! !! Adapted from David E. Kaufmann's Swifter routine: drift_kepu_lag.f90 !! Adapted from Hal Levison's Swift routine drift_kepu_lag.f - use swiftest implicit none integer(I4B), intent(out) :: iflag real(DP), intent(in) :: dt, r0, mu, alpha, u @@ -290,7 +283,6 @@ pure subroutine drift_kepu_new(s, dt, r0, mu, alpha, u, fp, c1, c2, c3, iflag) !! !! Adapted from David E. Kaufmann's Swifter routine: drift_kepu_new.f90 !! Adapted from Hal Levison's Swift routine drift_kepu_new.f - use swiftest implicit none integer(I4B), intent(out) :: iflag real(DP), intent(in) :: dt, r0, mu, alpha, u @@ -332,7 +324,6 @@ pure subroutine drift_kepu_p3solve(dt, r0, mu, alpha, u, s, iflag) !! !! Adapted from David E. Kaufmann's Swifter routine: drift_kepu_p3solve.f90 !! Adapted from Martin Duncan's Swift routine drift_kepu_p3solve.f - use swiftest implicit none integer(I4B), intent(out) :: iflag real(DP), intent(in) :: dt, r0, mu, alpha, u @@ -376,7 +367,6 @@ pure subroutine drift_kepu_stumpff(x, c0, c1, c2, c3) !! !! Adapted from David E. Kaufmann's Swifter routine: drift_kepu_stumpff.f90 !! Adapted from Hal Levison's Swift routine drift_kepu_stumpff.f - use swiftest implicit none real(DP), intent(inout) :: x real(DP), intent(out) :: c0, c1, c2, c3 diff --git a/src/driftkick/kick.f90 b/src/driftkick/kick.f90 index b4056afb3..32888c9a2 100644 --- a/src/driftkick/kick.f90 +++ b/src/driftkick/kick.f90 @@ -1,36 +1,42 @@ submodule(swiftest_classes) s_kick -contains - module procedure kick_vh_body - !! author: David A. Minton - !! - !! Kick heliocentric velocities of bodies - !! - !! Adapted from Martin Duncan and Hal Levison's Swift routine kickvh.f and kickvh_tp.f - !! Adapted from David E. Kaufmann's Swifter routine whm_kickvh.f90 and whm_kickvh_tp.f90 use swiftest - implicit none - integer(I4B) :: i - - associate(n => self%nbody, vh => self%vh, ah => self%ah, status => self%status) - if (n == 0) return - do concurrent(i = 1:n, status(i) == ACTIVE) - vh(:, i) = vh(:, i) + ah(:, i) * dt - end do - end associate +contains + module subroutine kick_vh_body(self, dt) + !! author: David A. Minton + !! + !! Kick heliocentric velocities of bodies + !! + !! Adapted from Martin Duncan and Hal Levison's Swift routine kickvh.f and kickvh_tp.f + !! Adapted from David E. Kaufmann's Swifter routine whm_kickvh.f90 and whm_kickvh_tp.f90 + implicit none + ! Arguments + class(swiftest_body), intent(inout) :: self !! Swiftest generic body object + real(DP), intent(in) :: dt !! Stepsize + ! Internals + integer(I4B) :: i - return + associate(n => self%nbody, vh => self%vh, ah => self%ah, status => self%status) + if (n == 0) return + do concurrent(i = 1:n, status(i) == ACTIVE) + vh(:, i) = vh(:, i) + ah(:, i) * dt + end do + end associate - end procedure kick_vh_body + return + end subroutine kick_vh_body - module procedure kick_vb_body + module subroutine kick_vb_body(self, dt) !! author: David A. Minton !! !! Kick barycentric velocities of bodies !! !! Adapted from Martin Duncan and Hal Levison's Swift routine kickvh.f and kickvh_tp.f !! Adapted from David E. Kaufmann's Swifter routine helio_kickvb.f90 and helio_kickvb_tp.f90 - use swiftest implicit none + ! Arguments + class(swiftest_body), intent(inout) :: self !! Swiftest generic body object + real(DP), intent(in) :: dt !! Stepsize + ! Internals integer(I4B) :: i associate(n => self%nbody, vb => self%vb, ah => self%ah, status => self%status) @@ -42,5 +48,5 @@ return - end procedure kick_vb_body + end subroutine kick_vb_body end submodule s_kick diff --git a/src/eucl/eucl.f90 b/src/eucl/eucl.f90 index 4e06ca8a2..49a72319c 100644 --- a/src/eucl/eucl.f90 +++ b/src/eucl/eucl.f90 @@ -1,6 +1,6 @@ submodule (swiftest_classes) s_eucl + use swiftest contains - module procedure eucl_dist_index_plpl !! author: Jacob R. Elliott and David A. Minton !! @@ -10,7 +10,6 @@ !! !! Mélodie Angeletti, Jean-Marie Bonny, Jonas Koko. Parallel Euclidean distance matrix computation on big datasets *. !! 2019. hal-0204751 - use swiftest implicit none integer(I4B) :: i, j, k, kp, p @@ -41,7 +40,6 @@ !! author: Jacob R. Elliott and David A. Minton !! !! Efficient parallel loop-blocking algrorithm for evaluating the Euclidean distance matrix for planet-planet - use swiftest implicit none integer(I4B) :: k, i, j real(DP), dimension(NDIM) :: dx diff --git a/src/gr/gr.f90 b/src/gr/gr.f90 index edfa205f0..134038d5a 100644 --- a/src/gr/gr.f90 +++ b/src/gr/gr.f90 @@ -1,6 +1,6 @@ submodule(swiftest_classes) s_gr + use swiftest contains - module procedure gr_getaccb_ns_body !! author: David A. Minton !! @@ -8,7 +8,6 @@ !! Based on Quinn, Tremaine, & Duncan 1998 !! !! Adapted from David A. Minton's Swifter routine routine gr_getaccb_ns.f90 - use swiftest implicit none real(DP), dimension(NDIM) :: xh, vh diff --git a/src/helio/helio_coord.f90 b/src/helio/helio_coord.f90 index b95308a92..e14ea4612 100644 --- a/src/helio/helio_coord.f90 +++ b/src/helio/helio_coord.f90 @@ -1,6 +1,6 @@ submodule (helio_classes) s_helio_coord + use swiftest contains - module subroutine helio_coord_vb2vh_pl(self, cb) !! author: David A. Minton !! @@ -8,7 +8,6 @@ module subroutine helio_coord_vb2vh_pl(self, cb) !! !! Adapted from David E. Kaufmann's Swifter routine coord_vb2vh.f90 !! Adapted from Hal Levison's Swift routine coord_vb2vh.f - use swiftest implicit none ! Arguments class(helio_pl), intent(inout) :: self !! Helio massive body object @@ -34,7 +33,6 @@ module subroutine helio_coord_vb2vh_tp(self, vbcb) !! !! Adapted from David E. Kaufmann's Swifter routine coord_vb2vh_tp.f90 !! Adapted from Hal Levison's Swift routine coord_vb2h_tp.f - use swiftest implicit none ! Arguments class(helio_tp), intent(inout) :: self !! Helio massive body object @@ -58,7 +56,6 @@ module subroutine helio_coord_vh2vb_pl(self, cb) !! !! Adapted from David E. Kaufmann's Swifter routine coord_vh2vb.f90 !! Adapted from Hal Levison's Swift routine coord_vh2b.f - use swiftest implicit none ! Arguments class(helio_pl), intent(inout) :: self !! Helio massive body object @@ -86,7 +83,6 @@ module subroutine helio_coord_vh2vb_tp(self, vbcb) !! !! Adapted from David E. Kaufmann's Swifter routine coord_vh2vb_tp.f90 !! Adapted from Hal Levison's Swift routine coord_vh2b_tp.f - use swiftest implicit none ! Arguments class(helio_tp), intent(inout) :: self !! Helio massive body object diff --git a/src/helio/helio_drift.f90 b/src/helio/helio_drift.f90 index ba7b798ae..244dd5f9e 100644 --- a/src/helio/helio_drift.f90 +++ b/src/helio/helio_drift.f90 @@ -1,60 +1,58 @@ submodule (helio_classes) s_helio_drift + use swiftest contains module subroutine helio_drift_pl(self, cb, config, dt) - !! author: David A. Minton - !! - !! Loop through massive bodies and call Danby drift routine - !! New vectorized version included - !! - !! Adapted from David E. Kaufmann's Swifter routine helio_drift.f90 - !! Adapted from Hal Levison's Swift routine drift.f - use swiftest - implicit none - ! Arguments - class(helio_pl), intent(inout) :: self !! Helio test particle data structure - class(swiftest_cb), intent(inout) :: cb !! Helio central body particle data structuree - class(swiftest_configuration), intent(in) :: config !! Input collection of - real(DP), intent(in) :: dt !! Stepsize - ! Internals - integer(I4B), dimension(:),allocatable :: iflag !! Vectorized error code flag - integer(I4B) :: i !! Loop counter - real(DP) :: rmag, vmag2, energy, dtp + !! author: David A. Minton + !! + !! Loop through massive bodies and call Danby drift routine + !! New vectorized version included + !! + !! Adapted from David E. Kaufmann's Swifter routine helio_drift.f90 + !! Adapted from Hal Levison's Swift routine drift.f + implicit none + ! Arguments + class(helio_pl), intent(inout) :: self !! Helio test particle data structure + class(swiftest_cb), intent(inout) :: cb !! Helio central body particle data structuree + class(swiftest_configuration), intent(in) :: config !! Input collection of + real(DP), intent(in) :: dt !! Stepsize + ! Internals + integer(I4B), dimension(:),allocatable :: iflag !! Vectorized error code flag + integer(I4B) :: i !! Loop counter + real(DP) :: rmag, vmag2, energy, dtp - associate(npl => self%nbody, & - xh => self%xh, & - vb => self%vb, & - status => self%status,& - mu => self%mu) + associate(npl => self%nbody, & + xh => self%xh, & + vb => self%vb, & + status => self%status,& + mu => self%mu) - if (npl == 0) return + if (npl == 0) return - allocate(iflag(npl)) - iflag(:) = 0 - - do i = 1, npl - if (config%lgr) then - rmag = norm2(xh(:, i)) - vmag2 = dot_product(vb(:, i), vb(:, i)) - energy = 0.5_DP * vmag2 - mu(i) / rmag - dtp = dt * (1.0_DP + 3 * config%inv_c2 * energy) - else - dtp = dt - end if - call drift_one(mu(i), xh(:, i), vb(:, i), dtp, iflag(i)) - end do - if (any(iflag(1:npl) /= 0)) then + allocate(iflag(npl)) + iflag(:) = 0 + do i = 1, npl - write(*, *) " Planet ", self%name(i), " is lost!!!!!!!!!!" - write(*, *) xh - write(*, *) vb - write(*, *) " stopping " - call util_exit(FAILURE) - end do - end if - end associate - - return - + if (config%lgr) then + rmag = norm2(xh(:, i)) + vmag2 = dot_product(vb(:, i), vb(:, i)) + energy = 0.5_DP * vmag2 - mu(i) / rmag + dtp = dt * (1.0_DP + 3 * config%inv_c2 * energy) + else + dtp = dt + end if + call drift_one(mu(i), xh(:, i), vb(:, i), dtp, iflag(i)) + end do + if (any(iflag(1:npl) /= 0)) then + do i = 1, npl + write(*, *) " Planet ", self%name(i), " is lost!!!!!!!!!!" + write(*, *) xh + write(*, *) vb + write(*, *) " stopping " + call util_exit(FAILURE) + end do + end if + end associate + return end subroutine helio_drift_pl module subroutine helio_drift_linear_pl(self, cb, dt, pt) @@ -64,7 +62,6 @@ module subroutine helio_drift_linear_pl(self, cb, dt, pt) !! !! Adapted from David E. Kaufmann's Swifter routine helio_lindrift.f90 !! Adapted from Hal Levison's Swift routine helio_lindrift.f - use swiftest implicit none ! Arguments class(helio_pl), intent(inout) :: self !! Helio massive body object @@ -99,7 +96,6 @@ module subroutine helio_drift_tp(self, cb, config, dt) !! !! Adapted from David E. Kaufmann's Swifter routine helio_drift_tp.f90 !! Adapted from Hal Levison's Swift routine drift_tp.f - use swiftest implicit none ! Arguments class(helio_tp), intent(inout) :: self !! Helio test particle data structure @@ -141,8 +137,6 @@ module subroutine helio_drift_tp(self, cb, config, dt) end associate return - - end subroutine helio_drift_tp module subroutine helio_drift_linear_tp(self, dt, pt) @@ -153,7 +147,6 @@ module subroutine helio_drift_linear_tp(self, dt, pt) !! !! Adapted from David E. Kaufmann's Swifter routine helio_lindrift_tp.f90 !! Adapted from Hal Levison's Swift routine helio_lindrift_tp.f - use swiftest implicit none ! Arguments class(helio_tp), intent(inout) :: self !! Helio test particle data structure @@ -169,7 +162,6 @@ module subroutine helio_drift_linear_tp(self, dt, pt) end associate return - end subroutine helio_drift_linear_tp end submodule s_helio_drift diff --git a/src/helio/helio_getacch.f90 b/src/helio/helio_getacch.f90 index 17453bc30..0ac8de225 100644 --- a/src/helio/helio_getacch.f90 +++ b/src/helio/helio_getacch.f90 @@ -1,144 +1,141 @@ submodule (helio_classes) s_helio_getacch -contains -module subroutine helio_getacch_pl(self, cb, config, t) - !! author: David A. Minton - !! - !! Compute heliocentric accelerations of massive bodies - !! - !! Adapted from David E. Kaufmann's Swifter routine helio_getacch.f90 - !! Adapted from Hal Levison's Swift routine helio_getacch.f use swiftest - implicit none - ! Arguments - class(helio_pl), intent(inout) :: self !! Helio massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structure - class(swiftest_configuration), intent(in) :: config !! Input collection of - real(DP), intent(in) :: t !! Current time - ! Internals - logical, save :: lmalloc = .true. - integer(I4B) :: i - real(DP) :: r2 - real(DP), dimension(:), allocatable, save :: irh - real(DP), dimension(:, :), allocatable, save :: xh_loc, aobl - - associate(npl => self%nbody) - !if (lflag) then - self%ahi(:,2:npl) = 0.0_DP - call helio_getacch_int_pl(self, t) - !end if - !if (config%loblatecb) call self%obl_acc(cb) TODO: Fix this - !else - self%ah(:,:) = self%ahi(:,:) - !end if - if (config%lextra_force) call self%user_getacch(cb, config, t) - end associate - - return - end subroutine helio_getacch_pl - - module subroutine helio_getacch_tp(self, cb, pl, config, t, xh) +contains + module subroutine helio_getacch_pl(self, cb, config, t) !! author: David A. Minton !! - !! Compute heliocentric accelerations of test particles + !! Compute heliocentric accelerations of massive bodies !! - !! Adapted from David E. Kaufmann's Swifter routine helio_getacch_tp.f90 - !! Adapted from Hal Levison's Swift routine helio_getacch_tp.f - use swiftest + !! Adapted from David E. Kaufmann's Swifter routine helio_getacch.f90 + !! Adapted from Hal Levison's Swift routine helio_getacch.f implicit none ! Arguments - class(helio_tp), intent(inout) :: self !! Helio test particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree - class(whm_pl), intent(inout) :: pl !! WHM massive body particle data structure. - class(swiftest_configuration), intent(in) :: config !! Input collection of - real(DP), intent(in) :: t !! Current time - real(DP), dimension(:,:), intent(in) :: xh !! Heliocentric positions of planets + class(helio_pl), intent(inout) :: self !! Helio massive body particle data structure + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structure + class(swiftest_configuration), intent(in) :: config !! Input collection of + real(DP), intent(in) :: t !! Current time ! Internals - logical, save :: lmalloc = .true. - integer(I4B) :: i - real(DP) :: r2, mu - real(DP), dimension(:), allocatable, save :: irh, irht - real(DP), dimension(:, :), allocatable, save :: aobl, xht, aoblt - - ! executable code - associate(ntp => self%nbody, npl => pl%nbody) + logical, save :: lmalloc = .true. + integer(I4B) :: i + real(DP) :: r2 + real(DP), dimension(:), allocatable, save :: irh + real(DP), dimension(:, :), allocatable, save :: xh_loc, aobl + + associate(npl => self%nbody) !if (lflag) then - self%ahi(:,:) = 0.0_DP - call helio_getacch_int_tp(self, pl, t, xh) + self%ahi(:,2:npl) = 0.0_DP + call helio_getacch_int_pl(self, t) !end if !if (config%loblatecb) call self%obl_acc(cb) TODO: Fix this + !else self%ah(:,:) = self%ahi(:,:) + !end if if (config%lextra_force) call self%user_getacch(cb, config, t) - if (config%lgr) call self%gr_getacch(cb, config) end associate + return - end subroutine helio_getacch_tp + end subroutine helio_getacch_pl - module subroutine helio_getacch_int_pl(self, t) - !! author: David A. Minton - !! - !! Compute direct cross term heliocentric accelerations of massive bodiese - !! - !! Adapted from David E. Kaufmann's Swifter routine helio_getacch_int.f90 - !! Adapted from Hal Levison's Swift routine getacch_ah3.f - use swiftest - implicit none - ! Arguments - class(helio_pl), intent(inout) :: self !! Helio massive body particle data structure - real(DP), intent(in) :: t !! Current time - ! Internals - integer(I4B) :: i, j - real(DP) :: rji2, irij3, faci, facj - real(DP), dimension(NDIM) :: dx - - associate(npl => self%nbody) - do i = 2, npl - 1 - do j = i + 1, npl - dx(:) = self%xh(:,j) - self%xh(:,i) - rji2 = dot_product(dx(:), dx(:)) - irij3 = 1.0_DP / (rji2 * sqrt(rji2)) - faci = self%Gmass(i) * irij3 - facj = self%Gmass(j) * irij3 - self%ahi(:,i) = self%ahi(:,i) + facj * dx(:) - self%ahi(:,i) = self%ahi(:,j) - faci * dx(:) + module subroutine helio_getacch_tp(self, cb, pl, config, t, xh) + !! author: David A. Minton + !! + !! Compute heliocentric accelerations of test particles + !! + !! Adapted from David E. Kaufmann's Swifter routine helio_getacch_tp.f90 + !! Adapted from Hal Levison's Swift routine helio_getacch_tp.f + implicit none + ! Arguments + class(helio_tp), intent(inout) :: self !! Helio test particle data structure + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree + class(whm_pl), intent(inout) :: pl !! WHM massive body particle data structure. + class(swiftest_configuration), intent(in) :: config !! Input collection of + real(DP), intent(in) :: t !! Current time + real(DP), dimension(:,:), intent(in) :: xh !! Heliocentric positions of planets + ! Internals + logical, save :: lmalloc = .true. + integer(I4B) :: i + real(DP) :: r2, mu + real(DP), dimension(:), allocatable, save :: irh, irht + real(DP), dimension(:, :), allocatable, save :: aobl, xht, aoblt + + ! executable code + associate(ntp => self%nbody, npl => pl%nbody) + !if (lflag) then + self%ahi(:,:) = 0.0_DP + call helio_getacch_int_tp(self, pl, t, xh) + !end if + !if (config%loblatecb) call self%obl_acc(cb) TODO: Fix this + self%ah(:,:) = self%ahi(:,:) + if (config%lextra_force) call self%user_getacch(cb, config, t) + if (config%lgr) call self%gr_getacch(cb, config) + end associate + return + end subroutine helio_getacch_tp + + module subroutine helio_getacch_int_pl(self, t) + !! author: David A. Minton + !! + !! Compute direct cross term heliocentric accelerations of massive bodiese + !! + !! Adapted from David E. Kaufmann's Swifter routine helio_getacch_int.f90 + !! Adapted from Hal Levison's Swift routine getacch_ah3.f + implicit none + ! Arguments + class(helio_pl), intent(inout) :: self !! Helio massive body particle data structure + real(DP), intent(in) :: t !! Current time + ! Internals + integer(I4B) :: i, j + real(DP) :: rji2, irij3, faci, facj + real(DP), dimension(NDIM) :: dx + + associate(npl => self%nbody) + do i = 2, npl - 1 + do j = i + 1, npl + dx(:) = self%xh(:,j) - self%xh(:,i) + rji2 = dot_product(dx(:), dx(:)) + irij3 = 1.0_DP / (rji2 * sqrt(rji2)) + faci = self%Gmass(i) * irij3 + facj = self%Gmass(j) * irij3 + self%ahi(:,i) = self%ahi(:,i) + facj * dx(:) + self%ahi(:,i) = self%ahi(:,j) - faci * dx(:) + end do end do - end do - end associate - - return - end subroutine helio_getacch_int_pl + end associate + + return + end subroutine helio_getacch_int_pl - module subroutine helio_getacch_int_tp(self, pl, t, xh) - !! author: David A. Minton - !! - !! Compute direct cross term heliocentric accelerations of test particles - !! - !! Adapted from David E. Kaufmann's Swifter routine helio_getacch_int_tp.f90 - !! Adapted from Hal Levison's Swift routine getacch_ah3_tp.f - use swiftest - implicit none - ! Arguments - class(helio_tp), intent(inout) :: self !! Helio test particle data structure - class(whm_pl), intent(inout) :: pl !! Helio massive body particle data structure - real(DP), intent(in) :: t !! Current time - real(DP), dimension(:,:), intent(in) :: xh !! Heliocentric positions of planets - ! Internals - integer(I4B) :: i, j - real(DP) :: r2, fac - real(DP), dimension(NDIM) :: dx - - associate(ntp => self%nbody, npl => pl%nbody) - do i = 1, ntp - if (self%status(i) == ACTIVE) then - do j = 2, npl - dx(:) = self%xh(:,i) - xh(:,j) - r2 = dot_product(dx(:), dx(:)) - fac = pl%Gmass(j) / (r2 * sqrt(r2)) - self%ahi(:,i) = self%ahi(:,i) - fac * dx(:) + module subroutine helio_getacch_int_tp(self, pl, t, xh) + !! author: David A. Minton + !! + !! Compute direct cross term heliocentric accelerations of test particles + !! + !! Adapted from David E. Kaufmann's Swifter routine helio_getacch_int_tp.f90 + !! Adapted from Hal Levison's Swift routine getacch_ah3_tp.f + implicit none + ! Arguments + class(helio_tp), intent(inout) :: self !! Helio test particle data structure + class(whm_pl), intent(inout) :: pl !! Helio massive body particle data structure + real(DP), intent(in) :: t !! Current time + real(DP), dimension(:,:), intent(in) :: xh !! Heliocentric positions of planets + ! Internals + integer(I4B) :: i, j + real(DP) :: r2, fac + real(DP), dimension(NDIM) :: dx + + associate(ntp => self%nbody, npl => pl%nbody) + do i = 1, ntp + if (self%status(i) == ACTIVE) then + do j = 2, npl + dx(:) = self%xh(:,i) - xh(:,j) + r2 = dot_product(dx(:), dx(:)) + fac = pl%Gmass(j) / (r2 * sqrt(r2)) + self%ahi(:,i) = self%ahi(:,i) - fac * dx(:) + end do + end if end do - end if - end do - end associate - return - end subroutine helio_getacch_int_tp + end associate + return + end subroutine helio_getacch_int_tp end submodule s_helio_getacch diff --git a/src/helio/helio_setup.f90 b/src/helio/helio_setup.f90 index 7a73c8145..b31b43b26 100644 --- a/src/helio/helio_setup.f90 +++ b/src/helio/helio_setup.f90 @@ -1,4 +1,5 @@ submodule(helio_classes) s_helio_setup + use swiftest contains module procedure helio_setup_pl !! author: David A. Minton & Carlisle A. Wishard @@ -6,7 +7,6 @@ !! Allocate Helio planet structure !! !! Equivalent in functionality to David E. Kaufmann's Swifter routine helio_setup.f90 - use swiftest implicit none !> Call allocation method for great-grandparent class (we don't need Jacobi variables from WHM/RMVS) @@ -24,7 +24,6 @@ !! Allocate Helio test particle structure !! !! Equivalent in functionality to David E. Kaufmann's Swifter routine helio_setup.f90 - use swiftest implicit none !> Call allocation method for great-grandparent class diff --git a/src/helio/helio_step.f90 b/src/helio/helio_step.f90 index 2a0099bed..1f8a6d96a 100644 --- a/src/helio/helio_step.f90 +++ b/src/helio/helio_step.f90 @@ -1,4 +1,5 @@ submodule(helio_classes) s_helio_step + use swiftest contains module subroutine helio_step_system(cb, pl, tp, config) !! author: David A. Minton @@ -7,7 +8,6 @@ module subroutine helio_step_system(cb, pl, tp, config) !! !! Adapted from Hal Levison's Swift routine step_kdk.f !! Adapted from David E. Kaufmann's Swifter routine helio_step.f90 - use swiftest implicit none ! Arguments class(helio_cb), intent(inout) :: cb !! Helio central body object @@ -33,7 +33,6 @@ module subroutine helio_step_pl(self, cb, config, t, dt) !! !! Adapted from David E. Kaufmann's Swifter helio_step_pl.f90 !! Adapted from Hal Levison's Swift routine helio_step_pl.f - use swiftest implicit none ! Arguments class(helio_pl), intent(inout) :: self !! WHM massive body particle data structure @@ -73,7 +72,6 @@ module subroutine helio_step_tp(self, cb, pl, config, t, dt) !! !! Adapted from David E. Kaufmann's Swifter routine helio_step_tp.f90 !! Adapted from Hal Levison's Swift routine helio_step_tp.f - use swiftest implicit none ! Arguments class(helio_tp), intent(inout) :: self !! Helio test particle data structure diff --git a/src/io/io.f90 b/src/io/io.f90 index 5c1ab239a..e0aa90e47 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -1,5 +1,5 @@ submodule (swiftest_classes) s_io -use swiftest + use swiftest contains module subroutine io_config_reader(self, unit, iotype, v_list, iostat, iomsg) !! author: The Purdue Swiftest Team - David A. Minton, Carlisle A. Wishard, Jennifer L.L. Pouplin, and Jacob R. Elliott diff --git a/src/modules/rmvs_classes.f90 b/src/modules/rmvs_classes.f90 index 81512b0d7..29fddbd3d 100644 --- a/src/modules/rmvs_classes.f90 +++ b/src/modules/rmvs_classes.f90 @@ -105,13 +105,11 @@ module rmvs_classes end type rmvs_nbody_system interface - module subroutine rmvs_discard_pl_tp(self, cb, pl, config, t, dt) + module subroutine rmvs_discard_pl_tp(self, pl, t, dt) use swiftest_classes, only : swiftest_cb, swiftest_pl, swiftest_configuration implicit none class(rmvs_tp), intent(inout) :: self - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object class(swiftest_pl), intent(inout) :: pl !! WHM massive body particle data structure. - class(swiftest_configuration), intent(in) :: config !! configuration parameters real(DP), intent(in) :: t !! Current simulation time real(DP), intent(in) :: dt !! Stepsize end subroutine rmvs_discard_pl_tp diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index 54af47849..c4363b2e9 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -520,22 +520,6 @@ module subroutine io_read_initialize_system(self, config) class(swiftest_configuration), intent(inout) :: config !! Input collection of configuration parameters end subroutine io_read_initialize_system - module function io_read_line_swifter(iu, id, d1, d2, d3, d4, d5, d6, out_type, MASS, RADIUS) - implicit none - integer(I4B) :: io_read_line_swifter - integer(I4B), intent(in) :: iu !! Unit number associated with input binary file - integer(I4B), intent(out) :: id !! Planet or test particle identifier - real(DP), intent(out) :: d1 !! First quantity (semimajor axis (pericentric distance for a parabola) or heliocentric x ) - real(DP), intent(out) :: d2 !! Second quantity (eccentricity or heliocentric y ) - real(DP), intent(out) :: d3 !! Third quantity (inclination or heliocentric z ) - real(DP), intent(out) :: d4 !! Fourth quantity (longitude of the ascending node or heliocentric vx) - real(DP), intent(out) :: d5 !! Fifth quantity (argument of pericenter or heliocentric vy) - real(DP), intent(out) :: d6 !! Sixth quantity (mean anomaly or heliocentric vz) - real(DP), intent(out), optional :: MASS !! Optional mass (omitted for massless test particle) - real(DP), intent(out), optional :: RADIUS !! Optional radius (omitted for massless test particle) - character(*), intent(in) :: out_type !! Format of input binary file - end function io_read_line_swifter - module subroutine io_write_discard(self, config, discards) implicit none class(swiftest_nbody_system), intent(inout) :: self !! Swiftest system object diff --git a/src/obl/obl.f90 b/src/obl/obl.f90 index 6c6bd03ef..1a6b77d09 100644 --- a/src/obl/obl.f90 +++ b/src/obl/obl.f90 @@ -1,4 +1,5 @@ submodule (swiftest_classes) s_obl + use swiftest contains module procedure obl_acc_body !! author: David A. Minton @@ -8,7 +9,6 @@ !! Adapted from David E. Kaufmann's Swifter routine: obl_acc.f90 and obl_acc_tp.f90 !! Adapted from Hal Levison's Swift routine obl_acc.f and obl_acc_tp.f - use swiftest implicit none integer(I4B) :: i real(DP) :: r2, irh, rinv2, t0, t1, t2, t3, fac1, fac2 diff --git a/src/operators/operator_cross.f90 b/src/operators/operator_cross.f90 index 57a500e33..d6fd82944 100644 --- a/src/operators/operator_cross.f90 +++ b/src/operators/operator_cross.f90 @@ -1,10 +1,11 @@ submodule(swiftest_operators) operator_cross_implementation + use swiftest !! author: David A. Minton !! !! Contains implementations for the .cross. operator for all defined integer and real types !! Single vector implementations: C(1:3) = A(1:3) .cross. B(1:3) !! Vector list implementations: C(1:3, :) = A(1:3, :) .cross. B(1:3, :) - contains +contains module procedure operator_cross_sp implicit none diff --git a/src/orbel/orbel_el2xv.f90 b/src/orbel/orbel_el2xv.f90 index ec4e9f5d5..3f1a81fe9 100644 --- a/src/orbel/orbel_el2xv.f90 +++ b/src/orbel/orbel_el2xv.f90 @@ -1,11 +1,10 @@ submodule (swiftest_classes) s_orbel_el2xv + use swiftest contains - module procedure orbel_el2xv_vec !! author: David A. Minton !! !! A wrapper method that converts all of the cartesian position and velocity vectors of a Swiftest body object to orbital elements. - use swiftest implicit none integer(I4B) :: i @@ -33,7 +32,6 @@ pure subroutine orbel_el2xv(mu, a, ie, inc, capom, omega, capm, x, v) !! DATE WRITTEN: May 11, 1992. !! REVISIONS: May 26 - now use better Kepler solver for ellipses !! and hyperbolae called EHYBRID.F and FHYBRID.F - use swiftest implicit none real(DP), intent(in) :: mu real(DP), intent(in) :: a, ie, inc, capom, omega, capm diff --git a/src/orbel/orbel_scget.f90 b/src/orbel/orbel_scget.f90 index 84c92ed81..0cdb67c72 100644 --- a/src/orbel/orbel_scget.f90 +++ b/src/orbel/orbel_scget.f90 @@ -1,4 +1,5 @@ submodule (swiftest_classes) s_orbel_scget + use swiftest contains module procedure orbel_scget !! author: David A. Minton @@ -8,7 +9,6 @@ !! !! Adapted from David E. Kaufmann's Swifter routine: orbel_scget.f90 !! Adapted from Hal Levison's Swift routine orbel_scget.f - use swiftest implicit none integer(I4B) :: nper real(DP) :: x diff --git a/src/orbel/orbel_xv2aeq.f90 b/src/orbel/orbel_xv2aeq.f90 index ed4524903..8338d6559 100644 --- a/src/orbel/orbel_xv2aeq.f90 +++ b/src/orbel/orbel_xv2aeq.f90 @@ -1,4 +1,5 @@ submodule (swiftest_classes) s_orbel_xv2aeq + use swiftest contains module procedure orbel_xv2aeq !! author: David A. Minton @@ -7,7 +8,6 @@ !! !! Adapted from David E. Kaufmann's Swifter routine: orbel_xv2aeq.f90 !! Adapted from Luke Dones' Swift routine orbel_xv2aeq.f - use swiftest implicit none integer(I4B) :: iorbit_type real(DP) :: r, v2, h2, energy, fac diff --git a/src/orbel/orbel_xv2aqt.f90 b/src/orbel/orbel_xv2aqt.f90 index ba06cf6a8..3c8bf3f3e 100644 --- a/src/orbel/orbel_xv2aqt.f90 +++ b/src/orbel/orbel_xv2aqt.f90 @@ -1,4 +1,5 @@ submodule (swiftest_classes) s_orbel_xv2aqt + use swiftest contains module procedure orbel_xv2aqt ! (mu, px, py, pz, vx, vy, vz, a, q, capm, tperi !! author: David A. Minton @@ -9,7 +10,6 @@ !! tperi < 0 means nearest pericenter passage is in the past !! !! Adapted from David E. Kaufmann's Swifter routine: orbel_xv2aqt.f90 - use swiftest implicit none integer(I4B) :: iorbit_type real(DP) :: r, v2, h2, rdotv, energy, fac, w, face, cape, e, tmpf, capf, mm diff --git a/src/orbel/orbel_xv2el.f90 b/src/orbel/orbel_xv2el.f90 index ed1b54251..434925c7d 100644 --- a/src/orbel/orbel_xv2el.f90 +++ b/src/orbel/orbel_xv2el.f90 @@ -1,11 +1,11 @@ submodule (swiftest_classes) s_orbel_xv2el + use swiftest contains module procedure orbel_xv2el_vec !! author: David A. Minton !! !! A wrapper method that converts all of the cartesian position and velocity vectors of a Swiftest body object to orbital elements. - use swiftest implicit none integer(I4B) :: i @@ -33,7 +33,6 @@ pure subroutine orbel_xv2el(mu, x, v, a, e, inc, capom, omega, capm) !! !! Adapted from David E. Kaufmann's Swifter routine: orbel_xv2el.f90 !! Adapted from Martin Duncan's Swift routine orbel_xv2el.f - use swiftest implicit none real(DP), intent(in) :: mu real(DP), dimension(:), intent(in) :: x, v diff --git a/src/rmvs/rmvs_discard.f90 b/src/rmvs/rmvs_discard.f90 index 26ef19c68..4f75cd645 100644 --- a/src/rmvs/rmvs_discard.f90 +++ b/src/rmvs/rmvs_discard.f90 @@ -1,7 +1,7 @@ submodule(rmvs_classes) s_rmvs_discard + use swiftest contains - - module subroutine rmvs_discard_pl_tp(self, cb, pl, config, t, dt) + module subroutine rmvs_discard_pl_tp(self, pl, t, dt) !! author: David A. Minton !! !! Check to see if test particles should be discarded based on pericenter passage distances with respect to @@ -9,13 +9,10 @@ module subroutine rmvs_discard_pl_tp(self, cb, pl, config, t, dt) !! !! Adapted from Hal Levison's Swift routine discard_pl.f !! Adapted from Hal Levison's Swift routine rmvs_discard_pl.f90 - use swiftest implicit none ! Arguments class(rmvs_tp), intent(inout) :: self - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object class(swiftest_pl), intent(inout) :: pl !! WHM massive body particle data structure. - class(swiftest_configuration), intent(in) :: config !! configuration parameters real(DP), intent(in) :: t !! Current simulation time real(DP), intent(in) :: dt !! Stepsize ! Internals @@ -33,7 +30,7 @@ module subroutine rmvs_discard_pl_tp(self, cb, pl, config, t, dt) end if end associate end do - call discard_pl_tp(tp, cb, pl, config, t, dt) + call discard_pl_tp(tp, pl, t, dt) end associate end subroutine rmvs_discard_pl_tp diff --git a/src/rmvs/rmvs_encounter_check.f90 b/src/rmvs/rmvs_encounter_check.f90 index 8da39833e..62d2b8281 100644 --- a/src/rmvs/rmvs_encounter_check.f90 +++ b/src/rmvs/rmvs_encounter_check.f90 @@ -1,4 +1,5 @@ submodule (rmvs_classes) s_rmvs_chk + use swiftest contains module function rmvs_encounter_check_tp(self, cb, pl, dt, rts) result(lencounter) !! author: David A. Minton @@ -7,7 +8,6 @@ module function rmvs_encounter_check_tp(self, cb, pl, dt, rts) result(lencounter !! !! Adapted from David E. Kaufmann's Swifter routine: rmvs_chk.f90 !! Adapted from Hal Levison's Swift routine rmvs3_chk.f - use swiftest implicit none ! Arguments class(rmvs_tp), intent(inout) :: self !! RMVS test particle object @@ -79,7 +79,6 @@ module function rmvs_chk_ind(xr, vr, dt, r2crit) result(lflag) !! !! Adapted from David E. Kaufmann's Swifter routine: rmvs_chk_ind.f90 !! Adapted from Hal Levison's Swift routine rmvs_chk_ind.f - use swiftest implicit none ! Arguments real(DP), intent(in) :: dt, r2crit diff --git a/src/rmvs/rmvs_getacch.f90 b/src/rmvs/rmvs_getacch.f90 index af285dd32..7005de895 100644 --- a/src/rmvs/rmvs_getacch.f90 +++ b/src/rmvs/rmvs_getacch.f90 @@ -1,6 +1,6 @@ submodule(rmvs_classes) s_rmvs_getacch + use swiftest contains - module subroutine rmvs_getacch_in_tp(self, cb, pl, config, t, xh) !! author: David A. Minton !! @@ -8,7 +8,6 @@ module subroutine rmvs_getacch_in_tp(self, cb, pl, config, t, xh) !! !! Performs a similar task as David E. Kaufmann's Swifter routine rmvs_getacch_tp.f90, but !! uses object polymorphism, and so is not directly adapted. - use swiftest implicit none ! Arguments class(rmvs_tp), intent(inout) :: self !! RMVS test particle data structure diff --git a/src/rmvs/rmvs_interp.f90 b/src/rmvs/rmvs_interp.f90 index a0cc8586b..a5e239ceb 100644 --- a/src/rmvs/rmvs_interp.f90 +++ b/src/rmvs/rmvs_interp.f90 @@ -1,4 +1,5 @@ submodule (rmvs_classes) s_rmvs_interp + use swiftest contains module subroutine rmvs_interp_in(self, cb, dt) !! author: David A. Minton @@ -8,7 +9,6 @@ module subroutine rmvs_interp_in(self, cb, dt) !! Adapted from David E. Kaufmann's Swifter routine rmvs_interp_in.f90 !! !! Adapted from Hal Levison's Swift routine rmvs3_interp.f - use swiftest implicit none ! Arguments class(rmvs_pl), intent(inout) :: self !! RMVS test particle object @@ -69,7 +69,6 @@ module subroutine rmvs_interp_out(self, cb, dt) !! Adapted from David E. Kaufmann's Swifter routine rmvs_interp_out.f90 !! !! Adapted from Hal Levison's Swift routine rmvs3_interp.f - use swiftest implicit none ! Arguments class(rmvs_pl), intent(inout) :: self !! RMVS test particle object diff --git a/src/rmvs/rmvs_obl.f90 b/src/rmvs/rmvs_obl.f90 index 0e3fa7561..80871fc72 100644 --- a/src/rmvs/rmvs_obl.f90 +++ b/src/rmvs/rmvs_obl.f90 @@ -1,13 +1,11 @@ submodule(rmvs_classes) s_rmvs_obl + use swiftest contains - module subroutine rmvs_obl_acc_in(self, cb) !! author: David A. Minton !! !! Compute the oblateness acceleration in the inner encounter region with planets !! - - use swiftest implicit none ! Arguments class(rmvs_pl), intent(inout) :: self !! RMVS massive body object diff --git a/src/rmvs/rmvs_peri.f90 b/src/rmvs/rmvs_peri.f90 index b2a24a314..7265c4693 100644 --- a/src/rmvs/rmvs_peri.f90 +++ b/src/rmvs/rmvs_peri.f90 @@ -1,4 +1,5 @@ submodule(rmvs_classes) s_rmvs_peri + use swiftest contains module subroutine rmvs_peri_tp(self, cb, pl, t, dt, lfirst, index, nenc, ipleP, config) !! author: David A. Minton @@ -7,7 +8,6 @@ module subroutine rmvs_peri_tp(self, cb, pl, t, dt, lfirst, index, nenc, ipleP, !! !! Adapted from Hal Levison's Swift routine Adapted from Hal Levison's Swift routine util_peri.f !! Adapted from David E. Kaufmann's Swifter routine rmvs_peri.f90 - use swiftest implicit none ! Arguments class(rmvs_tp), intent(inout) :: self !! RMVS test particle object diff --git a/src/rmvs/rmvs_setup.f90 b/src/rmvs/rmvs_setup.f90 index d651b9413..48382d866 100644 --- a/src/rmvs/rmvs_setup.f90 +++ b/src/rmvs/rmvs_setup.f90 @@ -1,4 +1,5 @@ submodule(rmvs_classes) s_rmvs_setup + use swiftest contains module subroutine rmvs_setup_pl(self,n) !! author: David A. Minton @@ -6,7 +7,6 @@ module subroutine rmvs_setup_pl(self,n) !! Allocate RMVS test particle structure !! !! Equivalent in functionality to David E. Kaufmann's Swifter routine rmvs_setup.f90 - use swiftest implicit none ! Arguments class(rmvs_pl), intent(inout) :: self !! RMVS test particle object @@ -41,7 +41,6 @@ module subroutine rmvs_setup_tp(self,n) !! Allocate WHM test particle structure !! !! Equivalent in functionality to David E. Kaufmann's Swifter routine whm_setup.f90 - use swiftest implicit none ! Arguments class(rmvs_tp), intent(inout) :: self !! RMVS test particle object @@ -69,7 +68,6 @@ module subroutine rmvs_setup_system(self, config) !! !! Wrapper method to initialize a basic Swiftest nbody system from files !! - use swiftest implicit none ! Arguments class(rmvs_nbody_system), intent(inout) :: self !! RMVS system object @@ -84,7 +82,6 @@ module subroutine rmvs_setup_set_beg_end(self, xbeg, xend, vbeg) !! author: David A. Minton !! !! Sets one or more of the values of xbeg, xend, and vbeg - use swiftest implicit none ! Arguments class(rmvs_tp), intent(inout) :: self !! RMVS test particle object diff --git a/src/rmvs/rmvs_spill_and_fill.f90 b/src/rmvs/rmvs_spill_and_fill.f90 index d29cbdbc8..6369375e0 100644 --- a/src/rmvs/rmvs_spill_and_fill.f90 +++ b/src/rmvs/rmvs_spill_and_fill.f90 @@ -1,111 +1,109 @@ submodule(rmvs_classes) s_rmvs_spill_and_fill -contains -module subroutine rmvs_spill_pl(self, discards, lspill_list) - !! author: David A. Minton - !! - !! Move spilled (discarded) RMVS test particle structure from active list to discard list - !! - !! Adapted from David E. Kaufmann's Swifter routine discard_discard_spill.f90 use swiftest - implicit none - ! Arguments - class(rmvs_pl), intent(inout) :: self !! Swiftest massive body body object - class(swiftest_body), intent(inout) :: discards !! Discarded object - logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards - ! Internals - integer(I4B) :: i - - associate(keeps => self) - select type(discards) - class is (rmvs_pl) - - discards%nenc(:) = pack(keeps%nenc(:), lspill_list(:)) - discards%tpenc1P(:) = pack(keeps%tpenc1P(:), lspill_list(:)) - if (count(.not.lspill_list(:)) > 0) then - keeps%nenc(:) = pack(keeps%nenc(:), .not. lspill_list(:)) - keeps%tpenc1P(:) = pack(keeps%tpenc1P(:), .not. lspill_list(:)) - end if - call whm_spill_pl(keeps, discards, lspill_list) - class default - write(*,*) 'Error! spill method called for incompatible return type on rmvs_pl' - end select - end associate - - return - - end subroutine rmvs_spill_pl - - module subroutine rmvs_fill_pl(self, inserts, lfill_list) +contains + module subroutine rmvs_spill_pl(self, discards, lspill_list) !! author: David A. Minton !! - !! Insert new RMVS massive body structure into an old one. - !! This is the inverse of a fill operation. + !! Move spilled (discarded) RMVS test particle structure from active list to discard list !! - use swiftest + !! Adapted from David E. Kaufmann's Swifter routine discard_discard_spill.f90 implicit none ! Arguments - class(rmvs_pl), intent(inout) :: self !! RMVS massive body object - class(swiftest_body), intent(inout) :: inserts !! Inserted object - logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps + class(rmvs_pl), intent(inout) :: self !! Swiftest massive body body object + class(swiftest_body), intent(inout) :: discards !! Discarded object + logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards ! Internals - integer(I4B) :: i - + integer(I4B) :: i + associate(keeps => self) - select type(inserts) + select type(discards) class is (rmvs_pl) - - keeps%nenc(:) = unpack(keeps%nenc(:), .not.lfill_list(:), keeps%nenc(:)) - keeps%nenc(:) = unpack(inserts%nenc(:), lfill_list(:), keeps%nenc(:)) - - keeps%tpenc1P(:) = unpack(keeps%tpenc1P(:), .not.lfill_list(:), keeps%tpenc1P(:)) - keeps%tpenc1P(:) = unpack(inserts%tpenc1P(:), lfill_list(:), keeps%tpenc1P(:)) - - call whm_fill_pl(keeps, inserts, lfill_list) + + discards%nenc(:) = pack(keeps%nenc(:), lspill_list(:)) + discards%tpenc1P(:) = pack(keeps%tpenc1P(:), lspill_list(:)) + if (count(.not.lspill_list(:)) > 0) then + keeps%nenc(:) = pack(keeps%nenc(:), .not. lspill_list(:)) + keeps%tpenc1P(:) = pack(keeps%tpenc1P(:), .not. lspill_list(:)) + end if + call whm_spill_pl(keeps, discards, lspill_list) class default write(*,*) 'Error! spill method called for incompatible return type on rmvs_pl' end select end associate - + return - + + end subroutine rmvs_spill_pl + + module subroutine rmvs_fill_pl(self, inserts, lfill_list) + !! author: David A. Minton + !! + !! Insert new RMVS massive body structure into an old one. + !! This is the inverse of a fill operation. + !! + implicit none + ! Arguments + class(rmvs_pl), intent(inout) :: self !! RMVS massive body object + class(swiftest_body), intent(inout) :: inserts !! Inserted object + logical, dimension(:), intent(in) :: lfill_list !! Logical array of bodies to merge into the keeps + ! Internals + integer(I4B) :: i + + associate(keeps => self) + select type(inserts) + class is (rmvs_pl) + + keeps%nenc(:) = unpack(keeps%nenc(:), .not.lfill_list(:), keeps%nenc(:)) + keeps%nenc(:) = unpack(inserts%nenc(:), lfill_list(:), keeps%nenc(:)) + + keeps%tpenc1P(:) = unpack(keeps%tpenc1P(:), .not.lfill_list(:), keeps%tpenc1P(:)) + keeps%tpenc1P(:) = unpack(inserts%tpenc1P(:), lfill_list(:), keeps%tpenc1P(:)) + + call whm_fill_pl(keeps, inserts, lfill_list) + class default + write(*,*) 'Error! spill method called for incompatible return type on rmvs_pl' + end select + end associate + + return + end subroutine rmvs_fill_pl module subroutine rmvs_spill_tp(self, discards, lspill_list) - !! author: David A. Minton - !! - !! Move spilled (discarded) RMVS test particle structure from active list to discard list - !! - !! Adapted from David E. Kaufmann's Swifter routine whm_discard_spill.f90 - use swiftest - implicit none - ! Arguments - class(rmvs_tp), intent(inout) :: self !! RMVS test particle object - class(swiftest_body), intent(inout) :: discards !! Discarded object - logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards - ! Internals - integer(I4B) :: i + !! author: David A. Minton + !! + !! Move spilled (discarded) RMVS test particle structure from active list to discard list + !! + !! Adapted from David E. Kaufmann's Swifter routine whm_discard_spill.f90 + implicit none + ! Arguments + class(rmvs_tp), intent(inout) :: self !! RMVS test particle object + class(swiftest_body), intent(inout) :: discards !! Discarded object + logical, dimension(:), intent(in) :: lspill_list !! Logical array of bodies to spill into the discards + ! Internals + integer(I4B) :: i - associate(keeps => self) - select type(discards) - class is (rmvs_tp) - discards%lperi(:) = pack(keeps%lperi(:), lspill_list(:)) - discards%plperP(:) = pack(keeps%plperP(:), lspill_list(:)) - discards%plencP(:) = pack(keeps%plencP(:), lspill_list(:)) - discards%tpencP(:) = pack(keeps%tpencP(:), lspill_list(:)) - if (count(.not.lspill_list(:)) > 0) then - keeps%lperi(:) = pack(keeps%lperi(:), .not. lspill_list(:)) - keeps%plperP(:) = pack(keeps%plperP(:), .not. lspill_list(:)) - keeps%plencP(:) = pack(keeps%plencP(:), .not. lspill_list(:)) - keeps%tpencP(:) = pack(keeps%tpencP(:), .not. lspill_list(:)) - end if + associate(keeps => self) + select type(discards) + class is (rmvs_tp) + discards%lperi(:) = pack(keeps%lperi(:), lspill_list(:)) + discards%plperP(:) = pack(keeps%plperP(:), lspill_list(:)) + discards%plencP(:) = pack(keeps%plencP(:), lspill_list(:)) + discards%tpencP(:) = pack(keeps%tpencP(:), lspill_list(:)) + if (count(.not.lspill_list(:)) > 0) then + keeps%lperi(:) = pack(keeps%lperi(:), .not. lspill_list(:)) + keeps%plperP(:) = pack(keeps%plperP(:), .not. lspill_list(:)) + keeps%plencP(:) = pack(keeps%plencP(:), .not. lspill_list(:)) + keeps%tpencP(:) = pack(keeps%tpencP(:), .not. lspill_list(:)) + end if - call util_spill_tp(keeps, discards, lspill_list) - class default - write(*,*) 'Error! spill method called for incompatible return type on rmvs_tp' - end select - end associate + call util_spill_tp(keeps, discards, lspill_list) + class default + write(*,*) 'Error! spill method called for incompatible return type on rmvs_tp' + end select + end associate - return + return end subroutine rmvs_spill_tp @@ -115,7 +113,6 @@ module subroutine rmvs_fill_tp(self, inserts, lfill_list) !! Insert new RMVS test particle structure into an old one. !! This is the inverse of a fill operation. !! - use swiftest implicit none ! Arguments class(rmvs_tp), intent(inout) :: self !! RMVS massive body object @@ -145,7 +142,7 @@ module subroutine rmvs_fill_tp(self, inserts, lfill_list) end associate return - - end subroutine rmvs_fill_tp + + end subroutine rmvs_fill_tp end submodule s_rmvs_spill_and_fill diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index f8165de31..f65c9337e 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -1,4 +1,5 @@ submodule(rmvs_classes) s_rmvs_step + use swiftest contains module subroutine rmvs_step_system(cb, pl, tp, config) !! author: David A. Minton @@ -7,7 +8,6 @@ module subroutine rmvs_step_system(cb, pl, tp, config) !! !! Adapted from Hal Levison's Swift routine rmvs3_step.f !! Adapted from David E. Kaufmann's Swifter routine rmvs_step.f90 - use swiftest implicit none ! Arguments class(rmvs_cb), intent(inout) :: cb !! RMVS central body object @@ -62,7 +62,6 @@ module subroutine rmvs_step_out(self, cb, tp, dt, config) !! !! Adapted from Hal Levison's Swift routine rmvs3_step_out.f !! Adapted from David E. Kaufmann's Swifter routine rmvs_step_out.f90 - use swiftest implicit none ! Arguments class(rmvs_pl), intent(inout) :: self !! RMVS massive body object @@ -73,7 +72,6 @@ module subroutine rmvs_step_out(self, cb, tp, dt, config) ! Internals integer(I4B) :: i, j, k, nenc, itpp real(DP) :: dto, time - type(rmvs_pl) :: rmvs_plep ! executable code associate(pl => self, npl => self%nbody, ntp => tp%nbody, t => config%t, xht => tp%xh, vht => tp%vh, & @@ -101,251 +99,246 @@ module subroutine rmvs_step_out(self, cb, tp, dt, config) end associate return - end subroutine rmvs_step_out + end subroutine rmvs_step_out - module subroutine rmvs_step_out2(cb, pl, tp, t, dt, index, config) - !! author: David A. Minton - !! - !! Step active test particles ahead in the outer encounter region, setting up and calling the inner region - !! integration if necessary - !! - !! Adapted from Hal Levison's Swift routine rmvs3_step_out.f - !! Adapted from David E. Kaufmann's Swifter routine rmvs_step_out.f90 - use swiftest - implicit none - ! Arguments - class(rmvs_cb), intent(inout) :: cb !! RMVS central body object - class(rmvs_pl), intent(inout) :: pl !! RMVS massive body object - class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object - real(DP), intent(in) :: t !! Simulation time - real(DP), intent(in) :: dt !! Step size - integer(I4B), intent(in) :: index !! outer substep number within current set - class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters - ! Internals - logical :: lfirsttp - integer(I4B) :: i - real(DP) :: rts - logical :: lencounter + module subroutine rmvs_step_out2(cb, pl, tp, t, dt, index, config) + !! author: David A. Minton + !! + !! Step active test particles ahead in the outer encounter region, setting up and calling the inner region + !! integration if necessary + !! + !! Adapted from Hal Levison's Swift routine rmvs3_step_out.f + !! Adapted from David E. Kaufmann's Swifter routine rmvs_step_out.f90 + implicit none + ! Arguments + class(rmvs_cb), intent(inout) :: cb !! RMVS central body object + class(rmvs_pl), intent(inout) :: pl !! RMVS massive body object + class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object + real(DP), intent(in) :: t !! Simulation time + real(DP), intent(in) :: dt !! Step size + integer(I4B), intent(in) :: index !! outer substep number within current set + class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters + ! Internals + logical :: lfirsttp + real(DP) :: rts + logical :: lencounter - associate(xht => tp%xh, vht => tp%vh, xbeg => tp%xbeg, xend => tp%xend, aht => tp%ah,& - xout1 => pl%xout(:, :, index-1), vout1 => pl%vout(:,:,index-1), xout2 => pl%xout(:,:,index)) + associate(xht => tp%xh, vht => tp%vh, xbeg => tp%xbeg, xend => tp%xend, aht => tp%ah,& + xout1 => pl%xout(:, :, index-1), vout1 => pl%vout(:,:,index-1), xout2 => pl%xout(:,:,index)) - call tp%set_beg_end(xbeg = pl%xout(:, :, index-1), & - vbeg = pl%vout(:, :, index-1), & - xend = pl%xout(:, :, index)) - rts = RHPSCALE - lencounter = tp%encounter_check(cb, pl, dt, rts) - if (lencounter) then - ! Interpolate planets in inner encounter region - pl%xin(:,:,0) = tp%xbeg(:,:) - pl%vin(:,:,0) = tp%vbeg(:,:) - pl%xin(:,:,NTPHENC) = tp%xend(:, :) - pl%vin(:,:,NTPHENC) = pl%vout(:, :, index) - call pl%interp_in(cb, dt) - call pl%step_in(cb, tp, config, t, dt) - lfirsttp = tp%lfirst - tp%lfirst = .true. - call tp%step(cb, pl, config, t, dt) - tp%lfirst = lfirsttp - else - call tp%step(cb, pl, config, t, dt) - end if - end associate - return - - end subroutine rmvs_step_out2 + call tp%set_beg_end(xbeg = pl%xout(:, :, index-1), & + vbeg = pl%vout(:, :, index-1), & + xend = pl%xout(:, :, index)) + rts = RHPSCALE + lencounter = tp%encounter_check(cb, pl, dt, rts) + if (lencounter) then + ! Interpolate planets in inner encounter region + pl%xin(:,:,0) = tp%xbeg(:,:) + pl%vin(:,:,0) = tp%vbeg(:,:) + pl%xin(:,:,NTPHENC) = tp%xend(:, :) + pl%vin(:,:,NTPHENC) = pl%vout(:, :, index) + call pl%interp_in(cb, dt) + call pl%step_in(cb, tp, config, t, dt) + lfirsttp = tp%lfirst + tp%lfirst = .true. + call tp%step(cb, pl, config, t, dt) + tp%lfirst = lfirsttp + else + call tp%step(cb, pl, config, t, dt) + end if + end associate + return - module subroutine rmvs_step_in_pl(self, cb, tp, config, t, dt) - !! author: David A. Minton - !! - !! Step active test particles ahead in the inner encounter region - !! - !! Adapted from Hal Levison's Swift routine rmvs3_step_in.f - !! Adapted from David E. Kaufmann's Swifter routine rmvs_step_in.f90 - use swiftest - implicit none - ! Arguments - class(rmvs_pl), intent(inout) :: self !! RMVS massive body object - class(rmvs_cb), intent(inout) :: cb !! RMVS central body object - class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object - class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters - real(DP), intent(in) :: t !! Current time - real(DP), intent(in) :: dt !! Step size - ! Internals - logical :: lfirsttp - integer(I4B) :: i, j, k, nenc, link - real(DP) :: mu, rhill, dti, time + end subroutine rmvs_step_out2 - dti = dt / NTPHENC - associate(pl => self, npl => self%nbody, xht => tp%xh, vht => tp%vh) - if (config%loblatecb) call pl%obl_acc_in(cb) - call pl%make_planetocentric(cb, tp, config) - do i = 1, npl - nenc = pl%nenc(i) - if (nenc > 0) then - ! There are inner encounters with this planet...switch to planetocentric coordinates to proceed - time = t - call pl%tpenc(i)%peri_pass(cb, pl, time, dti, .true., 0, nenc, i, config) - ! now step the encountering test particles fully through the inner encounter - lfirsttp = .true. - associate(index => pl%tpenc(i)%index, & - xpc => self%tpenc(i)%xh, vpc => self%tpenc(i)%vh, apc => self%tpenc(i)%ah) - pl%tpenc(i)%lfirst = .true. - do index = 1, NTPHENC ! Integrate over the encounter region, using the "substitute" planetocentric systems at each level - call pl%tpenc(i)%set_beg_end(xbeg = pl%plenc(i,index - 1)%xh, xend = pl%plenc(i, index)%xh) - call pl%tpenc(i)%step(pl%cbenc(i), pl%plenc(i, index), config, time, dti) - do j = 1, NDIM - pl%tpenc(i)%xheliocen(j, :) = pl%tpenc(i)%xh(j, :) + pl%xin(j, i, index) - end do - time = config%t + j * dti - call pl%tpenc(i)%peri_pass(cb, pl, time, dti, .false., index, nenc, i, config) + module subroutine rmvs_step_in_pl(self, cb, tp, config, t, dt) + !! author: David A. Minton + !! + !! Step active test particles ahead in the inner encounter region + !! + !! Adapted from Hal Levison's Swift routine rmvs3_step_in.f + !! Adapted from David E. Kaufmann's Swifter routine rmvs_step_in.f90 + implicit none + ! Arguments + class(rmvs_pl), intent(inout) :: self !! RMVS massive body object + class(rmvs_cb), intent(inout) :: cb !! RMVS central body object + class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object + class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters + real(DP), intent(in) :: t !! Current time + real(DP), intent(in) :: dt !! Step size + ! Internals + logical :: lfirsttp + integer(I4B) :: i, j, nenc + real(DP) :: dti, time + + dti = dt / NTPHENC + associate(pl => self, npl => self%nbody, xht => tp%xh, vht => tp%vh) + if (config%loblatecb) call pl%obl_acc_in(cb) + call pl%make_planetocentric(cb, tp, config) + do i = 1, npl + nenc = pl%nenc(i) + if (nenc > 0) then + ! There are inner encounters with this planet...switch to planetocentric coordinates to proceed + time = t + call pl%tpenc(i)%peri_pass(cb, pl, time, dti, .true., 0, nenc, i, config) + ! now step the encountering test particles fully through the inner encounter + lfirsttp = .true. + associate(index => pl%tpenc(i)%index, & + xpc => self%tpenc(i)%xh, vpc => self%tpenc(i)%vh, apc => self%tpenc(i)%ah) + pl%tpenc(i)%lfirst = .true. + do index = 1, NTPHENC ! Integrate over the encounter region, using the "substitute" planetocentric systems at each level + call pl%tpenc(i)%set_beg_end(xbeg = pl%plenc(i,index - 1)%xh, xend = pl%plenc(i, index)%xh) + call pl%tpenc(i)%step(pl%cbenc(i), pl%plenc(i, index), config, time, dti) + do j = 1, NDIM + pl%tpenc(i)%xheliocen(j, :) = pl%tpenc(i)%xh(j, :) + pl%xin(j, i, index) end do - where(pl%tpenc(i)%status(:) == ACTIVE) pl%tpenc(i)%status(:) = INACTIVE - end associate - end if - end do - call pl%end_planetocentric(cb,tp) + time = config%t + j * dti + call pl%tpenc(i)%peri_pass(cb, pl, time, dti, .false., index, nenc, i, config) + end do + where(pl%tpenc(i)%status(:) == ACTIVE) pl%tpenc(i)%status(:) = INACTIVE + end associate + end if + end do + call pl%end_planetocentric(cb,tp) + + end associate - end associate + return + end subroutine rmvs_step_in_pl - return - end subroutine rmvs_step_in_pl + module subroutine rmvs_step_make_planetocentric(self, cb, tp, config) + !! author: David A. Minton + !! + !! When encounters are detected, this method will call the interpolation methods for the planets and + !! creates a Swiftest test particle structure for each planet's encountering test particles to simplify the + !! planetocentric calculations. This subroutine is not based on an existing one from Swift and Swifter + !! + implicit none + ! Arguments + class(rmvs_pl), intent(inout) :: self !! RMVS test particle object + class(rmvs_cb), intent(inout) :: cb !! RMVS central body particle type + class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object + class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters + ! Internals + integer(I4B) :: i, j, k + type(rmvs_pl) :: cb_as_pl, tmp + logical, dimension(:), allocatable :: copyflag - module subroutine rmvs_step_make_planetocentric(self, cb, tp, config) - !! author: David A. Minton - !! - !! When encounters are detected, this method will call the interpolation methods for the planets and - !! creates a Swiftest test particle structure for each planet's encountering test particles to simplify the - !! planetocentric calculations. This subroutine is not based on an existing one from Swift and Swifter - !! - use swiftest - implicit none - ! Arguments - class(rmvs_pl), intent(inout) :: self !! RMVS test particle object - class(rmvs_cb), intent(inout) :: cb !! RMVS central body particle type - class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object - class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters - ! Internals - integer(I4B) :: i, j, k, link, nenc, ntp - type(rmvs_pl) :: cb_as_pl, tmp - logical, dimension(:), allocatable :: copyflag - - allocate(self%tpenc(self%nbody)) - allocate(self%plenc(self%nbody, 0:NTPHENC)) - allocate(self%cbenc(self%nbody)) - allocate(copyflag(self%nbody)) - associate(pl => self, npl => self%nbody, nenc => self%nenc, tpenc => self%tpenc, cbenc => self%cbenc, & - plenc => self%plenc) - ! Indicate that this is a planetocentric close encounter structor + allocate(self%tpenc(self%nbody)) + allocate(self%plenc(self%nbody, 0:NTPHENC)) + allocate(self%cbenc(self%nbody)) + allocate(copyflag(self%nbody)) + associate(pl => self, npl => self%nbody, nenc => self%nenc, tpenc => self%tpenc, cbenc => self%cbenc, & + plenc => self%plenc) + ! Indicate that this is a planetocentric close encounter structor - ! Save the original central body object so that it can be passed down as needed through the planetocentric structures - call cb_as_pl%setup(1) - call tmp%setup(1) - cb_as_pl%name(1) = 0 - cb_as_pl%Gmass(1) = cb%Gmass - cb_as_pl%mass(1) = cb%mass - cb_as_pl%radius(1) = cb%radius + ! Save the original central body object so that it can be passed down as needed through the planetocentric structures + call cb_as_pl%setup(1) + call tmp%setup(1) + cb_as_pl%name(1) = 0 + cb_as_pl%Gmass(1) = cb%Gmass + cb_as_pl%mass(1) = cb%mass + cb_as_pl%radius(1) = cb%radius - do i = 1, npl - if (nenc(i) > 0) then ! There are inner encounters with this planet - ! Make the planet a central body for the encounter - cbenc(i) = cb - cbenc(i)%Gmass = pl%Gmass(i) - cbenc(i)%mass = pl%mass(i) - cbenc(i)%radius = pl%radius(i) + do i = 1, npl + if (nenc(i) > 0) then ! There are inner encounters with this planet + ! Make the planet a central body for the encounter + cbenc(i) = cb + cbenc(i)%Gmass = pl%Gmass(i) + cbenc(i)%mass = pl%mass(i) + cbenc(i)%radius = pl%radius(i) - ! Create space for the heliocentric position values for those acceleration calculations that need them - allocate(tpenc(i)%xheliocen(NDIM, nenc(i))) + ! Create space for the heliocentric position values for those acceleration calculations that need them + allocate(tpenc(i)%xheliocen(NDIM, nenc(i))) - ! Save the index value of the planet corresponding to this encounter - tpenc(i)%ipleP = i - tpenc(i)%lplanetocentric = .true. - tpenc(i)%cb = cb - call pl%tpenc(i)%setup(nenc(i)) - tpenc(i)%status(:) = ACTIVE - call tp%spill(pl%tpenc(i), pl%encmask(:,i)) - ! Grab all the encountering test particles and convert them to a planetocentric frame - do j = 1, nenc(i) - tpenc(i)%xheliocen(:, j) = tpenc(i)%xh(:, j) - tpenc(i)%xh(:, j) = tpenc(i)%xh(:, j) - pl%xin(:, i, 0) - tpenc(i)%vh(:, j) = tpenc(i)%vh(:, j) - pl%vin(:, i, 0) - end do + ! Save the index value of the planet corresponding to this encounter + tpenc(i)%ipleP = i + tpenc(i)%lplanetocentric = .true. + tpenc(i)%cb = cb + call pl%tpenc(i)%setup(nenc(i)) + tpenc(i)%status(:) = ACTIVE + call tp%spill(pl%tpenc(i), pl%encmask(:,i)) + ! Grab all the encountering test particles and convert them to a planetocentric frame + do j = 1, nenc(i) + tpenc(i)%xheliocen(:, j) = tpenc(i)%xh(:, j) + tpenc(i)%xh(:, j) = tpenc(i)%xh(:, j) - pl%xin(:, i, 0) + tpenc(i)%vh(:, j) = tpenc(i)%vh(:, j) - pl%vin(:, i, 0) + end do - ! Make sure that the test particles get the planetocentric value of mu - call tpenc(i)%set_mu(pl%cbenc(i)) + ! Make sure that the test particles get the planetocentric value of mu + call tpenc(i)%set_mu(pl%cbenc(i)) - ! Save the heliocentric position values of the encountering planet - allocate(tpenc(i)%xh_pl(NDIM,0:NTPHENC)) - tpenc(i)%xh_pl(:,:) = pl%xin(:,i,0:NTPHENC) - ! Save the encountering planet's values of oblateness acceleration - if (config%loblatecb) then - allocate(tpenc(i)%aoblin_pl(NDIM,0:NTPHENC)) - tpenc(i)%aoblin_pl(:,:) = pl%aoblin(:,i,0:NTPHENC) - end if - ! Now create a planetocentric "planet" structure containing the *other* planets (plus the Sun) in it at each point along - ! the innner encounter trajectory of the planet - do k = 0, NTPHENC - call plenc(i, k)%setup(npl) - ! Copy all the basic planet parameters and positions - copyflag(:) = .true. - call plenc(i, k)%copy(pl,copyflag) - ! Give each planet a position and velocity vector that is planetocentric wrt planet i (the encountering planet) - do j = 1, npl - plenc(i, k)%xh(:, j) = pl%xin(:, j, k) - pl%xin(:, i, k) - plenc(i, k)%vh(:, j) = pl%vin(:, j, k) - pl%vin(:, i, k) - end do - cb_as_pl%xh(:, 1) = cb%xin(:) - pl%xin(:, i, k) - cb_as_pl%vh(:, 1) = cb%vin(:) - pl%vin(:, i, k) - ! Pull the encountering body out of the massive body list - copyflag(:) = .false. - copyflag(i) = .true. - call plenc(i, k)%spill(tmp, copyflag) - copyflag(:) = .false. - copyflag(1) = .true. ! Put the central body as planet 1 like in the original Swifter version - call plenc(i, k)%fill(cb_as_pl, copyflag) - end do + ! Save the heliocentric position values of the encountering planet + allocate(tpenc(i)%xh_pl(NDIM,0:NTPHENC)) + tpenc(i)%xh_pl(:,:) = pl%xin(:,i,0:NTPHENC) + ! Save the encountering planet's values of oblateness acceleration + if (config%loblatecb) then + allocate(tpenc(i)%aoblin_pl(NDIM,0:NTPHENC)) + tpenc(i)%aoblin_pl(:,:) = pl%aoblin(:,i,0:NTPHENC) end if - end do - end associate - end subroutine rmvs_step_make_planetocentric - - module subroutine rmvs_step_end_planetocentric(self, cb, tp) - !! author: David A. Minton - !! - !! Deallocates all of the encountering particle data structures for next time - !! - use swiftest - implicit none - ! Arguments - class(rmvs_pl), intent(inout) :: self !! RMVS test particle object - class(rmvs_cb), intent(inout) :: cb !! RMVS central body object - class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object - ! Internals - integer(I4B) :: i, j + ! Now create a planetocentric "planet" structure containing the *other* planets (plus the Sun) in it at each point along + ! the innner encounter trajectory of the planet + do k = 0, NTPHENC + call plenc(i, k)%setup(npl) + ! Copy all the basic planet parameters and positions + copyflag(:) = .true. + call plenc(i, k)%copy(pl,copyflag) + ! Give each planet a position and velocity vector that is planetocentric wrt planet i (the encountering planet) + do j = 1, npl + plenc(i, k)%xh(:, j) = pl%xin(:, j, k) - pl%xin(:, i, k) + plenc(i, k)%vh(:, j) = pl%vin(:, j, k) - pl%vin(:, i, k) + end do + cb_as_pl%xh(:, 1) = cb%xin(:) - pl%xin(:, i, k) + cb_as_pl%vh(:, 1) = cb%vin(:) - pl%vin(:, i, k) + ! Pull the encountering body out of the massive body list + copyflag(:) = .false. + copyflag(i) = .true. + call plenc(i, k)%spill(tmp, copyflag) + copyflag(:) = .false. + copyflag(1) = .true. ! Put the central body as planet 1 like in the original Swifter version + call plenc(i, k)%fill(cb_as_pl, copyflag) + end do + end if + end do + end associate + end subroutine rmvs_step_make_planetocentric - associate(pl => self, nenc => self%nenc, npl => self%nbody, name => tp%name, & - tpenc => self%tpenc, plenc => self%plenc, cbenc => self%cbenc, encmask => self%encmask) + module subroutine rmvs_step_end_planetocentric(self, cb, tp) + !! author: David A. Minton + !! + !! Deallocates all of the encountering particle data structures for next time + !! + implicit none + ! Arguments + class(rmvs_pl), intent(inout) :: self !! RMVS test particle object + class(rmvs_cb), intent(inout) :: cb !! RMVS central body object + class(rmvs_tp), intent(inout) :: tp !! RMVS test particle object + ! Internals + integer(I4B) :: i, j - do i = 1, npl - if (nenc(i) == 0) cycle - do j = 1, nenc(i) - ! Copy the results of the integration back over - tpenc(i)%xh(:, j) = tpenc(i)%xh(:, j) + pl%xin(:, i, NTPHENC) - tpenc(i)%vh(:, j) = tpenc(i)%vh(:, j) + pl%vin(:, i, NTPHENC) - ! Put back heliocentric value of mu - call tpenc(i)%set_mu(cb) - end do - ! Replace the old test particle with the new one - call tp%fill(tpenc(i), encmask(:,i)) + associate(pl => self, nenc => self%nenc, npl => self%nbody, name => tp%name, & + tpenc => self%tpenc, plenc => self%plenc, cbenc => self%cbenc, encmask => self%encmask) + + do i = 1, npl + if (nenc(i) == 0) cycle + do j = 1, nenc(i) + ! Copy the results of the integration back over + tpenc(i)%xh(:, j) = tpenc(i)%xh(:, j) + pl%xin(:, i, NTPHENC) + tpenc(i)%vh(:, j) = tpenc(i)%vh(:, j) + pl%vin(:, i, NTPHENC) + ! Put back heliocentric value of mu + call tpenc(i)%set_mu(cb) end do - end associate - deallocate(self%tpenc) - deallocate(self%cbenc) - deallocate(self%plenc) - deallocate(self%encmask) - - return + ! Replace the old test particle with the new one + call tp%fill(tpenc(i), encmask(:,i)) + end do + end associate + deallocate(self%tpenc) + deallocate(self%cbenc) + deallocate(self%plenc) + deallocate(self%encmask) + + return - end subroutine rmvs_step_end_planetocentric + end subroutine rmvs_step_end_planetocentric end submodule s_rmvs_step \ No newline at end of file diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 4fa5b1b84..7f2f26d7c 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -1,12 +1,11 @@ submodule (swiftest_classes) s_setup - + use swiftest contains module procedure setup_construct_system !! author: David A. Minton !! !! Constructor for a Swiftest nbody system. Creates the nbody system object based on the user-input integrator !! - use swiftest implicit none select case(config%integrator) case (BS) @@ -61,7 +60,6 @@ !! Constructor for base Swiftest particle class. Allocates space for all particles and !! initializes all components with a value. !! Note: Timing tests indicate that (NDIM, n) is more efficient than (NDIM, n) - use swiftest implicit none self%nbody = n @@ -113,7 +111,6 @@ !! !! Constructor for base Swiftest massive body class. Allocates space for all particles and !! initializes all components with a value. - use swiftest implicit none !> Call allocation method for parent class @@ -149,7 +146,6 @@ !! !! Constructor for base Swiftest test particle particle class. Allocates space for !! all particles and initializes all components with a value. - use swiftest implicit none !> Call allocation method for parent class @@ -181,7 +177,6 @@ !! author: David A. Minton !! !! Computes G * (M + m) for each massive body - use swiftest implicit none if (self%nbody > 0) self%mu(:) = cb%Gmass + self%Gmass(:) @@ -193,7 +188,6 @@ !! author: David A. Minton !! !! Converts certain scalar values to arrays so that they can be used in elemental functions - use swiftest implicit none if (self%nbody > 0) self%mu(:) = cb%Gmass @@ -205,7 +199,6 @@ !! author: David A. Minton !! !! Sets the value of the Hill's radius - use swiftest implicit none @@ -221,7 +214,6 @@ !! author: David A. Minton !! !! Sets the inverse heliocentric radius term (1/rh**3) for all bodies in a structure - use swiftest implicit none integer(I4B) :: i diff --git a/src/step/step.f90 b/src/step/step.f90 index 9c961ddd1..7fddc8225 100644 --- a/src/step/step.f90 +++ b/src/step/step.f90 @@ -1,5 +1,5 @@ submodule (swiftest_classes) s_step - + use swiftest contains module procedure step_system !! author: David A. Minton @@ -10,7 +10,6 @@ !! called internally. Before this point, the actual types of cb, pl, and tp are ambiguous. The select type constructs remove the !! ambiguity. !! - use swiftest implicit none select type(system => self) diff --git a/src/user/user_getacch.f90 b/src/user/user_getacch.f90 index 41f41aefd..d77f107b6 100644 --- a/src/user/user_getacch.f90 +++ b/src/user/user_getacch.f90 @@ -1,20 +1,20 @@ submodule(swiftest_classes) s_user_getacch + use swiftest contains module subroutine user_getacch_body(self, cb, config, t) - !! author: David A. Minton - !! - !! Add user-supplied heliocentric accelerations to planets - !! - !! Adapted from David E. Kaufmann's Swifter routine whm_user_getacch.f90 - use swiftest - implicit none - ! Arguments - class(swiftest_body), intent(inout) :: self !! Swiftest massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree - class(swiftest_configuration), intent(in) :: config !! Input collection of user configuration parameters - real(DP), intent(in) :: t !! Current time + !! author: David A. Minton + !! + !! Add user-supplied heliocentric accelerations to planets + !! + !! Adapted from David E. Kaufmann's Swifter routine whm_user_getacch.f90 + implicit none + ! Arguments + class(swiftest_body), intent(inout) :: self !! Swiftest massive body particle data structure + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structuree + class(swiftest_configuration), intent(in) :: config !! Input collection of user configuration parameters + real(DP), intent(in) :: t !! Current time - return + return end subroutine user_getacch_body end submodule s_user_getacch diff --git a/src/util/util_coord.f90 b/src/util/util_coord.f90 index affb895e8..f1c39b203 100644 --- a/src/util/util_coord.f90 +++ b/src/util/util_coord.f90 @@ -1,4 +1,5 @@ submodule(swiftest_classes) s_util_coord + use swiftest contains module subroutine util_coord_h2b_pl(self, cb) !! author: David A. Minton @@ -7,7 +8,6 @@ module subroutine util_coord_h2b_pl(self, cb) !! !! Adapted from David E. Kaufmann's Swifter routine coord_h2b.f90 !! Adapted from Hal Levison's Swift routine coord_h2b.f - use swiftest implicit none ! Arguments class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object @@ -38,10 +38,9 @@ module subroutine util_coord_h2b_tp(self, cb) !! !! Adapted from David E. Kaufmann's Swifter routine coord_h2b_tp.f90 !! Adapted from Hal Levison's Swift routine coord_h2b_tp.f - use swiftest implicit none ! Arguments - class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object + class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object class(swiftest_cb), intent(in) :: cb !! Swiftest central body object associate(ntp => self%nbody, xbcb => cb%xb, vbcb => cb%vb, status => self%status, & @@ -68,7 +67,6 @@ module subroutine util_coord_b2h_pl(self, cb) !! !! Adapted from David E. Kaufmann's Swifter routine coord_b2h.f90 !! Adapted from Hal Levison's Swift routine coord_b2h.f - use swiftest implicit none ! Arguments class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object @@ -94,7 +92,6 @@ module subroutine util_coord_b2h_tp(self, cb) !! !! Adapted from David E. Kaufmann's Swifter routine coord_b2h_tp.f90 !! Adapted from Hal Levison's Swift routine coord_b2h_tp.f - use swiftest implicit none ! Arguments class(swiftest_tp), intent(inout) :: self !! Swiftest massive body object diff --git a/src/util/util_copy.f90 b/src/util/util_copy.f90 index 83e676604..9731e9eff 100644 --- a/src/util/util_copy.f90 +++ b/src/util/util_copy.f90 @@ -1,12 +1,11 @@ submodule (swiftest_classes) s_util_copy - + use swiftest contains module procedure util_copy_cb !! author: David A. Minton !! !! Copies elements of one Swiftest central body object to another. The mask is ignored for this case !! - use swiftest implicit none select type(src) class is(swiftest_cb) @@ -35,7 +34,6 @@ !! !! Copies elements of one Swiftest generic body object to another. The elements are determined by a mask. !! - use swiftest implicit none select type(src) @@ -83,8 +81,7 @@ !! !! Copies elements of one Swiftest massive body object to another. The elements are determined by a mask. !! - use swiftest - + implicit none select type(src) class is(swiftest_pl) where (mask(:)) @@ -114,7 +111,6 @@ !! author: David A. Minton !! !! Copies elements of one swiftest test particle object to another. The elements are determined by a mask. - use swiftest implicit none select type(src) @@ -139,7 +135,6 @@ !! Copies elements of one Swiftest nbody system object to another. The elements are determined by a mask. !! In this case the mask contains the pl elements followed by the tp elements. !! - use swiftest implicit none associate(pl => self%pl, tp => self%tp, npl => self%pl%nbody, ntp => self%tp%nbody) diff --git a/src/util/util_exit.f90 b/src/util/util_exit.f90 index 410c81acc..ecf8d096d 100644 --- a/src/util/util_exit.f90 +++ b/src/util/util_exit.f90 @@ -1,4 +1,5 @@ submodule (util) s_util_exit + use swiftest contains module procedure util_exit !! author: David A. Minton @@ -7,7 +8,6 @@ !! !! Adapted from David E. Kaufmann's Swifter routine: util_exit.f90 !! Adapted from Hal Levison's Swift routine util_exit.f - use swiftest character(*), parameter :: BAR = '("------------------------------------------------")' select case(code) diff --git a/src/util/util_hills.f90 b/src/util/util_hills.f90 index f69e157e8..b743e1f31 100644 --- a/src/util/util_hills.f90 +++ b/src/util/util_hills.f90 @@ -1,4 +1,5 @@ submodule (util) s_util_hills + use swiftest contains module procedure util_hills !! author: David A. Minton @@ -7,7 +8,6 @@ !! !! Adapted from David E. Kaufmann's Swifter routine: util_hills.f90 !! Adapted from Hal Levison's Swift routine util_hills.f - use swiftest integer(I4B) :: i real(DP) :: msun, mp, mu, energy, ap, r, v2 diff --git a/src/util/util_index.f90 b/src/util/util_index.f90 index 90d54933a..4bc760ec2 100644 --- a/src/util/util_index.f90 +++ b/src/util/util_index.f90 @@ -1,4 +1,5 @@ submodule (util) s_util_index + use swiftest contains module procedure util_index !! author: David A. Minton @@ -8,7 +9,6 @@ !! Adapted from David E. Kaufmann's Swifter routine: util_index.f90 !! Adapted from Numerical Recipes in Fortran 90: The Art of Parallel Scientific Computing, by Press, Teukolsky, !! Vetterling, and Flannery, 2nd ed., pp. 1173-4 - use swiftest integer(I4B), parameter :: nn = 15, nstack = 50 integer(I4B) :: n, k, i, j, indext, jstack, l, r, dum integer(I4B), dimension(nstack) :: istack diff --git a/src/util/util_peri.f90 b/src/util/util_peri.f90 index bcd31719a..bcface016 100644 --- a/src/util/util_peri.f90 +++ b/src/util/util_peri.f90 @@ -1,4 +1,5 @@ submodule (util) s_util_peri + use swiftest contains module procedure util_peri !! author: David A. Minton @@ -9,7 +10,6 @@ !! !! Adapted from David E. Kaufmann's Swifter routine: util_peri.f90 !! Adapted from Hal Levison's Swift routine util_peri.f - use swiftest implicit none integer(I4B) :: i diff --git a/src/util/util_reverse_status.f90 b/src/util/util_reverse_status.f90 index b1d64cf45..253d10fde 100644 --- a/src/util/util_reverse_status.f90 +++ b/src/util/util_reverse_status.f90 @@ -1,11 +1,10 @@ submodule (swiftest_classes) s_util_reverse_status - + use swiftest contains module procedure util_reverse_status !! author: David A. Minton !! !! Reverses the active/inactive status of all particles in a structure - use swiftest implicit none where (self%status(:) == ACTIVE) diff --git a/src/util/util_sort_dp.f90 b/src/util/util_sort_dp.f90 index ddc39facd..5d02a82c1 100644 --- a/src/util/util_sort_dp.f90 +++ b/src/util/util_sort_dp.f90 @@ -1,4 +1,5 @@ submodule (util) s_util_sort_dp + use swiftest contains module procedure util_sort_dp !! author: David A. Minton @@ -8,7 +9,6 @@ !! Adapted from David E. Kaufmann's Swifter routine: util_sort_dp.f90 !! Adapted from Numerical Recipes in Fortran 90: The Art of Parallel Scientific Computing, by Press, Teukolsky, !! Vetterling, and Flannery, 2nd ed., pp. 1169-70 - use swiftest integer(I4B), parameter :: NN = 15, NSTACK = 50 real(DP) :: a, dum integer(I4B) :: n, k, i, j, jstack, l, r diff --git a/src/util/util_sort_i4b.f90 b/src/util/util_sort_i4b.f90 index 1af5ad42a..7c1e0d08b 100644 --- a/src/util/util_sort_i4b.f90 +++ b/src/util/util_sort_i4b.f90 @@ -1,4 +1,5 @@ submodule (util) s_util_sort_i4b + use swiftest contains module procedure util_sort_i4b !! author: David A. Minton @@ -8,7 +9,6 @@ !! Adapted from David E. Kaufmann's Swifter routine: util_sort_i4b.f90 !! Adapted from Numerical Recipes in Fortran 90: The Art of Parallel Scientific Computing, by Press, Teukolsky, !! Vetterling, and Flannery, 2nd ed., pp. 1169-70 - use swiftest integer(I4B), parameter :: NN = 15, NSTACK = 50 integer(I4B) :: a, n, k, i, j, jstack, l, r, dum integer(I4B), dimension(NSTACK) :: istack diff --git a/src/util/util_sort_sp.f90 b/src/util/util_sort_sp.f90 index 240063448..baa82e941 100644 --- a/src/util/util_sort_sp.f90 +++ b/src/util/util_sort_sp.f90 @@ -1,4 +1,5 @@ submodule (util) s_util_sort_sp + use swiftest contains module procedure util_sort_sp !! author: David A. Minton @@ -8,7 +9,6 @@ !! Adapted from David E. Kaufmann's Swifter routine: util_sort_DP.f90 !! Adapted from Numerical Recipes in Fortran 90: The Art of Parallel Scientific Computing, by Press, Teukolsky, !! Vetterling, and Flannery, 2nd ed., pp. 1169-70 - use swiftest integer(I4B), parameter :: NN = 15, NSTACK = 50 real(SP) :: a, dum integer(I4B) :: n, k, i, j, jstack, l, r diff --git a/src/util/util_spill_and_fill.f90 b/src/util/util_spill_and_fill.f90 index 8cf3c178c..25b269329 100644 --- a/src/util/util_spill_and_fill.f90 +++ b/src/util/util_spill_and_fill.f90 @@ -1,12 +1,11 @@ submodule (swiftest_classes) s_util_spill_and_fill -!! This submodule contains the methods spill, spill_body, spill_pl, and spill_tp for basic swiftest particles + use swiftest contains module subroutine util_spill_body(self, discards, lspill_list) !! author: David A. Minton !! !! Move spilled (discarded) Swiftest generic particle structure from active list to discard list !! Adapted from David E. Kaufmann's Swifter routine whm_discard_spill.f90 - use swiftest implicit none ! Arguments class(swiftest_body), intent(inout) :: self !! Swiftest generic body object @@ -69,7 +68,6 @@ module subroutine util_fill_body(self, inserts, lfill_list) !! !! Insert new Swiftest generic particle structure into an old one. !! This is the inverse of a fill operation. - use swiftest implicit none ! Arguments class(swiftest_body), intent(inout) :: self !! Swiftest generic body object @@ -145,7 +143,6 @@ end subroutine util_fill_body !! !! Move spilled (discarded) Swiftest massive body structure from active list to discard list !! Adapted from David E. Kaufmann's Swifter routine whm_discard_spill.f90 - use swiftest implicit none integer(I4B) :: i @@ -194,7 +191,6 @@ end subroutine util_fill_body !! !! Insert new Swiftest massive body structure into an old one. !! This is the inverse of a fill operation. - use swiftest implicit none integer(I4B) :: i @@ -247,7 +243,6 @@ end subroutine util_fill_body !! !! Move spilled (discarded) Swiftest test particle structure from active list to discard list !! Adapted from David E. Kaufmann's Swifter routine whm_discard_spill.f90 - use swiftest implicit none associate(keeps => self, ntp => self%nbody) @@ -270,34 +265,33 @@ end subroutine util_fill_body return end procedure util_spill_tp - module procedure util_fill_tp - !! author: David A. Minton - !! - !! Insert new Swiftest test particle structure into an old one. - !! This is the inverse of a fill operation. - use swiftest - implicit none - - associate(keeps => self) - select type(inserts) - class is (swiftest_tp) - !> Spill components specific to the test particle class - keeps%isperi(:) = unpack(keeps%isperi(:), .not.lfill_list(:), keeps%isperi(:)) - keeps%isperi(:) = unpack(inserts%isperi(:), lfill_list(:), keeps%isperi(:)) - - keeps%peri(:) = unpack(keeps%peri(:), .not.lfill_list(:), keeps%peri(:)) - keeps%peri(:) = unpack(inserts%peri(:), lfill_list(:), keeps%peri(:)) - - keeps%atp(:) = unpack(keeps%atp(:), .not.lfill_list(:), keeps%atp(:)) - keeps%atp(:) = unpack(inserts%atp(:), lfill_list(:), keeps%atp(:)) - - call util_fill_body(keeps, inserts, lfill_list) - class default - write(*,*) 'Error! fill method called for incompatible return type on swiftest_tp' - end select - end associate - return - end procedure util_fill_tp + module procedure util_fill_tp + !! author: David A. Minton + !! + !! Insert new Swiftest test particle structure into an old one. + !! This is the inverse of a fill operation. + implicit none + + associate(keeps => self) + select type(inserts) + class is (swiftest_tp) + !> Spill components specific to the test particle class + keeps%isperi(:) = unpack(keeps%isperi(:), .not.lfill_list(:), keeps%isperi(:)) + keeps%isperi(:) = unpack(inserts%isperi(:), lfill_list(:), keeps%isperi(:)) + + keeps%peri(:) = unpack(keeps%peri(:), .not.lfill_list(:), keeps%peri(:)) + keeps%peri(:) = unpack(inserts%peri(:), lfill_list(:), keeps%peri(:)) + + keeps%atp(:) = unpack(keeps%atp(:), .not.lfill_list(:), keeps%atp(:)) + keeps%atp(:) = unpack(inserts%atp(:), lfill_list(:), keeps%atp(:)) + + call util_fill_body(keeps, inserts, lfill_list) + class default + write(*,*) 'Error! fill method called for incompatible return type on swiftest_tp' + end select + end associate + return + end procedure util_fill_tp end submodule s_util_spill_and_fill diff --git a/src/util/util_toupper.f90 b/src/util/util_toupper.f90 index aac9f6840..f5bf8b979 100644 --- a/src/util/util_toupper.f90 +++ b/src/util/util_toupper.f90 @@ -1,4 +1,5 @@ submodule (util) s_util_toupper + use swiftest contains module procedure util_toupper !! author: David A. Minton @@ -6,7 +7,6 @@ !! Convert string to uppercase !! !! Adapted from David E. Kaufmann's Swifter routine: util_toupper.f90 - use swiftest integer(I4B) :: i, length, idx length = len(string) diff --git a/src/util/util_valid.f90 b/src/util/util_valid.f90 index afca40eb2..7dd755646 100644 --- a/src/util/util_valid.f90 +++ b/src/util/util_valid.f90 @@ -1,4 +1,5 @@ submodule (util) s_util_valid + use swiftest contains module procedure util_valid !! author: David A. Minton @@ -7,7 +8,6 @@ !! Subroutine causes program to exit with error if any ids are not unique !! !! Adapted from David E. Kaufmann's Swifter routine: util_valid.f90 - use swiftest integer(I4B) :: i integer(I4B), dimension(:), allocatable :: idarr diff --git a/src/util/util_version.f90 b/src/util/util_version.f90 index ed29c8307..0c0888f21 100644 --- a/src/util/util_version.f90 +++ b/src/util/util_version.f90 @@ -1,52 +1,52 @@ submodule (util) s_util_version + use swiftest contains module procedure util_version - !! author: David A. Minton - !! - !! Print program version information to terminale - !! - !! Adapted from David E. Kaufmann's Swifter routine: util_version.f90 - use swiftest - - write(*, 200) VERSION_NUMBER + !! author: David A. Minton + !! + !! Print program version information to terminale + !! + !! Adapted from David E. Kaufmann's Swifter routine: util_version.f90 + implicit none + write(*, 200) VERSION_NUMBER 200 format(/, "************* Swiftest: Version ", f3.1, " *************", //, & - "Based off of Swifter:", //, & - "Authors:", //, & - " The Purdue University Swiftest Development team ", /, & - " Lead by David A. Minton ", /, & - " Single loop blocking by Jacob R. Elliott", /, & - " Fragmentation by Carlisle A. Wishard and", //, & - " Jennifer L. L. Poutplin ", //, & - "Please address comments and questions to:", //, & - " David A. Minton", /, & - " Department Earth, Atmospheric, & Planetary Sciences ",/, & - " Purdue University", /, & - " 550 Stadium Mall Drive", /, & - " West Lafayette, Indiana 47907", /, & - " 765-250-8034 ", /, & - " daminton@purdue.edu", /, & - "Special thanks to Hal Levison and Martin Duncan for the original",/,& - "SWIFTER and SWIFT codes that made this possible.", //, & - "************************************************", /) + "Based off of Swifter:", //, & + "Authors:", //, & + " The Purdue University Swiftest Development team ", /, & + " Lead by David A. Minton ", /, & + " Single loop blocking by Jacob R. Elliott", /, & + " Fragmentation by Carlisle A. Wishard and", //, & + " Jennifer L. L. Poutplin ", //, & + "Please address comments and questions to:", //, & + " David A. Minton", /, & + " Department Earth, Atmospheric, & Planetary Sciences ",/, & + " Purdue University", /, & + " 550 Stadium Mall Drive", /, & + " West Lafayette, Indiana 47907", /, & + " 765-250-8034 ", /, & + " daminton@purdue.edu", /, & + "Special thanks to Hal Levison and Martin Duncan for the original",/,& + "SWIFTER and SWIFT codes that made this possible.", //, & + "************************************************", /) -100 FORMAT(/, "************* SWIFTER: Version ", F3.1, " *************", //, & - "Authors:", //, & - " Martin Duncan: Queen's University", /, & - " Hal Levison : Southwest Research Institute", //, & - "Please address comments and questions to:", //, & - " Hal Levison or David Kaufmann", /, & - " Department of Space Studies", /, & - " Southwest Research Institute", /, & - " 1050 Walnut Street, Suite 400", /, & - " Boulder, Colorado 80302", /, & - " 303-546-0290 (HFL), 720-240-0119 (DEK)", /, & - " 303-546-9687 (fax)", /, & - " hal@gort.boulder.swri.edu (HFL)", /, & - " kaufmann@boulder.swri.edu (DEK)", //, & - "************************************************", /) - - return + 100 FORMAT(/, "************* SWIFTER: Version ", F3.1, " *************", //, & + "Authors:", //, & + " Martin Duncan: Queen's University", /, & + " Hal Levison : Southwest Research Institute", //, & + "Please address comments and questions to:", //, & + " Hal Levison or David Kaufmann", /, & + " Department of Space Studies", /, & + " Southwest Research Institute", /, & + " 1050 Walnut Street, Suite 400", /, & + " Boulder, Colorado 80302", /, & + " 303-546-0290 (HFL), 720-240-0119 (DEK)", /, & + " 303-546-9687 (fax)", /, & + " hal@gort.boulder.swri.edu (HFL)", /, & + " kaufmann@boulder.swri.edu (DEK)", //, & + "************************************************", /) + return end procedure util_version + end submodule s_util_version diff --git a/src/whm/whm_coord.f90 b/src/whm/whm_coord.f90 index 69ffbaf0d..af5368aa8 100644 --- a/src/whm/whm_coord.f90 +++ b/src/whm/whm_coord.f90 @@ -1,4 +1,5 @@ submodule (whm_classes) s_whm_coord + use swiftest contains module subroutine whm_coord_h2j_pl(self, cb) !! author: David A. Minton @@ -10,7 +11,6 @@ module subroutine whm_coord_h2j_pl(self, cb) !! Adapted from David E. Kaufmann's Swifter routine coord_h2j.f90 !! !! Adapted from Hal Levison's Swift routine coord_h2j.f - use swiftest implicit none ! Arguments class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure @@ -50,7 +50,6 @@ module subroutine whm_coord_j2h_pl(self, cb) !! Adapted from David E. Kaufmann's Swifter routine coord_j2h.f90 !! !! Adapted from Hal Levison's Swift routine coord_j2h.f - use swiftest implicit none ! Arguments class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure @@ -87,7 +86,6 @@ module subroutine whm_coord_vh2vj_pl(self, cb) !! Adapted from David E. Kaufmann's Swifter routine coord_vh2vj.f90 !! !! Adapted from Hal Levison's Swift routine coord_vh2vj.f - use swiftest implicit none ! Arguments class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure diff --git a/src/whm/whm_drift.f90 b/src/whm/whm_drift.f90 index 4e7e7365d..e595edb7e 100644 --- a/src/whm/whm_drift.f90 +++ b/src/whm/whm_drift.f90 @@ -1,4 +1,5 @@ submodule(whm_classes) whm_drift + use swiftest contains module subroutine whm_drift_pl(self, cb, config, dt) !! author: David A. Minton @@ -7,7 +8,6 @@ module subroutine whm_drift_pl(self, cb, config, dt) !! !! Adapted from Hal Levison's Swift routine drift.f !! Adapted from David E. Kaufmann's Swifter routine whm_drift.f90 - use swiftest implicit none ! Arguments class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure @@ -64,7 +64,6 @@ module subroutine whm_drift_tp(self, cb, config, dt) !! Adapted from Hal Levison's Swift routine drift_tp.f !! Includes !! Adapted from David E. Kaufmann's Swifter routine whm_drift_tp.f90 - use swiftest implicit none ! Arguments class(whm_tp), intent(inout) :: self !! WHM test particle data structure diff --git a/src/whm/whm_getacch.f90 b/src/whm/whm_getacch.f90 index 9362a1a37..fb21789c0 100644 --- a/src/whm/whm_getacch.f90 +++ b/src/whm/whm_getacch.f90 @@ -1,44 +1,42 @@ submodule(whm_classes) s_whm_getacch + use swiftest contains module subroutine whm_getacch_pl(self, cb, config, t) - !! author: David A. Minton - !! - !! Compute heliocentric accelerations of planets - !! - !! Adapted from Hal Levison's Swift routine getacch.f - !! Adapted from David E. Kaufmann's Swifter routine whm_getacch.f90 - use swiftest - implicit none - ! Arguments - class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure - class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structure - class(swiftest_configuration), intent(in) :: config !! Input collection of - real(DP), intent(in) :: t !! Current time - ! Internals - integer(I4B) :: i - real(DP), dimension(:), allocatable, save :: fac - real(DP), dimension(NDIM) :: ah0 - real(DP) :: r2 - - associate(pl => self, npl => self%nbody, j2rp2 => cb%j2rp2, & - ah => self%ah, xh => self%xh, xj => self%xj, vh => self%vh, vj => self%vj) - if (npl == 0) return - call pl%set_ir3() - - ah0 = whm_getacch_ah0(pl%Gmass(2:npl), pl%xh(:,2:npl)) - do i = 1, NDIM - pl%ah(i, 1:npl) = ah0(i) - end do - call whm_getacch_ah1(cb, pl) - call whm_getacch_ah2(cb, pl) - call whm_getacch_ah3(pl) + !! author: David A. Minton + !! + !! Compute heliocentric accelerations of planets + !! + !! Adapted from Hal Levison's Swift routine getacch.f + !! Adapted from David E. Kaufmann's Swifter routine whm_getacch.f90 + implicit none + ! Arguments + class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure + class(swiftest_cb), intent(inout) :: cb !! Swiftest central body particle data structure + class(swiftest_configuration), intent(in) :: config !! Input collection of + real(DP), intent(in) :: t !! Current time + ! Internals + integer(I4B) :: i + real(DP), dimension(NDIM) :: ah0 + + associate(pl => self, npl => self%nbody, j2rp2 => cb%j2rp2, & + ah => self%ah, xh => self%xh, xj => self%xj, vh => self%vh, vj => self%vj) + if (npl == 0) return + call pl%set_ir3() + + ah0 = whm_getacch_ah0(pl%Gmass(2:npl), pl%xh(:,2:npl)) + do i = 1, NDIM + pl%ah(i, 1:npl) = ah0(i) + end do + call whm_getacch_ah1(cb, pl) + call whm_getacch_ah2(cb, pl) + call whm_getacch_ah3(pl) - if (config%loblatecb) call pl%obl_acc(cb) - if (config%lextra_force) call pl%user_getacch(cb, config, t) - if (config%lgr) call pl%gr_getacch(cb, config) + if (config%loblatecb) call pl%obl_acc(cb) + if (config%lextra_force) call pl%user_getacch(cb, config, t) + if (config%lgr) call pl%gr_getacch(cb, config) - end associate - return + end associate + return end subroutine whm_getacch_pl module subroutine whm_getacch_tp(self, cb, pl, config, t, xh) @@ -48,7 +46,6 @@ module subroutine whm_getacch_tp(self, cb, pl, config, t, xh) !! !! Adapted from Hal Levison's Swift routine getacch_tp.f !! Adapted from David E. Kaufmann's Swifter routine whm_getacch_tp.f90 - use swiftest implicit none ! Arguments class(whm_tp), intent(inout) :: self !! WHM test particle data structure @@ -59,9 +56,7 @@ module subroutine whm_getacch_tp(self, cb, pl, config, t, xh) real(DP), dimension(:,:), intent(in) :: xh !! Heliocentric positions of planets ! Internals integer(I4B) :: i - real(DP), dimension(:), allocatable, save :: fac real(DP), dimension(NDIM) :: ah0 - real(DP) :: r2 associate(tp => self, ntp => self%nbody, npl => pl%nbody, j2rp2 => cb%j2rp2, aht => self%ah, & ir3h => pl%ir3h, GMpl => pl%Gmass) @@ -83,7 +78,6 @@ pure function whm_getacch_ah0(mu, xh) result(ah0) !! author: David A. Minton !! !! Compute zeroth term heliocentric accelerations of planets - use swiftest implicit none ! Arguments real(DP), dimension(:), intent(in) :: mu @@ -114,7 +108,6 @@ pure subroutine whm_getacch_ah1(cb, pl) !! !! Adapted from Hal Levison's Swift routine getacch_ah1.f !! Adapted from David E. Kaufmann's Swifter routine whm_getacch_ah1.f90 - use swiftest implicit none ! Arguments class(swiftest_cb), intent(in) :: cb !! Swiftest central body object @@ -142,7 +135,6 @@ pure subroutine whm_getacch_ah2(cb, pl) !! !! Adapted from Hal Levison's Swift routine getacch_ah2.f !! Adapted from David E. Kaufmann's Swifter routine whm_getacch_ah2.f90 - use swiftest implicit none ! Arguments class(swiftest_cb), intent(in) :: cb !! Swiftest central body object @@ -175,7 +167,6 @@ pure subroutine whm_getacch_ah3(pl) !! !! Adapted from Hal Levison's Swift routine getacch_ah3.f !! Adapted from David E. Kaufmann's Swifter routine whm_getacch_ah3.f90 - use swiftest implicit none class(whm_pl), intent(inout) :: pl @@ -215,7 +206,6 @@ pure subroutine whm_getacch_ah3_tp(cb, pl, tp, xh) !! !! Adapted from Hal Levison's Swift routine getacch_ah3_tp.f !! Adapted from David E. Kaufmann's Swifter routine whm_getacch_ah3.f90 - use swiftest implicit none ! Arguments class(swiftest_cb), intent(in) :: cb !! Swiftest central body object diff --git a/src/whm/whm_gr.f90 b/src/whm/whm_gr.f90 index a6f8901cf..edd6ee26d 100644 --- a/src/whm/whm_gr.f90 +++ b/src/whm/whm_gr.f90 @@ -1,4 +1,5 @@ submodule(whm_classes) s_whm_gr + use swiftest contains module subroutine whm_gr_getacch_pl(self, cb, config) !! author: David A. Minton @@ -7,7 +8,6 @@ module subroutine whm_gr_getacch_pl(self, cb, config) !! Based on Saha & Tremaine (1994) Eq. 28 !! !! Adapted from David A. Minton's Swifter routine routine gr_whm_getacch.f90 - use swiftest implicit none ! Arguments class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure @@ -46,7 +46,6 @@ module subroutine whm_gr_getacch_tp(self, cb, config) !! Based on Saha & Tremaine (1994) Eq. 28 !! !! Adapted from David A. Minton's Swifter routine routine gr_whm_getacch.f90 - use swiftest implicit none ! Arguments class(whm_tp), intent(inout) :: self !! WHM massive body particle data structure @@ -76,7 +75,6 @@ module pure subroutine whm_gr_p4_pl(self, config, dt) !! Based on Saha & Tremaine (1994) Eq. 28 !! !! Adapted from David A. Minton's Swifter routine routine gr_whm_p4.f90 - use swiftest implicit none ! Arguments class(whm_pl), intent(inout) :: self !! Swiftest particle object @@ -103,7 +101,6 @@ module pure subroutine whm_gr_p4_tp(self, config, dt) !! Based on Saha & Tremaine (1994) Eq. 28 !! !! Adapted from David A. Minton's Swifter routine routine gr_whm_p4.f90 - use swiftest implicit none ! Arguments class(whm_tp), intent(inout) :: self !! Swiftest particle object @@ -128,14 +125,12 @@ module pure subroutine whm_gr_pv2vh_pl(self, config) !! !! Wrapper function that converts from pseudovelocity to heliocentric velocity for massive bodies !! in a WHM object - use swiftest implicit none ! Arguments class(whm_pl), intent(inout) :: self !! Swiftest particle object class(swiftest_configuration), intent(in) :: config !! Input collection of on parameters ! Internals integer(I4B) :: i - real(DP), dimension(:), allocatable :: mu real(DP), dimension(:,:), allocatable :: vh !! Temporary holder of pseudovelocity for in-place conversion associate(n => self%nbody, xh => self%xh, pv => self%vh, status => self%status, mu => self%muj) @@ -156,14 +151,12 @@ module pure subroutine whm_gr_pv2vh_tp(self, config) !! !! Wrapper function that converts from pseudovelocity to heliocentric velocity for test particles bodies !! in a WHM object - use swiftest implicit none ! Arguments class(whm_tp), intent(inout) :: self !! Swiftest particle object class(swiftest_configuration), intent(in) :: config !! Input collection of on parameters ! Internals integer(I4B) :: i - real(DP), dimension(:), allocatable :: mu real(DP), dimension(:,:), allocatable :: vh !! Temporary holder of pseudovelocity for in-place conversion associate(n => self%nbody, xh => self%xh, pv => self%vh, status => self%status, mu => self%mu) @@ -184,14 +177,12 @@ module pure subroutine whm_gr_vh2pv_pl(self, config) !! !! Wrapper function that converts from heliocentric velocity to pseudovelocity for massive bodies !! in a WHM object - use swiftest implicit none ! Arguments class(whm_pl), intent(inout) :: self !! Swiftest particle object class(swiftest_configuration), intent(in) :: config !! Input collection of on parameters ! Internals integer(I4B) :: i - real(DP), dimension(:), allocatable :: mu real(DP), dimension(:,:), allocatable :: pv !! Temporary holder of pseudovelocity for in-place conversion associate(n => self%nbody, xh => self%xh, vh => self%vh, status => self%status, mu => self%muj) @@ -212,14 +203,12 @@ module pure subroutine whm_gr_vh2pv_tp(self, config) !! !! Wrapper function that converts from heliocentric velocity to pseudovelocity for teset particles !! in a WHM object - use swiftest implicit none ! Arguments class(whm_tp), intent(inout) :: self !! Swiftest particle object class(swiftest_configuration), intent(in) :: config !! Input collection of on parameters ! Internals integer(I4B) :: i - real(DP), dimension(:), allocatable :: mu real(DP), dimension(:,:), allocatable :: pv !! Temporary holder of pseudovelocity for in-place conversion associate(n => self%nbody, xh => self%xh, vh => self%vh, status => self%status, mu => self%mu) @@ -243,7 +232,6 @@ pure subroutine gr_vel2pseudovel(config, mu, xh, vh, pv) !! this is only done once per run). !! !! Adapted from David A. Minton's Swifter routine gr_vel2pseudovel.f90 - use swiftest implicit none class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters @@ -317,7 +305,6 @@ pure subroutine gr_pseudovel2vel(config, mu, xh, pv, vh) !! Based on Saha & Tremaine (1994) Eq. 32 !! !! Adapted from David A. Minton's Swifter routine gr_pseudovel2vel.f90 - use swiftest implicit none class(swiftest_configuration), intent(in) :: config !! Input collection of configuration parameters real(DP), intent(in) :: mu !! G * (Mcb + m), G = gravitational constant, Mcb = mass of central body, m = mass of body @@ -337,14 +324,12 @@ pure subroutine gr_pseudovel2vel(config, mu, xh, pv, vh) end subroutine gr_pseudovel2vel pure subroutine p4_func(x, v, dt, c2) - use swiftest implicit none real(DP), dimension(:), intent(inout) :: x real(DP), dimension(:), intent(in) :: v real(DP), intent(in) :: dt, c2 real(DP), dimension(NDIM) :: dr real(DP) :: vmag2 - integer(I4B) :: i vmag2 = dot_product(v(:), v(:)) dr(:) = -2 * c2 * vmag2 * v(:) diff --git a/src/whm/whm_setup.f90 b/src/whm/whm_setup.f90 index 4ef7397eb..3db8ee067 100644 --- a/src/whm/whm_setup.f90 +++ b/src/whm/whm_setup.f90 @@ -1,4 +1,5 @@ submodule(whm_classes) s_whm_setup + use swiftest contains module subroutine whm_setup_pl(self,n) !! author: David A. Minton @@ -6,7 +7,6 @@ module subroutine whm_setup_pl(self,n) !! Allocate WHM planet structure !! !! Equivalent in functionality to David E. Kaufmann's Swifter routine whm_setup.f90 - use swiftest implicit none ! Arguments class(whm_pl), intent(inout) :: self !! Swiftest test particle object @@ -36,7 +36,6 @@ module subroutine whm_setup_tp(self,n) !! Allocate WHM test particle structure !! !! Equivalent in functionality to David E. Kaufmann's Swifter routine whm_setup.f90 - use swiftest implicit none ! Arguments class(whm_tp), intent(inout) :: self !! WHM test particle data structure @@ -52,7 +51,6 @@ module subroutine whm_setup_set_mu_eta_pl(self, cb) !! author: David A. Minton !! !! Sets the Jacobi mass value eta for all massive bodies - use swiftest implicit none ! Arguments class(whm_pl), intent(inout) :: self !! Swiftest system object @@ -79,7 +77,6 @@ module subroutine whm_setup_system(self, config) !! !! Wrapper method to initialize a basic Swiftest nbody system from files !! - use swiftest implicit none ! Arguments class(whm_nbody_system), intent(inout) :: self !! Swiftest system object @@ -111,7 +108,6 @@ module subroutine whm_setup_set_ir3j(self) !! author: David A. Minton !! !! Sets the inverse Jacobi and heliocentric radii cubed (1/rj**3 and 1/rh**3) - use swiftest implicit none ! Arguments class(whm_pl), intent(inout) :: self !! WHM massive body object @@ -131,12 +127,10 @@ module subroutine whm_setup_set_ir3j(self) end if end subroutine whm_setup_set_ir3j - module subroutine whm_setup_set_beg_end(self, xbeg, xend, vbeg) !! author: David A. Minton !! !! Sets one or more of the values of xbeg and xend - use swiftest implicit none ! Arguments class(whm_tp), intent(inout) :: self !! Swiftest test particle object diff --git a/src/whm/whm_spill_and_fill.f90 b/src/whm/whm_spill_and_fill.f90 index 0f2f30b7a..f5edf894c 100644 --- a/src/whm/whm_spill_and_fill.f90 +++ b/src/whm/whm_spill_and_fill.f90 @@ -1,4 +1,5 @@ submodule(whm_classes) s_whm_spill_and_fill + use swiftest contains module subroutine whm_spill_pl(self, discards, lspill_list) !! author: David A. Minton @@ -6,7 +7,6 @@ module subroutine whm_spill_pl(self, discards, lspill_list) !! Move spilled (discarded) WHM test particle structure from active list to discard list !! !! Adapted from David E. Kaufmann's Swifter routine whm_discard_spill.f90 - use swiftest implicit none ! Arguments class(whm_pl), intent(inout) :: self !! WHM massive body object @@ -51,7 +51,6 @@ module subroutine whm_fill_pl(self, inserts, lfill_list) !! This is the inverse of a fill operation. !! !! Adapted from David E. Kaufmann's Swifter routine whm_discard_spill.f90 - use swiftest implicit none ! Arguments class(whm_pl), intent(inout) :: self !! WHM massive body object diff --git a/src/whm/whm_step.f90 b/src/whm/whm_step.f90 index 83c637d00..120a7194b 100644 --- a/src/whm/whm_step.f90 +++ b/src/whm/whm_step.f90 @@ -1,4 +1,5 @@ submodule(whm_classes) s_whm_step + use swiftest contains module subroutine whm_step_system(cb, pl, tp, config) !! author: David A. Minton @@ -7,7 +8,6 @@ module subroutine whm_step_system(cb, pl, tp, config) !! !! Adapted from Hal Levison's Swift routine step_kdk.f !! Adapted from David E. Kaufmann's Swifter routine whm_step.f90 - use swiftest implicit none ! Arguments class(whm_cb), intent(inout) :: cb !! Swiftest central body object @@ -34,7 +34,6 @@ module subroutine whm_step_pl(self, cb, config, t, dt) !! Adapted from Hal Levison's Swift routine step_kdk_pl.f !! Adapted from David E. Kaufmann's Swifter routine whm_step_pl.f90 !logical, save :: lfirst = .true. - use swiftest implicit none ! Arguments class(whm_pl), intent(inout) :: self !! WHM massive body particle data structure @@ -74,7 +73,6 @@ module subroutine whm_step_tp(self, cb, pl, config, t, dt) !! !! Adapted from Hal Levison's Swift routine step_kdk_tp.f !! Adapted from David E. Kaufmann's Swifter routine whm_step_tp.f90 - use swiftest implicit none ! Arguments class(whm_tp), intent(inout) :: self !! WHM test particle data structure