diff --git a/src/discard/discard.f90 b/src/discard/discard.f90 index 9c8044c61..2998e8804 100644 --- a/src/discard/discard.f90 +++ b/src/discard/discard.f90 @@ -53,7 +53,7 @@ module subroutine discard_pl(self, system, param) class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameter if (self%nbody == 0) return - self%ldiscard(:) = .false. + self%ldiscard(1:self%nbody) = .false. return end subroutine discard_pl @@ -86,9 +86,9 @@ module subroutine discard_tp(self, system, param) if ((param%rmin >= 0.0_DP) .or. (param%rmax >= 0.0_DP) .or. (param%rmaxu >= 0.0_DP)) call discard_cb_tp(tp, system, param) if (param%qmin >= 0.0_DP) call discard_peri_tp(tp, system, param) if (param%lclose) call discard_pl_tp(tp, system, param) - if (any(tp%ldiscard)) then + if (any(tp%ldiscard(1:ntp))) then allocate(ldiscard, source=tp%ldiscard) - call tp%spill(system%tp_discards, ldiscard, ldestructive=.true.) + call tp%spill(system%tp_discards, ldiscard(1:ntp), ldestructive=.true.) end if end associate diff --git a/src/helio/helio_kick.f90 b/src/helio/helio_kick.f90 index eebd17f53..8e7571a0b 100644 --- a/src/helio/helio_kick.f90 +++ b/src/helio/helio_kick.f90 @@ -65,9 +65,9 @@ module subroutine helio_kick_getacch_tp(self, system, param, t, lbeg) associate(tp => self, cb => system%cb, pl => system%pl, npl => system%pl%nbody) system%lbeg = lbeg if (system%lbeg) then - call tp%accel_int(pl%Gmass(:), pl%xbeg(:,:), npl) + call tp%accel_int(pl%Gmass(1:npl), pl%xbeg(:,1:npl), npl) else - call tp%accel_int(pl%Gmass(:), pl%xend(:,:), npl) + call tp%accel_int(pl%Gmass(1:npl), pl%xend(:,1:npl), npl) end if if (param%loblatecb) call tp%accel_obl(system) if (param%lextra_force) call tp%accel_user(system, param, t, lbeg) @@ -99,7 +99,7 @@ module subroutine helio_kick_vb_pl(self, system, param, t, dt, lbeg) if (self%nbody == 0) return associate(pl => self, npl => self%nbody) - pl%ah(:,:) = 0.0_DP + pl%ah(:, 1:npl) = 0.0_DP call pl%accel(system, param, t, lbeg) if (lbeg) then call pl%set_beg_end(xbeg = pl%xh) @@ -136,7 +136,7 @@ module subroutine helio_kick_vb_tp(self, system, param, t, dt, lbeg) if (self%nbody == 0) return associate(tp => self, ntp => self%nbody) - tp%ah(:,:) = 0.0_DP + tp%ah(:, 1:ntp) = 0.0_DP call tp%accel(system, param, t, lbeg) do concurrent(i = 1:ntp, tp%lmask(i)) tp%vb(:, i) = tp%vb(:, i) + tp%ah(:, i) * dt diff --git a/src/kick/kick.f90 b/src/kick/kick.f90 index c41bf7649..af704c13c 100644 --- a/src/kick/kick.f90 +++ b/src/kick/kick.f90 @@ -64,6 +64,8 @@ module pure subroutine kick_getacch_int_tp(self, GMpl, xhp, npl) real(DP) :: rji2, irij3, fac, r2 real(DP), dimension(NDIM) :: dx + if ((self%nbody == 0) .or. (npl == 0)) return + associate(tp => self, ntp => self%nbody) do concurrent(i = 1:ntp, tp%lmask(i)) do j = 1, npl diff --git a/src/rmvs/rmvs_encounter_check.f90 b/src/rmvs/rmvs_encounter_check.f90 index e4c441472..4cb119e48 100644 --- a/src/rmvs/rmvs_encounter_check.f90 +++ b/src/rmvs/rmvs_encounter_check.f90 @@ -23,13 +23,14 @@ module function rmvs_encounter_check_tp(self, system, dt) result(lencounter) real(DP), dimension(system%pl%nbody) :: r2crit logical :: lflag + lencounter = .false. if (self%nbody == 0) return select type(pl => system%pl) class is (rmvs_pl) associate(tp => self, ntp => self%nbody, npl => pl%nbody, rts => system%rts) - r2crit(:) = (rts * pl%rhill(:))**2 - tp%plencP(:) = 0 + r2crit(1:npl) = (rts * pl%rhill(1:npl))**2 + tp%plencP(1:ntp) = 0 do j = 1, npl do i = 1, ntp if ((.not.tp%lmask(i)).or.(tp%plencP(i) /= 0)) cycle @@ -41,9 +42,9 @@ module function rmvs_encounter_check_tp(self, system, dt) result(lencounter) lflag = rmvs_chk_ind(r2, v2, vdotr, dt, r2crit(j)) if (lflag) tp%plencP(i) = j end do - pl%nenc(j) = count(tp%plencP(:) == j) + pl%nenc(j) = count(tp%plencP(1:ntp) == j) end do - lencounter = any(pl%nenc(:) > 0) + lencounter = any(pl%nenc(1:npl) > 0) end associate end select return diff --git a/src/rmvs/rmvs_kick.f90 b/src/rmvs/rmvs_kick.f90 index 018ada8f3..5d9514763 100644 --- a/src/rmvs/rmvs_kick.f90 +++ b/src/rmvs/rmvs_kick.f90 @@ -61,7 +61,9 @@ module subroutine rmvs_kick_getacch_tp(self, system, param, t, lbeg) end if ! Swap the planetocentric and heliocentric position vectors and central body masses - tp%xh(:,:) = tp%xheliocentric(:,:) + do concurrent(i = 1:ntp, tp%lmask(i)) + tp%xh(:, i) = tp%xheliocentric(:, i) + end do GMcb_original = cb%Gmass cb%Gmass = tp%cb_heliocentric%Gmass @@ -71,7 +73,7 @@ module subroutine rmvs_kick_getacch_tp(self, system, param, t, lbeg) if (param%lgr) call tp%accel_gr(param) ! Put everything back the way we found it - tp%xh(:,:) = xh_original(:,:) + call move_alloc(xh_original, tp%xh) cb%Gmass = GMcb_original end associate diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index c801ee16f..b8cfcd688 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -21,6 +21,8 @@ module subroutine rmvs_step_system(self, param, t, dt) real(DP), dimension(:,:), allocatable :: xbeg, xend, vbeg integer(I4B) :: i + if (self%tp%nbody == 0) call whm_step_system(self, param, t, dt) + select type(cb => self%cb) class is (rmvs_cb) select type(pl => self%pl) @@ -36,11 +38,11 @@ module subroutine rmvs_step_system(self, param, t, dt) lencounter = tp%encounter_check(system, dt) if (lencounter) then lfirstpl = pl%lfirst - pl%outer(0)%x(:,:) = xbeg(:,:) - pl%outer(0)%v(:,:) = vbeg(:,:) + pl%outer(0)%x(:, 1:npl) = xbeg(:, 1:npl) + pl%outer(0)%v(:, 1:npl) = vbeg(:, 1:npl) call pl%step(system, param, t, dt) - pl%outer(NTENC)%x(:,:) = pl%xh(:,:) - pl%outer(NTENC)%v(:,:) = pl%vh(:,:) + pl%outer(NTENC)%x(:, 1:npl) = pl%xh(:, 1:npl) + pl%outer(NTENC)%v(:, 1:npl) = pl%vh(:, 1:npl) call rmvs_interp_out(cb, pl, dt) call rmvs_step_out(cb, pl, tp, system, param, t, dt) tp%lmask(1:ntp) = .not. tp%lmask(1:ntp) @@ -89,10 +91,10 @@ subroutine rmvs_interp_out(cb, pl, dt) allocate(GMcb(npl)) allocate(dto(npl)) allocate(iflag(npl)) - dto(:) = dt / dntenc - GMcb(:) = cb%Gmass - xtmp(:,:) = pl%outer(0)%x(:, :) - vtmp(:,:) = pl%outer(0)%v(:, :) + dto(1:npl) = dt / dntenc + GMcb(1:npl) = cb%Gmass + xtmp(:,1:npl) = pl%outer(0)%x(:, 1:npl) + vtmp(:,1:npl) = pl%outer(0)%v(:, 1:npl) do outer_index = 1, NTENC - 1 call drift_one(GMcb(1:npl), xtmp(1,1:npl), xtmp(2,1:npl), xtmp(3,1:npl), & vtmp(1,1:npl), vtmp(2,1:npl), vtmp(3,1:npl), & @@ -110,11 +112,11 @@ subroutine rmvs_interp_out(cb, pl, dt) end do end if frac = 1.0_DP - outer_index / dntenc - pl%outer(outer_index)%x(:, :) = frac * xtmp(:,:) - pl%outer(outer_index)%v(:, :) = frac * vtmp(:,:) + pl%outer(outer_index)%x(:, 1:npl) = frac * xtmp(:, 1:npl) + pl%outer(outer_index)%v(:, 1:npl) = frac * vtmp(:, 1:npl) end do - xtmp(:,:) = pl%outer(NTENC)%x(:, :) - vtmp(:,:) = pl%outer(NTENC)%v(:, :) + xtmp(:, 1:npl) = pl%outer(NTENC)%x(:, 1:npl) + vtmp(:, 1:npl) = pl%outer(NTENC)%v(:, 1:npl) do outer_index = NTENC - 1, 1, -1 call drift_one(GMcb(1:npl), xtmp(1,1:npl), xtmp(2,1:npl), xtmp(3,1:npl), & vtmp(1,1:npl), vtmp(2,1:npl), vtmp(3,1:npl), & @@ -132,8 +134,8 @@ subroutine rmvs_interp_out(cb, pl, dt) end do end if frac = outer_index / dntenc - pl%outer(outer_index)%x(:, :) = pl%outer(outer_index)%x(:, :) + frac * xtmp(:,:) - pl%outer(outer_index)%v(:, :) = pl%outer(outer_index)%v(:, :) + frac * vtmp(:,:) + pl%outer(outer_index)%x(:, 1:npl) = pl%outer(outer_index)%x(:, 1:npl) + frac * xtmp(:, 1:npl) + pl%outer(outer_index)%v(:, 1:npl) = pl%outer(outer_index)%v(:, 1:npl) + frac * vtmp(:, 1:npl) end do end associate @@ -165,16 +167,16 @@ subroutine rmvs_step_out(cb, pl, tp, system, param, t, dt) associate(npl => pl%nbody, ntp => tp%nbody) dto = dt / NTENC - where(tp%plencP(:) == 0) - tp%lmask(:) = .false. + where(tp%plencP(1:ntp) == 0) + tp%lmask(1:ntp) = .false. elsewhere - tp%lperi(:) = .false. + tp%lperi(1:ntp) = .false. end where do outer_index = 1, NTENC outer_time = t + (outer_index - 1) * dto - call pl%set_beg_end(xbeg = pl%outer(outer_index - 1)%x(:, :), & - vbeg = pl%outer(outer_index - 1)%v(:, :), & - xend = pl%outer(outer_index )%x(:, :)) + call pl%set_beg_end(xbeg = pl%outer(outer_index - 1)%x(:, 1:npl), & + vbeg = pl%outer(outer_index - 1)%v(:, 1:npl), & + xend = pl%outer(outer_index )%x(:, 1:npl)) system%rts = RHPSCALE lencounter = tp%encounter_check(system, dto) if (lencounter) then @@ -192,8 +194,8 @@ subroutine rmvs_step_out(cb, pl, tp, system, param, t, dt) do j = 1, npl if (pl%nenc(j) == 0) cycle tp%lfirst = .true. - where((tp%plencP(:) == j) .and. (.not.tp%lmask(:))) - tp%lmask(:) = .true. + where((tp%plencP(1:ntp) == j) .and. (.not.tp%lmask(1:ntp))) + tp%lmask(1:ntp) = .true. end where end do end do @@ -222,7 +224,7 @@ subroutine rmvs_interp_in(cb, pl, system, param, dt, outer_index) ! Internals integer(I4B) :: i, inner_index real(DP) :: frac, dntphenc - real(DP), dimension(:,:), allocatable :: xtmp, vtmp, xh_original + real(DP), dimension(:,:), allocatable :: xtmp, vtmp, xh_original, ah_original real(DP), dimension(:), allocatable :: GMcb, dti integer(I4B), dimension(:), allocatable :: iflag @@ -230,32 +232,33 @@ subroutine rmvs_interp_in(cb, pl, system, param, dt, outer_index) dntphenc = real(NTPHENC, kind=DP) ! Set the endpoints of the inner region from the outer region values in the current outer step index - pl%inner(0)%x(:,:) = pl%outer(outer_index - 1)%x(:, :) - pl%inner(0)%v(:,:) = pl%outer(outer_index - 1)%v(:, :) - pl%inner(NTPHENC)%x(:,:) = pl%outer(outer_index)%x(:, :) - pl%inner(NTPHENC)%v(:,:) = pl%outer(outer_index)%v(:, :) + pl%inner(0)%x(:, 1:npl) = pl%outer(outer_index - 1)%x(:, 1:npl) + pl%inner(0)%v(:, 1:npl) = pl%outer(outer_index - 1)%v(:, 1:npl) + pl%inner(NTPHENC)%x(:, 1:npl) = pl%outer(outer_index)%x(:, 1:npl) + pl%inner(NTPHENC)%v(:, 1:npl) = pl%outer(outer_index)%v(:, 1:npl) allocate(xtmp,mold=pl%xh) allocate(vtmp,mold=pl%vh) allocate(GMcb(npl)) allocate(dti(npl)) allocate(iflag(npl)) - dti(:) = dt / dntphenc - GMcb(:) = cb%Gmass - xtmp(:, :) = pl%inner(0)%x(:, :) - vtmp(:, :) = pl%inner(0)%v(:, :) + dti(1:npl) = dt / dntphenc + GMcb(1:npl) = cb%Gmass + xtmp(:, 1:npl) = pl%inner(0)%x(:, 1:npl) + vtmp(:, 1:npl) = pl%inner(0)%v(:, 1:npl) if ((param%loblatecb) .or. (param%ltides)) then allocate(xh_original, source=pl%xh) - pl%xh(:, :) = xtmp(:, :) ! Temporarily replace heliocentric position with inner substep values to calculate the oblateness terms + allocate(ah_original, source=pl%ah) + pl%xh(:, 1:npl) = xtmp(:, 1:npl) ! Temporarily replace heliocentric position with inner substep values to calculate the oblateness terms end if if (param%loblatecb) then call pl%accel_obl(system) - pl%inner(0)%aobl(:, :) = pl%aobl(:, :) ! Save the oblateness acceleration on the planet for this substep + pl%inner(0)%aobl(:, 1:npl) = pl%aobl(:, 1:npl) ! Save the oblateness acceleration on the planet for this substep end if if (param%ltides) then call pl%accel_tides(system) - pl%inner(0)%atide(:, :) = pl%atide(:, :) ! Save the oblateness acceleration on the planet for this substep + pl%inner(0)%atide(:, 1:npl) = pl%atide(:, 1:npl) ! Save the oblateness acceleration on the planet for this substep end if do inner_index = 1, NTPHENC - 1 @@ -275,12 +278,12 @@ subroutine rmvs_interp_in(cb, pl, system, param, dt, outer_index) end do end if frac = 1.0_DP - inner_index / dntphenc - pl%inner(inner_index)%x(:, :) = frac * xtmp(:,:) - pl%inner(inner_index)%v(:, :) = frac * vtmp(:,:) + pl%inner(inner_index)%x(:, 1:npl) = frac * xtmp(:, 1:npl) + pl%inner(inner_index)%v(:, 1:npl) = frac * vtmp(:, 1:npl) end do - xtmp(:,:) = pl%inner(NTPHENC)%x(:, :) - vtmp(:,:) = pl%inner(NTPHENC)%v(:, :) + xtmp(:, 1:npl) = pl%inner(NTPHENC)%x(:, 1:npl) + vtmp(:, 1:npl) = pl%inner(NTPHENC)%v(:, 1:npl) do inner_index = NTPHENC - 1, 1, -1 call drift_one(GMcb(1:npl), xtmp(1,1:npl), xtmp(2,1:npl), xtmp(3,1:npl), & @@ -299,31 +302,32 @@ subroutine rmvs_interp_in(cb, pl, system, param, dt, outer_index) end do end if frac = inner_index / dntphenc - pl%inner(inner_index)%x(:, :) = pl%inner(inner_index)%x(:, :) + frac * xtmp(:, :) - pl%inner(inner_index)%v(:, :) = pl%inner(inner_index)%v(:, :) + frac * vtmp(:, :) + pl%inner(inner_index)%x(:, 1:npl) = pl%inner(inner_index)%x(:, 1:npl) + frac * xtmp(:, 1:npl) + pl%inner(inner_index)%v(:, 1:npl) = pl%inner(inner_index)%v(:, 1:npl) + frac * vtmp(:, 1:npl) if (param%loblatecb) then - pl%xh(:,:) = pl%inner(inner_index)%x(:, :) + pl%xh(:,1:npl) = pl%inner(inner_index)%x(:, 1:npl) call pl%accel_obl(system) - pl%inner(inner_index)%aobl(:, :) = pl%aobl(:, :) + pl%inner(inner_index)%aobl(:, 1:npl) = pl%aobl(:, 1:npl) end if if (param%ltides) then call pl%accel_tides(system) - pl%inner(inner_index)%atide(:, :) = pl%atide(:, :) + pl%inner(inner_index)%atide(:, 1:npl) = pl%atide(:, 1:npl) end if end do if (param%loblatecb) then ! Calculate the final value of oblateness accelerations at the final inner substep - pl%xh(:,:) = pl%inner(NTPHENC)%x(:, :) + pl%xh(:, 1:npl) = pl%inner(NTPHENC)%x(:, 1:npl) call pl%accel_obl(system) - pl%inner(NTPHENC)%aobl(:, :) = pl%aobl(:, :) + pl%inner(NTPHENC)%aobl(:, 1:npl) = pl%aobl(:, 1:npl) end if if (param%ltides) then call pl%accel_tides(system) - pl%inner(NTPHENC)%atide(:, :) = pl%atide(:, :) + pl%inner(NTPHENC)%atide(:, 1:npl) = pl%atide(:, 1:npl) end if - ! Put the planet positions back into place + ! Put the planet positions and accelerations back into place if (allocated(xh_original)) call move_alloc(xh_original, pl%xh) + if (allocated(ah_original)) call move_alloc(ah_original, pl%ah) end associate return @@ -371,7 +375,7 @@ subroutine rmvs_step_in(cb, pl, tp, param, outer_time, dto) ! now step the encountering test particles fully through the inner encounter lfirsttp = .true. do inner_index = 1, NTPHENC ! Integrate over the encounter region, using the "substitute" planetocentric systems at each level - plenci%xh(:,:) = plenci%inner(inner_index - 1)%x(:,:) + plenci%xh(:, 1:npl) = plenci%inner(inner_index - 1)%x(:, 1:npl) call plenci%set_beg_end(xbeg = plenci%inner(inner_index - 1)%x, & xend = plenci%inner(inner_index)%x) @@ -428,7 +432,7 @@ subroutine rmvs_make_planetocentric(param, cb, pl, tp) ! There are inner encounters with this planet if (allocated(encmask)) deallocate(encmask) allocate(encmask(ntp)) - encmask(:) = tp%plencP(:) == i + encmask(1:ntp) = tp%plencP(1:ntp) == i allocate(rmvs_tp :: pl%planetocentric(i)%tp) ! Create encountering test particle structure select type(cbenci => pl%planetocentric(i)%cb) @@ -438,51 +442,53 @@ subroutine rmvs_make_planetocentric(param, cb, pl, tp) select type(tpenci => pl%planetocentric(i)%tp) class is (rmvs_tp) tpenci%lplanetocentric = .true. - call tpenci%setup(pl%nenc(i), param) - tpenci%cb_heliocentric = cb - tpenci%ipleP = i - tpenci%lmask(:) = .true. - tpenci%status(:) = ACTIVE - ! Grab all the encountering test particles and convert them to a planetocentric frame - tpenci%id(:) = pack(tp%id(:), encmask(:)) - do j = 1, NDIM - tpenci%xheliocentric(j, :) = pack(tp%xh(j,:), encmask(:)) - tpenci%xh(j, :) = tpenci%xheliocentric(j, :) - pl%inner(0)%x(j, i) - tpenci%vh(j, :) = pack(tp%vh(j,:), encmask(:)) - pl%inner(0)%v(j, i) - end do - tpenci%lperi(:) = pack(tp%lperi(:), encmask(:)) - tpenci%plperP(:) = pack(tp%plperP(:), encmask(:)) - ! Make sure that the test particles get the planetocentric value of mu - allocate(cbenci%inner(0:NTPHENC)) - do inner_index = 0, NTPHENC - allocate(plenci%inner(inner_index)%x, mold=pl%inner(inner_index)%x) - allocate(plenci%inner(inner_index)%v, mold=pl%inner(inner_index)%x) - allocate(cbenci%inner(inner_index)%x(NDIM,1)) - allocate(cbenci%inner(inner_index)%v(NDIM,1)) - cbenci%inner(inner_index)%x(:,1) = pl%inner(inner_index)%x(:, i) - cbenci%inner(inner_index)%v(:,1) = pl%inner(inner_index)%v(:, i) - plenci%inner(inner_index)%x(:,1) = -cbenci%inner(inner_index)%x(:,1) - plenci%inner(inner_index)%v(:,1) = -cbenci%inner(inner_index)%v(:,1) - - if (param%loblatecb) then - allocate(plenci%inner(inner_index)%aobl, mold=pl%inner(inner_index)%aobl) - allocate(cbenci%inner(inner_index)%aobl(NDIM,1)) - cbenci%inner(inner_index)%aobl(:,1) = pl%inner(inner_index)%aobl(:, i) - end if - - if (param%ltides) then - allocate(plenci%inner(inner_index)%atide, mold=pl%inner(inner_index)%atide) - allocate(cbenci%inner(inner_index)%atide(NDIM,1)) - cbenci%inner(inner_index)%atide(:,1) = pl%inner(inner_index)%atide(:, i) - end if - - do j = 2, npl - ipc2hc = plenci%plind(j) - plenci%inner(inner_index)%x(:,j) = pl%inner(inner_index)%x(:, ipc2hc) - cbenci%inner(inner_index)%x(:,1) - plenci%inner(inner_index)%v(:,j) = pl%inner(inner_index)%v(:, ipc2hc) - cbenci%inner(inner_index)%v(:,1) + associate(nenci => pl%nenc(i)) + call tpenci%setup(nenci, param) + tpenci%cb_heliocentric = cb + tpenci%ipleP = i + tpenci%lmask(1:nenci) = .true. + tpenci%status(1:nenci) = ACTIVE + ! Grab all the encountering test particles and convert them to a planetocentric frame + tpenci%id(1:nenci) = pack(tp%id(1:ntp), encmask(1:ntp)) + do j = 1, NDIM + tpenci%xheliocentric(j, 1:nenci) = pack(tp%xh(j,1:ntp), encmask(:)) + tpenci%xh(j, 1:nenci) = tpenci%xheliocentric(j, 1:nenci) - pl%inner(0)%x(j, i) + tpenci%vh(j, 1:nenci) = pack(tp%vh(j, 1:ntp), encmask(1:ntp)) - pl%inner(0)%v(j, i) + end do + tpenci%lperi(1:nenci) = pack(tp%lperi(1:ntp), encmask(1:ntp)) + tpenci%plperP(1:nenci) = pack(tp%plperP(1:ntp), encmask(1:ntp)) + ! Make sure that the test particles get the planetocentric value of mu + allocate(cbenci%inner(0:NTPHENC)) + do inner_index = 0, NTPHENC + allocate(plenci%inner(inner_index)%x, mold=pl%inner(inner_index)%x) + allocate(plenci%inner(inner_index)%v, mold=pl%inner(inner_index)%x) + allocate(cbenci%inner(inner_index)%x(NDIM,1)) + allocate(cbenci%inner(inner_index)%v(NDIM,1)) + cbenci%inner(inner_index)%x(:,1) = pl%inner(inner_index)%x(:, i) + cbenci%inner(inner_index)%v(:,1) = pl%inner(inner_index)%v(:, i) + plenci%inner(inner_index)%x(:,1) = -cbenci%inner(inner_index)%x(:,1) + plenci%inner(inner_index)%v(:,1) = -cbenci%inner(inner_index)%v(:,1) + + if (param%loblatecb) then + allocate(plenci%inner(inner_index)%aobl, mold=pl%inner(inner_index)%aobl) + allocate(cbenci%inner(inner_index)%aobl(NDIM,1)) + cbenci%inner(inner_index)%aobl(:,1) = pl%inner(inner_index)%aobl(:, i) + end if + + if (param%ltides) then + allocate(plenci%inner(inner_index)%atide, mold=pl%inner(inner_index)%atide) + allocate(cbenci%inner(inner_index)%atide(NDIM,1)) + cbenci%inner(inner_index)%atide(:,1) = pl%inner(inner_index)%atide(:, i) + end if + + do j = 2, npl + ipc2hc = plenci%plind(j) + plenci%inner(inner_index)%x(:,j) = pl%inner(inner_index)%x(:, ipc2hc) - cbenci%inner(inner_index)%x(:,1) + plenci%inner(inner_index)%v(:,j) = pl%inner(inner_index)%v(:, ipc2hc) - cbenci%inner(inner_index)%v(:,1) + end do end do - end do - call tpenci%set_mu(cbenci) + call tpenci%set_mu(cbenci) + end associate end select end select end select @@ -599,31 +605,33 @@ subroutine rmvs_end_planetocentric(pl, tp) class is (rmvs_pl) select type(tpenci => pl%planetocentric(i)%tp) class is (rmvs_tp) - if (allocated(tpind)) deallocate(tpind) - allocate(tpind(pl%nenc(i))) - ! Index array of encountering test particles - if (allocated(encmask)) deallocate(encmask) - allocate(encmask(ntp)) - encmask(:) = tp%plencP(:) == i - tpind(:) = pack([(j,j=1,ntp)], encmask(:)) - - ! Copy the results of the integration back over and shift back to heliocentric reference - tp%status(tpind(1:pl%nenc(i))) = tpenci%status(1:pl%nenc(i)) - tp%lmask(tpind(1:pl%nenc(i))) = tpenci%lmask(1:pl%nenc(i)) - do j = 1, NDIM - tp%xh(j, tpind(1:pl%nenc(i))) = tpenci%xh(j,1:pl%nenc(i)) + pl%inner(NTPHENC)%x(j, i) - tp%vh(j, tpind(1:pl%nenc(i))) = tpenci%vh(j,1:pl%nenc(i)) + pl%inner(NTPHENC)%v(j, i) - end do - tp%lperi(tpind(1:pl%nenc(i))) = tpenci%lperi(1:pl%nenc(i)) - tp%plperP(tpind(1:pl%nenc(i))) = tpenci%plperP(1:pl%nenc(i)) - deallocate(pl%planetocentric(i)%tp) - deallocate(cbenci%inner) - do inner_index = 0, NTPHENC - deallocate(plenci%inner(inner_index)%x) - deallocate(plenci%inner(inner_index)%v) - if (allocated(plenci%inner(inner_index)%aobl)) deallocate(plenci%inner(inner_index)%aobl) - if (allocated(plenci%inner(inner_index)%atide)) deallocate(plenci%inner(inner_index)%atide) - end do + associate(nenci => pl%nenc(i)) + if (allocated(tpind)) deallocate(tpind) + allocate(tpind(nenci)) + ! Index array of encountering test particles + if (allocated(encmask)) deallocate(encmask) + allocate(encmask(ntp)) + encmask(1:ntp) = tp%plencP(1:ntp) == i + tpind(1:nenci) = pack([(j,j=1,ntp)], encmask(1:ntp)) + + ! Copy the results of the integration back over and shift back to heliocentric reference + tp%status(tpind(1:nenci)) = tpenci%status(1:nenci) + tp%lmask(tpind(1:nenci)) = tpenci%lmask(1:nenci) + do j = 1, NDIM + tp%xh(j, tpind(1:nenci)) = tpenci%xh(j,1:nenci) + pl%inner(NTPHENC)%x(j, i) + tp%vh(j, tpind(1:nenci)) = tpenci%vh(j,1:nenci) + pl%inner(NTPHENC)%v(j, i) + end do + tp%lperi(tpind(1:nenci)) = tpenci%lperi(1:nenci) + tp%plperP(tpind(1:nenci)) = tpenci%plperP(1:nenci) + deallocate(pl%planetocentric(i)%tp) + deallocate(cbenci%inner) + do inner_index = 0, NTPHENC + deallocate(plenci%inner(inner_index)%x) + deallocate(plenci%inner(inner_index)%v) + if (allocated(plenci%inner(inner_index)%aobl)) deallocate(plenci%inner(inner_index)%aobl) + if (allocated(plenci%inner(inner_index)%atide)) deallocate(plenci%inner(inner_index)%atide) + end do + end associate end select end select end select diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 7c49009f2..b00a4a8e6 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -841,32 +841,32 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, ibiggest = maxloc(pl%Gmass(family(:)), dim=1) ! Copy over identification, information, and physical properties of the new bodies from the fragment list - plnew%id(:) = id_frag(:) + plnew%id(1:nfrag) = id_frag(1:nfrag) system%maxid = system%maxid + nfrag - plnew%xb(:,:) = xb_frag(:, :) - plnew%vb(:,:) = vb_frag(:, :) + plnew%xb(:, 1:nfrag) = xb_frag(:, 1:nfrag) + plnew%vb(:, 1:nfrag) = vb_frag(:, 1:nfrag) do i = 1, nfrag plnew%xh(:,i) = xb_frag(:, i) - cb%xb(:) plnew%vh(:,i) = vb_frag(:, i) - cb%vb(:) end do - plnew%mass(:) = m_frag(:) - plnew%Gmass(:) = param%GU * m_frag(:) - plnew%radius(:) = rad_frag(:) - plnew%density(:) = m_frag(:) / rad_frag(:) + plnew%mass(1:nfrag) = m_frag(1:nfrag) + plnew%Gmass(1:nfrag) = param%GU * m_frag(1:nfrag) + plnew%radius(1:nfrag) = rad_frag(1:nfrag) + plnew%density(1:nfrag) = m_frag(1:nfrag) / rad_frag(1:nfrag) select case(status) case(DISRUPTION) - plnew%info(:)%origin_type = "Disruption" - plnew%status(:) = NEW_PARTICLE - plnew%info(:)%origin_time = param%t + plnew%info(1:nfrag)%origin_type = "Disruption" + plnew%status(1:nfrag) = NEW_PARTICLE + plnew%info(1:nfrag)%origin_time = param%t do i = 1, nfrag plnew%info(i)%origin_xh(:) = plnew%xh(:,i) plnew%info(i)%origin_vh(:) = plnew%vh(:,i) end do case(SUPERCATASTROPHIC) - plnew%info(:)%origin_type = "Supercatastrophic" - plnew%status(:) = NEW_PARTICLE - plnew%info(:)%origin_time = param%t + plnew%info(1:nfrag)%origin_type = "Supercatastrophic" + plnew%status(1:nfrag) = NEW_PARTICLE + plnew%info(1:nfrag)%origin_time = param%t do i = 1, nfrag plnew%info(i)%origin_xh(:) = plnew%xh(:,i) plnew%info(i)%origin_vh(:) = plnew%vh(:,i) @@ -887,8 +887,8 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, end select if (param%lrotation) then - plnew%Ip(:,:) = Ip_frag(:,:) - plnew%rot(:,:) = rot_frag(:,:) + plnew%Ip(:, 1:nfrag) = Ip_frag(:, 1:nfrag) + plnew%rot(:, 1:nfrag) = rot_frag(:, 1:nfrag) end if if (param%ltides) then @@ -899,11 +899,11 @@ subroutine symba_collision_mergeaddsub(system, param, family, id_frag, Ip_frag, call plnew%set_mu(cb) !Copy over or set integration parameters for new bodies - plnew%lcollision(:) = .false. - plnew%ldiscard(:) = .false. - plnew%lmtiny(:) = plnew%Gmass(:) > param%GMTINY - plnew%levelg(:) = pl%levelg(ibiggest) - plnew%levelm(:) = pl%levelm(ibiggest) + plnew%lcollision(1:nfrag) = .false. + plnew%ldiscard(1:nfrag) = .false. + plnew%lmtiny(1:nfrag) = plnew%Gmass(1:nfrag) > param%GMTINY + plnew%levelg(1:nfrag) = pl%levelg(ibiggest) + plnew%levelm(1:nfrag) = pl%levelm(ibiggest) ! Append the new merged body to the list and record how many we made nstart = pl_adds%nbody + 1 diff --git a/src/symba/symba_drift.f90 b/src/symba/symba_drift.f90 index c4efee05f..e3bbc3660 100644 --- a/src/symba/symba_drift.f90 +++ b/src/symba/symba_drift.f90 @@ -14,13 +14,14 @@ module subroutine symba_drift_pl(self, system, param, dt) real(DP), intent(in) :: dt !! Stepsize if (self%nbody == 0) return - - select type(system) - class is (symba_nbody_system) - self%lmask(:) = self%status(:) /= INACTIVE .and. self%levelg(:) == system%irec - call helio_drift_body(self, system, param, dt) - self%lmask(:) = self%status(:) /= INACTIVE - end select + associate(pl => self, npl => self%nbody) + select type(system) + class is (symba_nbody_system) + pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE .and. pl%levelg(1:npl) == system%irec + call helio_drift_body(pl, system, param, dt) + pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE + end select + end associate return end subroutine symba_drift_pl @@ -38,13 +39,14 @@ module subroutine symba_drift_tp(self, system, param, dt) real(DP), intent(in) :: dt !! Stepsize if (self%nbody == 0) return - - select type(system) - class is (symba_nbody_system) - self%lmask(:) = self%status(:) /= INACTIVE .and. self%levelg(:) == system%irec - call helio_drift_body(self, system, param, dt) - self%lmask(:) = self%status(:) /= INACTIVE - end select + associate (tp => self, ntp => self%nbody) + select type(system) + class is (symba_nbody_system) + tp%lmask(1:ntp) = tp%status(1:ntp) /= INACTIVE .and. tp%levelg(1:ntp) == system%irec + call helio_drift_body(tp, system, param, dt) + tp%lmask(1:ntp) = tp%status(1:ntp) /= INACTIVE + end select + end associate return end subroutine symba_drift_tp diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index eb230b7e0..fe593bf7c 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -177,14 +177,14 @@ module function symba_encounter_check_tp(self, system, dt, irec) result(lany_enc call pltpenc_list%resize(nenc) pltpenc_list%status(1:nenc) = ACTIVE pltpenc_list%level(1:nenc) = irec - pltpenc_list%lvdotr(1:nenc) = pack(loc_lvdotr(:,:), lencounter(:,:)) - pltpenc_list%index1(1:nenc) = pack(spread([(i, i = 1, npl)], dim=1, ncopies=ntp), lencounter(:,:)) - pltpenc_list%index2(1:nenc) = pack(spread([(i, i = 1, ntp)], dim=2, ncopies=npl), lencounter(:,:)) + pltpenc_list%lvdotr(1:nenc) = pack(loc_lvdotr(1:ntp, 1:npl), lencounter(1:ntp, 1:npl)) + pltpenc_list%index1(1:nenc) = pack(spread([(i, i = 1, npl)], dim=1, ncopies=ntp), lencounter(1:ntp, 1:npl)) + pltpenc_list%index2(1:nenc) = pack(spread([(i, i = 1, ntp)], dim=2, ncopies=npl), lencounter(1:ntp, 1:npl)) pltpenc_list%id1(1:nenc) = pl%id(pltpenc_list%index1(1:nenc)) pltpenc_list%id2(1:nenc) = tp%id(pltpenc_list%index2(1:nenc)) select type(pl) class is (symba_pl) - pl%lencounter(:) = .false. + pl%lencounter(1:npl) = .false. do k = 1, nenc pl%lencounter(pltpenc_list%index1(k)) = .true. end do diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 775efa962..5ab525e23 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -273,34 +273,36 @@ module subroutine symba_io_read_particle(system, param) class is (symba_pl) select type(tp => system%tp) class is (symba_tp) - do - lmatch = .false. - read(LUN, err = 667, iomsg = errmsg) id - - if (idx == cb%id) then - read(LUN, err = 667, iomsg = errmsg) cb%info - lmatch = .true. - else - if (pl%nbody > 0) then - idx = findloc(pl%id(:), id, dim=1) - if (idx /= 0) then - read(LUN, err = 667, iomsg = errmsg) pl%info(idx) - lmatch = .true. + associate(npl => pl%nbody, ntp => tp%nbody) + do + lmatch = .false. + read(LUN, err = 667, iomsg = errmsg) id + + if (idx == cb%id) then + read(LUN, err = 667, iomsg = errmsg) cb%info + lmatch = .true. + else + if (npl > 0) then + idx = findloc(pl%id(1:npl), id, dim=1) + if (idx /= 0) then + read(LUN, err = 667, iomsg = errmsg) pl%info(idx) + lmatch = .true. + end if end if - end if - if (.not.lmatch .and. tp%nbody > 0) then - idx = findloc(tp%id(:), id, dim=1) - if (idx /= 0) then - read(LUN, err = 667, iomsg = errmsg) tp%info(idx) - lmatch = .true. + if (.not.lmatch .and. ntp > 0) then + idx = findloc(tp%id(1:ntp), id, dim=1) + if (idx /= 0) then + read(LUN, err = 667, iomsg = errmsg) tp%info(idx) + lmatch = .true. + end if end if end if - end if - if (.not.lmatch) then - write(*,*) 'Particle id ',id,' not found. Skipping' - read(LUN, err = 667, iomsg = errmsg) tmpinfo - end if - end do + if (.not.lmatch) then + write(*,*) 'Particle id ',id,' not found. Skipping' + read(LUN, err = 667, iomsg = errmsg) tmpinfo + end if + end do + end associate close(unit = LUN, err = 667, iomsg = errmsg) end select end select diff --git a/src/symba/symba_kick.f90 b/src/symba/symba_kick.f90 index ccf56e039..1eabf9b6d 100644 --- a/src/symba/symba_kick.f90 +++ b/src/symba/symba_kick.f90 @@ -126,9 +126,9 @@ module subroutine symba_kick_encounter(self, system, dt, irec, sgn) class is (symba_pl) select type(tp => system%tp) class is (symba_tp) - associate(ind1 => self%index1, ind2 => self%index2) - if (pl%nbody > 0) pl%lmask(:) = pl%status(:) /= INACTIVE - if (tp%nbody > 0) tp%lmask(:) = tp%status(:) /= INACTIVE + associate(ind1 => self%index1, ind2 => self%index2, npl => pl%nbody, ntp => tp%nbody, nenc => self%nenc) + if (npl > 0) pl%lmask(1:npl) = pl%status(1:npl) /= INACTIVE + if (ntp > 0) tp%lmask(1:ntp) = tp%status(1:ntp) /= INACTIVE irm1 = irec - 1 @@ -139,13 +139,13 @@ module subroutine symba_kick_encounter(self, system, dt, irec, sgn) end if if (isplpl) then - pl%ah(:,ind1(1:self%nenc)) = 0.0_DP - pl%ah(:,ind2(1:self%nenc)) = 0.0_DP + pl%ah(:,ind1(1:nenc)) = 0.0_DP + pl%ah(:,ind2(1:nenc)) = 0.0_DP else - tp%ah(:,ind2(1:self%nenc)) = 0.0_DP + tp%ah(:,ind2(1:nenc)) = 0.0_DP end if - do k = 1, self%nenc + do k = 1, nenc if (isplpl) then lgoodlevel = (pl%levelg(ind1(k)) >= irm1) .and. (pl%levelg(ind2(k)) >= irm1) else @@ -176,15 +176,15 @@ module subroutine symba_kick_encounter(self, system, dt, irec, sgn) faci = fac * pl%Gmass(ind1(k)) if (isplpl) then facj = fac * pl%Gmass(ind2(k)) - pl%ah(:,ind1(k)) = pl%ah(:,ind1(k)) + facj * dx(:) - pl%ah(:,ind2(k)) = pl%ah(:,ind2(k)) - faci * dx(:) + pl%ah(:, ind1(k)) = pl%ah(:, ind1(k)) + facj * dx(:) + pl%ah(:, ind2(k)) = pl%ah(:, ind2(k)) - faci * dx(:) else - tp%ah(:,ind2(k)) = tp%ah(:,ind2(k)) - faci * dx(:) + tp%ah(:, ind2(k)) = tp%ah(:, ind2(k)) - faci * dx(:) end if end if end do if (isplpl) then - do k = 1, self%nenc + do k = 1, nenc pl%vb(:,ind1(k)) = pl%vb(:,ind1(k)) + sgn * dt * pl%ah(:,ind1(k)) pl%vb(:,ind2(k)) = pl%vb(:,ind2(k)) + sgn * dt * pl%ah(:,ind2(k)) pl%ah(:,ind1(k)) = 0.0_DP diff --git a/src/whm/whm_kick.f90 b/src/whm/whm_kick.f90 index 2da00c332..f831ce15b 100644 --- a/src/whm/whm_kick.f90 +++ b/src/whm/whm_kick.f90 @@ -80,17 +80,17 @@ module subroutine whm_kick_getacch_tp(self, system, param, t, lbeg) system%lbeg = lbeg if (lbeg) then - ah0(:) = whm_kick_getacch_ah0(pl%Gmass(:), pl%xbeg(:,:), npl) + ah0(:) = whm_kick_getacch_ah0(pl%Gmass(1:npl), pl%xbeg(:, 1:npl), npl) do concurrent(i = 1:ntp, tp%lmask(i)) tp%ah(:, i) = tp%ah(:, i) + ah0(:) end do - call tp%accel_int(pl%Gmass(:), pl%xbeg(:,:), npl) + call tp%accel_int(pl%Gmass(1:npl), pl%xbeg(:, 1:npl), npl) else - ah0(:) = whm_kick_getacch_ah0(pl%Gmass(:), pl%xend(:,:), npl) + ah0(:) = whm_kick_getacch_ah0(pl%Gmass(1:npl), pl%xend(:, 1:npl), npl) do concurrent(i = 1:ntp, tp%lmask(i)) tp%ah(:, i) = tp%ah(:, i) + ah0(:) end do - call tp%accel_int(pl%Gmass(:), pl%xend(:,:), npl) + call tp%accel_int(pl%Gmass(1:npl), pl%xend(:, 1:npl), npl) end if if (param%loblatecb) call tp%accel_obl(system) @@ -213,13 +213,13 @@ module subroutine whm_kick_vh_pl(self, system, param, t, dt, lbeg) if (lbeg) then if (pl%lfirst) then call pl%h2j(cb) - pl%ah(:,:) = 0.0_DP + pl%ah(:, 1:npl) = 0.0_DP call pl%accel(system, param, t, lbeg) pl%lfirst = .false. end if call pl%set_beg_end(xbeg = pl%xh) else - pl%ah(:,:) = 0.0_DP + pl%ah(:, 1:npl) = 0.0_DP call pl%accel(system, param, t, lbeg) call pl%set_beg_end(xend = pl%xh) end if @@ -254,20 +254,16 @@ module subroutine whm_kick_vh_tp(self, system, param, t, dt, lbeg) associate(tp => self, ntp => self%nbody) if (tp%lfirst) then - where(tp%lmask(1:ntp)) - tp%ah(1,1:ntp) = 0.0_DP - tp%ah(2,1:ntp) = 0.0_DP - tp%ah(3,1:ntp) = 0.0_DP - end where + do concurrent(i = 1:ntp, tp%lmask(i)) + tp%ah(:, i) = 0.0_DP + end do call tp%accel(system, param, t, lbeg=.true.) tp%lfirst = .false. end if if (.not.lbeg) then - where(tp%lmask(1:ntp)) - tp%ah(1,1:ntp) = 0.0_DP - tp%ah(2,1:ntp) = 0.0_DP - tp%ah(3,1:ntp) = 0.0_DP - end where + do concurrent(i = 1:ntp, tp%lmask(i)) + tp%ah(:, i) = 0.0_DP + end do call tp%accel(system, param, t, lbeg) end if do concurrent(i = 1:ntp, tp%lmask(i))