diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90 index dadde1dfd..befc4ba8a 100644 --- a/src/collision/collision_resolve.f90 +++ b/src/collision/collision_resolve.f90 @@ -350,7 +350,7 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status) class is (swiftest_nbody_system) select type(param) class is (swiftest_parameters) - associate(pl => nbody_system%pl, pl_discards => nbody_system%pl_discards, info => nbody_system%pl%info, pl_adds => nbody_system%pl_adds, cb => nbody_system%cb, npl => nbody_system%pl%nbody, & + associate(pl => nbody_system%pl, pl_discards => nbody_system%pl_discards, info => nbody_system%pl%info, pl_adds => nbody_system%pl_adds, cb => nbody_system%cb, & collider => nbody_system%collider, impactors => nbody_system%collider%impactors,fragments => nbody_system%collider%fragments) ! Add the impactors%id bodies to the subtraction list diff --git a/src/collision/collision_util.f90 b/src/collision/collision_util.f90 index f4f8e341b..45f4f516d 100644 --- a/src/collision/collision_util.f90 +++ b/src/collision/collision_util.f90 @@ -536,7 +536,7 @@ module subroutine collision_util_set_coordinate_collider(self) ! Arguments class(collision_basic), intent(inout) :: self !! Collisional nbody_system - associate(fragments => self%fragments, impactors => self%impactors, nfrag => self%fragments%nbody) + associate(fragments => self%fragments, impactors => self%impactors) call impactors%set_coordinate_system() if (.not.allocated(self%fragments)) return diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 93b4c3431..f1ba8d606 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -350,7 +350,7 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu real(DP), dimension(2) :: fragment_cloud_radius logical, dimension(collider%fragments%nbody) :: loverlap real(DP), dimension(collider%fragments%nbody) :: mass_rscale, phi, theta, u - integer(I4B) :: i, j, loop, istart + integer(I4B) :: i, j, loop, istart, nfrag, npl, ntp logical :: lsupercat, lhitandrun integer(I4B), parameter :: MAXLOOP = 10000 real(DP), parameter :: cloud_size_scale_factor = 3.0_DP ! Scale factor to apply to the size of the cloud relative to the distance from the impact point. @@ -358,8 +358,10 @@ module subroutine fraggle_generate_pos_vec(collider, nbody_system, param, lfailu real(DP), parameter :: rbuffer = 1.01_DP ! Body radii are inflated by this scale factor to prevent secondary collisions real(DP), parameter :: pack_density = 0.5236_DP ! packing density of loose spheres - associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody, & - pl => nbody_system%pl, tp => nbody_system%tp, npl => nbody_system%pl%nbody, ntp => nbody_system%tp%nbody) + associate(fragments => collider%fragments, impactors => collider%impactors, pl => nbody_system%pl, tp => nbody_system%tp) + nfrag = collider%fragments%nbody + npl = nbody_system%pl%nbody + ntp = nbody_system%tp%nbody lsupercat = (impactors%regime == COLLRESOLVE_REGIME_SUPERCATASTROPHIC) lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) @@ -508,7 +510,7 @@ module subroutine fraggle_generate_rot_vec(collider, nbody_system, param) class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters ! Internals - integer(I4B) :: i + integer(I4B) :: i, nfrag real(DP), parameter :: frag_rot_fac = 0.1_DP ! Fraction of projectile rotation magnitude to add as random noise to fragment rotation real(DP), parameter :: hitandrun_momentum_transfer = 0.01_DP ! Fraction of projectile momentum transfered to target in a hit and run real(DP) :: mass_fac @@ -516,7 +518,8 @@ module subroutine fraggle_generate_rot_vec(collider, nbody_system, param) integer(I4B), parameter :: MAXLOOP = 10 logical :: lhitandrun - associate(fragments => collider%fragments, impactors => collider%impactors, nfrag => collider%fragments%nbody) + associate(fragments => collider%fragments, impactors => collider%impactors) + nfrag = collider%fragments%nbody lhitandrun = (impactors%regime == COLLRESOLVE_REGIME_HIT_AND_RUN) ! Initialize fragment rotations and velocities to be pre-impact rotation for body 1, and randomized for bodies >1 and scaled to the original rotation. @@ -805,7 +808,7 @@ module subroutine fraggle_generate_vel_vec(collider, nbody_system, param, lfailu L_residual(:) = (collider%L_total(:,2) - collider%L_total(:,1)) call collision_util_velocity_torque(-L_residual(:), collider%fragments%mtot, impactors%rbcom, impactors%vbcom) - do concurrent(i = 1:collider%fragments%nbody) + do concurrent(i = 1:nfrag) collider%fragments%vb(:,i) = collider%fragments%vc(:,i) + impactors%vbcom(:) end do diff --git a/src/helio/helio_gr.f90 b/src/helio/helio_gr.f90 index 2092fca6a..def25224d 100644 --- a/src/helio/helio_gr.f90 +++ b/src/helio/helio_gr.f90 @@ -71,11 +71,12 @@ pure module subroutine helio_gr_p4_pl(self, nbody_system, param, dt) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters real(DP), intent(in) :: dt !! Step size ! Internals - integer(I4B) :: i + integer(I4B) :: i, npl if (self%nbody == 0) return - associate(pl => self, npl => self%nbody) + associate(pl => self) + npl = self%nbody do concurrent(i = 1:npl, pl%lmask(i)) call swiftest_gr_p4_pos_kick(param, pl%rh(:, i), pl%vb(:, i), dt) end do @@ -99,11 +100,12 @@ pure module subroutine helio_gr_p4_tp(self, nbody_system, param, dt) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters real(DP), intent(in) :: dt !! Step size ! Internals - integer(I4B) :: i + integer(I4B) :: i, ntp if (self%nbody == 0) return - associate(tp => self, ntp => self%nbody) + associate(tp => self) + ntp = self%nbody do concurrent(i = 1:ntp, tp%lmask(i)) call swiftest_gr_p4_pos_kick(param, tp%rh(:, i), tp%vb(:, i), dt) end do diff --git a/src/helio/helio_kick.f90 b/src/helio/helio_kick.f90 index 6fb4fad43..89c93a7ed 100644 --- a/src/helio/helio_kick.f90 +++ b/src/helio/helio_kick.f90 @@ -104,11 +104,12 @@ module subroutine helio_kick_vb_pl(self, nbody_system, param, t, dt, lbeg) real(DP), intent(in) :: dt !! Stepsize logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. ! Internals - integer(I4B) :: i + integer(I4B) :: i, npl if (self%nbody == 0) return - associate(pl => self, npl => self%nbody) + associate(pl => self) + npl = self%nbody pl%ah(:, 1:npl) = 0.0_DP call pl%accel(nbody_system, param, t, lbeg) if (lbeg) then @@ -143,11 +144,12 @@ module subroutine helio_kick_vb_tp(self, nbody_system, param, t, dt, lbeg) real(DP), intent(in) :: dt !! Stepsize logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. ! Internals - integer(I4B) :: i + integer(I4B) :: i, ntp if (self%nbody == 0) return - associate(tp => self, ntp => self%nbody) + associate(tp => self) + ntp = self%nbody tp%ah(:, 1:ntp) = 0.0_DP call tp%accel(nbody_system, param, t, lbeg) do concurrent(i = 1:ntp, tp%lmask(i)) diff --git a/src/rmvs/rmvs_kick.f90 b/src/rmvs/rmvs_kick.f90 index 7d113f863..8fb6c14ce 100644 --- a/src/rmvs/rmvs_kick.f90 +++ b/src/rmvs/rmvs_kick.f90 @@ -29,11 +29,13 @@ module subroutine rmvs_kick_getacch_tp(self, nbody_system, param, t, lbeg) class(swiftest_parameters), allocatable :: param_planetocen real(DP), dimension(:, :), allocatable :: rh_original real(DP) :: GMcb_original - integer(I4B) :: i + integer(I4B) :: i, ntp, inner_index if (self%nbody == 0) return - associate(tp => self, ntp => self%nbody, ipleP => self%ipleP, inner_index => self%index) + associate(tp => self, ipleP => self%ipleP) + ntp = self%nbody + inner_index = self%index select type(nbody_system) class is (rmvs_nbody_system) if (nbody_system%lplanetocentric) then ! This is a close encounter step, so any accelerations requiring heliocentric position values diff --git a/src/swiftest/swiftest_obl.f90 b/src/swiftest/swiftest_obl.f90 index 6bd7480fb..7a6463677 100644 --- a/src/swiftest/swiftest_obl.f90 +++ b/src/swiftest/swiftest_obl.f90 @@ -73,11 +73,12 @@ module subroutine swiftest_obl_acc_pl(self, nbody_system) class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object ! Internals - integer(I4B) :: i + integer(I4B) :: i, npl if (self%nbody == 0) return - associate(pl => self, npl => self%nbody, cb => nbody_system%cb) + associate(pl => self, cb => nbody_system%cb) + npl = self%nbody call swiftest_obl_acc(npl, cb%Gmass, cb%j2rp2, cb%j4rp4, pl%rh, pl%lmask, pl%aobl, pl%Gmass, cb%aobl) do concurrent(i = 1:npl, pl%lmask(i)) @@ -103,11 +104,12 @@ module subroutine swiftest_obl_acc_tp(self, nbody_system) class(swiftest_nbody_system), intent(inout) :: nbody_system !! Swiftest nbody system object ! Internals real(DP), dimension(NDIM) :: aoblcb - integer(I4B) :: i + integer(I4B) :: i, ntp if (self%nbody == 0) return - associate(tp => self, ntp => self%nbody, cb => nbody_system%cb) + associate(tp => self, cb => nbody_system%cb) + ntp = self%nbody call swiftest_obl_acc(ntp, cb%Gmass, cb%j2rp2, cb%j4rp4, tp%rh, tp%lmask, tp%aobl) if (nbody_system%lbeg) then aoblcb = cb%aoblbeg @@ -139,10 +141,11 @@ module subroutine swiftest_obl_pot_system(self) ! Arguments class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object ! Internals - integer(I4B) :: i + integer(I4B) :: i, npl real(DP), dimension(self%pl%nbody) :: oblpot_arr - associate(nbody_system => self, pl => self%pl, npl => self%pl%nbody, cb => self%cb) + associate(nbody_system => self, pl => self%pl, cb => self%cb) + npl = self%pl%nbody if (.not. any(pl%lmask(1:npl))) return do concurrent (i = 1:npl, pl%lmask(i)) oblpot_arr(i) = swiftest_obl_pot_one(cb%Gmass, pl%Gmass(i), cb%j2rp2, cb%j4rp4, pl%rh(3,i), 1.0_DP / norm2(pl%rh(:,i))) diff --git a/src/swiftest/swiftest_util.f90 b/src/swiftest/swiftest_util.f90 index b620a7f9d..1784a1252 100644 --- a/src/swiftest/swiftest_util.f90 +++ b/src/swiftest/swiftest_util.f90 @@ -277,10 +277,11 @@ module subroutine swiftest_util_coord_h2b_tp(self, cb) class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object class(swiftest_cb), intent(in) :: cb !! Swiftest central body object ! Internals - integer(I4B) :: i + integer(I4B) :: i, ntp if (self%nbody == 0) return - associate(tp => self, ntp => self%nbody) + associate(tp => self) + ntp = self%nbody do concurrent (i = 1:ntp, tp%status(i) /= INACTIVE) tp%rb(:, i) = tp%rh(:, i) + cb%rb(:) tp%vb(:, i) = tp%vh(:, i) + cb%vb(:) @@ -303,11 +304,12 @@ module subroutine swiftest_util_coord_b2h_pl(self, cb) class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object ! Internals - integer(I4B) :: i + integer(I4B) :: i, npl if (self%nbody == 0) return - associate(pl => self, npl => self%nbody) + associate(pl => self) + npl = self%nbody do concurrent (i = 1:npl, pl%status(i) /= INACTIVE) pl%rh(:, i) = pl%rb(:, i) - cb%rb(:) pl%vh(:, i) = pl%vb(:, i) - cb%vb(:) @@ -330,11 +332,12 @@ module subroutine swiftest_util_coord_b2h_tp(self, cb) class(swiftest_tp), intent(inout) :: self !! Swiftest massive body object class(swiftest_cb), intent(in) :: cb !! Swiftest central body object ! Internals - integer(I4B) :: i + integer(I4B) :: i, ntp if (self%nbody == 0) return - associate(tp => self, ntp => self%nbody) + associate(tp => self) + ntp = self%nbody do concurrent(i = 1:ntp, tp%status(i) /= INACTIVE) tp%rh(:, i) = tp%rb(:, i) - cb%rb(:) tp%vh(:, i) = tp%vb(:, i) - cb%vb(:) @@ -357,11 +360,12 @@ module subroutine swiftest_util_coord_vb2vh_pl(self, cb) class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object ! Internals - integer(I4B) :: i + integer(I4B) :: i, npl if (self%nbody == 0) return - associate(pl => self, npl => self%nbody) + associate(pl => self) + npl = self%nbody cb%vb(:) = 0.0_DP do i = npl, 1, -1 if (pl%status(i) /= INACTIVE) cb%vb(:) = cb%vb(:) - pl%Gmass(i) * pl%vb(:, i) / cb%Gmass @@ -413,12 +417,13 @@ module subroutine swiftest_util_coord_vh2vb_pl(self, cb) class(swiftest_pl), intent(inout) :: self !! Swiftest massive body object class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object ! Internals - integer(I4B) :: i + integer(I4B) :: i, npl real(DP) :: Gmtot if (self%nbody == 0) return - associate(pl => self, npl => self%nbody) + associate(pl => self) + npl = self%nbody Gmtot = cb%Gmass + sum(pl%Gmass(1:npl)) cb%vb(:) = 0.0_DP do i = 1, npl @@ -508,10 +513,11 @@ module subroutine swiftest_util_coord_rh2rb_tp(self, cb) class(swiftest_tp), intent(inout) :: self !! Swiftest test particle object class(swiftest_cb), intent(in) :: cb !! Swiftest central body object ! Internals - integer(I4B) :: i + integer(I4B) :: i, ntp if (self%nbody == 0) return - associate(tp => self, ntp => self%nbody) + associate(tp => self) + ntp = self%nbody do concurrent (i = 1:ntp, tp%status(i) /= INACTIVE) tp%rb(:, i) = tp%rh(:, i) + cb%rb(:) end do @@ -1145,7 +1151,7 @@ module subroutine swiftest_util_get_energy_and_momentum_system(self, param) class(swiftest_nbody_system), intent(inout) :: self !! Swiftest nbody system object class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters ! Internals - integer(I4B) :: i,j + integer(I4B) :: i,j, npl real(DP) :: kecb, kespincb real(DP), dimension(self%pl%nbody) :: kepl, kespinpl real(DP), dimension(NDIM,self%pl%nbody) :: Lplorbit @@ -1153,7 +1159,8 @@ module subroutine swiftest_util_get_energy_and_momentum_system(self, param) real(DP), dimension(NDIM) :: Lcborbit, Lcbspin real(DP), dimension(NDIM) :: h - associate(nbody_system => self, pl => self%pl, npl => self%pl%nbody, cb => self%cb) + associate(nbody_system => self, pl => self%pl, cb => self%cb) + npl = self%pl%nbody nbody_system%L_orbit(:) = 0.0_DP nbody_system%L_spin(:) = 0.0_DP nbody_system%L_total(:) = 0.0_DP diff --git a/src/whm/whm_gr.f90 b/src/whm/whm_gr.f90 index b0891f006..8eb69abc6 100644 --- a/src/whm/whm_gr.f90 +++ b/src/whm/whm_gr.f90 @@ -84,10 +84,12 @@ pure module subroutine whm_gr_p4_pl(self, nbody_system, param, dt) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters real(DP), intent(in) :: dt !! Step size ! Internals - integer(I4B) :: i + integer(I4B) :: i, npl - associate(pl => self, npl => self%nbody) - if (npl == 0) return + if (self%nbody == 0) return + + associate(pl => self) + npl = self%nbody do concurrent(i = 1:npl, pl%lmask(i)) call swiftest_gr_p4_pos_kick(param, pl%xj(:, i), pl%vj(:, i), dt) end do @@ -111,9 +113,10 @@ pure module subroutine whm_gr_p4_tp(self, nbody_system, param, dt) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters real(DP), intent(in) :: dt !! Step size ! Internals - integer(I4B) :: i + integer(I4B) :: i, ntp - associate(tp => self, ntp => self%nbody) + associate(tp => self) + ntp = self%nbody if (ntp == 0) return do concurrent(i = 1:ntp, tp%lmask(i)) call swiftest_gr_p4_pos_kick(param, tp%rh(:, i), tp%vh(:, i), dt) diff --git a/src/whm/whm_kick.f90 b/src/whm/whm_kick.f90 index 403678ed6..a1d4b0dc6 100644 --- a/src/whm/whm_kick.f90 +++ b/src/whm/whm_kick.f90 @@ -82,10 +82,12 @@ module subroutine whm_kick_getacch_tp(self, nbody_system, param, t, lbeg) real(DP), intent(in) :: t !! Current time logical, intent(in) :: lbeg !! Logical flag that determines whether or not this is the beginning or end of the step ! Internals - integer(I4B) :: i + integer(I4B) :: i, npl, ntp real(DP), dimension(NDIM) :: ah0 - associate(tp => self, ntp => self%nbody, pl => nbody_system%pl, cb => nbody_system%cb, npl => nbody_system%pl%nbody) + associate(tp => self, pl => nbody_system%pl, cb => nbody_system%cb) + npl = nbody_system%pl%nbody + ntp = self%nbody if (ntp == 0 .or. npl == 0) return nbody_system%lbeg = lbeg @@ -151,16 +153,15 @@ pure subroutine whm_kick_getacch_ah1(cb, pl) class(swiftest_cb), intent(in) :: cb !! WHM central body object class(whm_pl), intent(inout) :: pl !! WHM massive body object ! Internals - integer(I4B) :: i + integer(I4B) :: i, npl real(DP), dimension(NDIM) :: ah1h, ah1j - associate(npl => pl%nbody) - do concurrent (i = 2:npl, pl%lmask(i)) - ah1j(:) = pl%xj(:, i) * pl%ir3j(i) - ah1h(:) = pl%rh(:, i) * pl%ir3h(i) - pl%ah(:, i) = pl%ah(:, i) + cb%Gmass * (ah1j(:) - ah1h(:)) - end do - end associate + npl = pl%nbody + do concurrent (i = 2:npl, pl%lmask(i)) + ah1j(:) = pl%xj(:, i) * pl%ir3j(i) + ah1h(:) = pl%rh(:, i) * pl%ir3h(i) + pl%ah(:, i) = pl%ah(:, i) + cb%Gmass * (ah1j(:) - ah1h(:)) + end do return end subroutine whm_kick_getacch_ah1 @@ -178,22 +179,21 @@ pure subroutine whm_kick_getacch_ah2(cb, pl) class(swiftest_cb), intent(in) :: cb !! Swiftest central body object class(whm_pl), intent(inout) :: pl !! WHM massive body object ! Internals - integer(I4B) :: i + integer(I4B) :: i, npl real(DP) :: etaj, fac real(DP), dimension(NDIM) :: ah2, ah2o - associate(npl => pl%nbody) - ah2(:) = 0.0_DP - ah2o(:) = 0.0_DP - etaj = cb%Gmass - do concurrent(i = 2:npl, pl%lmask(i)) - etaj = etaj + pl%Gmass(i - 1) - fac = pl%Gmass(i) * cb%Gmass * pl%ir3j(i) / etaj - ah2(:) = ah2o + fac * pl%xj(:, i) - pl%ah(:,i) = pl%ah(:, i) + ah2(:) - ah2o(:) = ah2(:) - end do - end associate + npl = pl%nbody + ah2(:) = 0.0_DP + ah2o(:) = 0.0_DP + etaj = cb%Gmass + do concurrent(i = 2:npl, pl%lmask(i)) + etaj = etaj + pl%Gmass(i - 1) + fac = pl%Gmass(i) * cb%Gmass * pl%ir3j(i) / etaj + ah2(:) = ah2o + fac * pl%xj(:, i) + pl%ah(:,i) = pl%ah(:, i) + ah2(:) + ah2o(:) = ah2(:) + end do return end subroutine whm_kick_getacch_ah2 @@ -215,9 +215,10 @@ module subroutine whm_kick_vh_pl(self, nbody_system, param, t, dt, lbeg) real(DP), intent(in) :: dt !! Stepsize logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. ! Internals - integer(I4B) :: i + integer(I4B) :: i, npl - associate(pl => self, npl => self%nbody, cb => nbody_system%cb) + associate(pl => self, cb => nbody_system%cb) + npl = self%nbody if (npl == 0) return if (lbeg) then if (pl%lfirst) then @@ -257,11 +258,12 @@ module subroutine whm_kick_vh_tp(self, nbody_system, param, t, dt, lbeg) real(DP), intent(in) :: dt !! Stepsize logical, intent(in) :: lbeg !! Logical flag indicating whether this is the beginning of the half step or not. ! Internals - integer(I4B) :: i + integer(I4B) :: i, ntp if (self%nbody == 0) return - associate(tp => self, ntp => self%nbody) + associate(tp => self) + ntp = self%nbody if (tp%lfirst) then do concurrent(i = 1:ntp, tp%lmask(i)) tp%ah(:, i) = 0.0_DP