From 27ae7369614f412140fb02040973e9b112ccd8b2 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 9 Dec 2021 11:48:19 -0500 Subject: [PATCH] Added in override to gr_pos_kick for SyMBA, which allows the position kick operation to be masked during recursive steps so only bodies at the current recursion depth are kicked --- src/helio/helio_gr.f90 | 23 +++++++----- src/helio/helio_step.f90 | 12 +++--- src/modules/fraggle_classes.f90 | 2 +- src/modules/helio_classes.f90 | 24 ++++++------ src/modules/symba_classes.f90 | 23 ++++++++++++ src/modules/whm_classes.f90 | 18 +++++---- src/symba/symba_gr.f90 | 66 +++++++++++++++++++++++++++++++++ src/symba/symba_step.f90 | 17 ++++++--- src/whm/whm_gr.f90 | 22 ++++++----- src/whm/whm_step.f90 | 8 ++-- 10 files changed, 160 insertions(+), 55 deletions(-) create mode 100644 src/symba/symba_gr.f90 diff --git a/src/helio/helio_gr.f90 b/src/helio/helio_gr.f90 index 5d2936324..1ff8eb5d5 100644 --- a/src/helio/helio_gr.f90 +++ b/src/helio/helio_gr.f90 @@ -56,7 +56,7 @@ module pure subroutine helio_gr_kick_getacch_tp(self, param) end subroutine helio_gr_kick_getacch_tp - module pure subroutine helio_gr_p4_pl(self, param, dt) + module pure subroutine helio_gr_p4_pl(self, system, param, dt) !! author: David A. Minton !! !! Position kick to massive bodies due to p**4 term in the post-Newtonian correction @@ -65,11 +65,12 @@ module pure subroutine helio_gr_p4_pl(self, param, dt) !! Adapted from David A. Minton's Swifter routine routine gr_helio_p4.f90 implicit none ! Arguments - class(helio_pl), intent(inout) :: self !! Swiftest particle object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: dt !! Step size + class(helio_pl), intent(inout) :: self !! Swiftest particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size ! Internals - integer(I4B) :: i + integer(I4B) :: i if (self%nbody == 0) return @@ -82,7 +83,8 @@ module pure subroutine helio_gr_p4_pl(self, param, dt) return end subroutine helio_gr_p4_pl - module pure subroutine helio_gr_p4_tp(self, param, dt) + + module pure subroutine helio_gr_p4_tp(self, system, param, dt) !! author: David A. Minton !! !! Position kick to test particles due to p**4 term in the post-Newtonian correction @@ -91,11 +93,12 @@ module pure subroutine helio_gr_p4_tp(self, param, dt) !! Adapted from David A. Minton's Swifter routine routine gr_helio_p4.f90 implicit none ! Arguments - class(helio_tp), intent(inout) :: self !! Swiftest particle object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: dt !! Step size + class(helio_tp), intent(inout) :: self !! Swiftest particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size ! Internals - integer(I4B) :: i + integer(I4B) :: i if (self%nbody == 0) return diff --git a/src/helio/helio_step.f90 b/src/helio/helio_step.f90 index 039884596..fb14c9e75 100644 --- a/src/helio/helio_step.f90 +++ b/src/helio/helio_step.f90 @@ -35,7 +35,7 @@ module subroutine helio_step_pl(self, system, param, t, dt) implicit none ! Arguments class(helio_pl), intent(inout) :: self !! Helio massive body particle data structure - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nboody system + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time real(DP), intent(in) :: dt !! Stepsize @@ -54,10 +54,10 @@ module subroutine helio_step_pl(self, system, param, t, dt) end if call pl%lindrift(cb, dth, lbeg=.true.) call pl%kick(system, param, t, dth, lbeg=.true.) - if (param%lgr) call pl%gr_pos_kick(param, dth) + if (param%lgr) call pl%gr_pos_kick(system, param, dth) call pl%drift(system, param, dt) call pl%kick(system, param, t + dt, dth, lbeg=.false.) - if (param%lgr) call pl%gr_pos_kick(param, dth) + if (param%lgr) call pl%gr_pos_kick(system, param, dth) call pl%lindrift(cb, dth, lbeg=.false.) call pl%vb2vh(cb) end select @@ -78,7 +78,7 @@ module subroutine helio_step_tp(self, system, param, t, dt) implicit none ! Arguments class(helio_tp), intent(inout) :: self !! Helio test particle data structure - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nboody system + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time real(DP), intent(in) :: dt !! Stepsize @@ -97,10 +97,10 @@ module subroutine helio_step_tp(self, system, param, t, dt) end if call tp%lindrift(cb, dth, lbeg=.true.) call tp%kick(system, param, t, dth, lbeg=.true.) - if (param%lgr) call tp%gr_pos_kick(param, dth) + if (param%lgr) call tp%gr_pos_kick(system, param, dth) call tp%drift(system, param, dt) call tp%kick(system, param, t + dt, dth, lbeg=.false.) - if (param%lgr) call tp%gr_pos_kick(param, dth) + if (param%lgr) call tp%gr_pos_kick(system, param, dth) call tp%lindrift(cb, dth, lbeg=.false.) call tp%vb2vh(vbcb = -cb%ptend) end select diff --git a/src/modules/fraggle_classes.f90 b/src/modules/fraggle_classes.f90 index cd648c04b..042ff5a78 100644 --- a/src/modules/fraggle_classes.f90 +++ b/src/modules/fraggle_classes.f90 @@ -150,7 +150,7 @@ module subroutine fraggle_placeholder_step(self, system, param, t, dt) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none class(fraggle_fragments), intent(inout) :: self !! Helio massive body particle object - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nboody system + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time real(DP), intent(in) :: dt !! Stepsiz diff --git a/src/modules/helio_classes.f90 b/src/modules/helio_classes.f90 index 2d15565b2..78a4bdc34 100644 --- a/src/modules/helio_classes.f90 +++ b/src/modules/helio_classes.f90 @@ -122,20 +122,22 @@ module pure subroutine helio_gr_kick_getacch_tp(self, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine helio_gr_kick_getacch_tp - module pure subroutine helio_gr_p4_pl(self, param, dt) - use swiftest_classes, only : swiftest_parameters + module pure subroutine helio_gr_p4_pl(self, system, param, dt) + use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system implicit none - class(helio_pl), intent(inout) :: self !! Swiftest particle object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: dt !! Step size + class(helio_pl), intent(inout) :: self !! Swiftest particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size end subroutine helio_gr_p4_pl - module pure subroutine helio_gr_p4_tp(self, param, dt) - use swiftest_classes, only : swiftest_parameters + module pure subroutine helio_gr_p4_tp(self, system, param, dt) + use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system implicit none - class(helio_tp), intent(inout) :: self !! Swiftest particle object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: dt !! Step size + class(helio_tp), intent(inout) :: self !! Swiftest particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size end subroutine helio_gr_p4_tp module subroutine helio_kick_getacch_pl(self, system, param, t, lbeg) @@ -191,7 +193,7 @@ module subroutine helio_step_pl(self, system, param, t, dt) use swiftest_classes, only : swiftest_nbody_system, swiftest_parameters implicit none class(helio_pl), intent(inout) :: self !! Helio massive body particle object - class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nboody system + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system class(swiftest_parameters), intent(inout) :: param !! Current run configuration parameters real(DP), intent(in) :: t !! Current simulation time real(DP), intent(in) :: dt !! Stepsize diff --git a/src/modules/symba_classes.f90 b/src/modules/symba_classes.f90 index 4a584eda9..99c527143 100644 --- a/src/modules/symba_classes.f90 +++ b/src/modules/symba_classes.f90 @@ -75,6 +75,7 @@ module symba_classes procedure :: discard => symba_discard_pl !! Process massive body discards procedure :: drift => symba_drift_pl !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level procedure :: encounter_check => symba_encounter_check_pl !! Checks if massive bodies are going through close encounters with each other + procedure :: gr_pos_kick => symba_gr_p4_pl !! Position kick due to p**4 term in the post-Newtonian correction procedure :: accel_int => symba_kick_getacch_int_pl !! Compute direct cross (third) term heliocentric accelerations of massive bodiess, with no mutual interactions between bodies below GMTINY procedure :: accel => symba_kick_getacch_pl !! Compute heliocentric accelerations of massive bodies procedure :: setup => symba_setup_pl !! Constructor method - Allocates space for the input number of bodies @@ -113,6 +114,7 @@ module symba_classes contains procedure :: drift => symba_drift_tp !! Method for Danby drift in Democratic Heliocentric coordinates. Sets the mask to the current recursion level procedure :: encounter_check => symba_encounter_check_tp !! Checks if any test particles are undergoing a close encounter with a massive body + procedure :: gr_pos_kick => symba_gr_p4_tp !! Position kick due to p**4 term in the post-Newtonian correction procedure :: accel => symba_kick_getacch_tp !! Compute heliocentric accelerations of test particles procedure :: setup => symba_setup_tp !! Constructor method - Allocates space for the input number of bodies procedure :: append => symba_util_append_tp !! Appends elements from one structure to another @@ -272,6 +274,7 @@ module subroutine symba_drift_tp(self, system, param, dt) end subroutine symba_drift_tp module function symba_encounter_check_pl(self, param, system, dt, irec) result(lany_encounter) + use swiftest_classes, only : swiftest_nbody_system implicit none class(symba_pl), intent(inout) :: self !! SyMBA test particle object class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters @@ -282,6 +285,7 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l end function symba_encounter_check_pl module function symba_encounter_check(self, param, system, dt, irec) result(lany_encounter) + use swiftest_classes, only : swiftest_parameters implicit none class(symba_encounter), intent(inout) :: self !! SyMBA pl-pl encounter list object class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters @@ -292,6 +296,7 @@ module function symba_encounter_check(self, param, system, dt, irec) result(lany end function symba_encounter_check module function symba_encounter_check_tp(self, param, system, dt, irec) result(lany_encounter) + use swiftest_classes, only : swiftest_parameters implicit none class(symba_tp), intent(inout) :: self !! SyMBA test particle object class(swiftest_parameters), intent(inout) :: param !! Current swiftest run configuration parameters @@ -301,6 +306,24 @@ module function symba_encounter_check_tp(self, param, system, dt, irec) result(l logical :: lany_encounter !! Returns true if there is at least one close encounter end function symba_encounter_check_tp + module pure subroutine symba_gr_p4_pl(self, system, param, dt) + use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system + implicit none + class(symba_pl), intent(inout) :: self !! SyMBA massive body object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size + end subroutine symba_gr_p4_pl + + module pure subroutine symba_gr_p4_tp(self, system, param, dt) + use swiftest_classes, only : swiftest_parameters, swiftest_nbody_system + implicit none + class(symba_tp), intent(inout) :: self !! SyMBA test particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size + end subroutine symba_gr_p4_tp + module function symba_collision_casedisruption(system, param, colliders, frag) result(status) use fraggle_classes, only : fraggle_colliders, fraggle_fragments implicit none diff --git a/src/modules/whm_classes.f90 b/src/modules/whm_classes.f90 index 066975d41..a7cd2f49c 100644 --- a/src/modules/whm_classes.f90 +++ b/src/modules/whm_classes.f90 @@ -170,20 +170,22 @@ module pure subroutine whm_gr_kick_getacch_tp(self, param) class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters end subroutine whm_gr_kick_getacch_tp - module pure subroutine whm_gr_p4_pl(self, param, dt) + module pure subroutine whm_gr_p4_pl(self, system, param, dt) use swiftest_classes, only : swiftest_parameters implicit none - class(whm_pl), intent(inout) :: self !! WHM massive body object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: dt !! Step size + class(whm_pl), intent(inout) :: self !! WHM massive body object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size end subroutine whm_gr_p4_pl - module pure subroutine whm_gr_p4_tp(self, param, dt) + module pure subroutine whm_gr_p4_tp(self, system, param, dt) use swiftest_classes, only : swiftest_parameters implicit none - class(whm_tp), intent(inout) :: self !! WHM test particle object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: dt !! Step size + class(whm_tp), intent(inout) :: self !! WHM test particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size end subroutine whm_gr_p4_tp !> Reads WHM massive body object in from file diff --git a/src/symba/symba_gr.f90 b/src/symba/symba_gr.f90 new file mode 100644 index 000000000..1b0246aa8 --- /dev/null +++ b/src/symba/symba_gr.f90 @@ -0,0 +1,66 @@ +submodule(symba_classes) s_symba_gr + use swiftest +contains + + module pure subroutine symba_gr_p4_pl(self, system, param, dt) + !! author: David A. Minton + !! + !! Position kick to massive bodies due to p**4 term in the post-Newtonian correction + !! Based on Saha & Tremaine (1994) Eq. 28 + !! + !! Adapted from David A. Minton's Swifter routine routine gr_symba_p4.f90 + implicit none + ! Arguments + class(symba_pl), intent(inout) :: self !! SyMBA massive body object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size + ! Internals + integer(I4B) :: i + + if (self%nbody == 0) return + + associate(pl => self, npl => self%nbody) + select type(system) + class is (symba_nbody_system) + do concurrent(i = 1:npl, pl%lmask(i) .and. pl%levelg(i) == system%irec ) + call gr_p4_pos_kick(param, pl%xh(:, i), pl%vb(:, i), dt) + end do + end select + end associate + + return + end subroutine symba_gr_p4_pl + + + module pure subroutine symba_gr_p4_tp(self, system, param, dt) + !! author: David A. Minton + !! + !! Position kick to test particles due to p**4 term in the post-Newtonian correction + !! Based on Saha & Tremaine (1994) Eq. 28 + !! + !! Adapted from David A. Minton's Swifter routine routine gr_symba_p4.f90 + implicit none + ! Arguments + class(symba_tp), intent(inout) :: self !! SyMBA test particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system object + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size + ! Internals + integer(I4B) :: i + + if (self%nbody == 0) return + + associate(tp => self, ntp => self%nbody) + select type(system) + class is (symba_nbody_system) + do concurrent(i = 1:ntp, tp%lmask(i) .and. tp%levelg(i) == system%irec) + call gr_p4_pos_kick(param, tp%xh(:, i), tp%vb(:, i), dt) + end do + end select + end associate + + return + end subroutine symba_gr_p4_tp + +end submodule s_symba_gr diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 330f47a56..f10850e7f 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -71,24 +71,24 @@ module subroutine symba_step_interp_system(self, param, t, dt) call pl%vh2vb(cb) call pl%lindrift(cb, dth, lbeg=.true.) call pl%kick(system, param, t, dth, lbeg=.true.) - if (param%lgr) call pl%gr_pos_kick(param, dth) + if (param%lgr) call pl%gr_pos_kick(system, param, dth) call pl%drift(system, param, dt) call tp%vh2vb(vbcb = -cb%ptbeg) call tp%lindrift(cb, dth, lbeg=.true.) call tp%kick(system, param, t, dth, lbeg=.true.) - if (param%lgr) call tp%gr_pos_kick(param, dth) + if (param%lgr) call tp%gr_pos_kick(system, param, dth) call tp%drift(system, param, dt) call system%recursive_step(param, t, 0) call pl%kick(system, param, t, dth, lbeg=.false.) - if (param%lgr) call pl%gr_pos_kick(param, dth) + if (param%lgr) call pl%gr_pos_kick(system, param, dth) call pl%lindrift(cb, dth, lbeg=.false.) call pl%vb2vh(cb) call tp%kick(system, param, t, dth, lbeg=.false.) - if (param%lgr) call tp%gr_pos_kick(param, dth) + if (param%lgr) call tp%gr_pos_kick(system, param, dth) call tp%lindrift(cb, dth, lbeg=.false.) call tp%vb2vh(vbcb = -cb%ptend) end select @@ -185,7 +185,10 @@ module recursive subroutine symba_step_recur_system(self, param, t, ireci) call plplenc_list%kick(system, dth, irecp, -1) call pltpenc_list%kick(system, dth, irecp, -1) end if - + if (param%lgr) then + call pl%gr_pos_kick(system, param, dth) + call tp%gr_pos_kick(system, param, dth) + end if call pl%drift(system, param, dtl) call tp%drift(system, param, dtl) @@ -198,6 +201,10 @@ module recursive subroutine symba_step_recur_system(self, param, t, ireci) call plplenc_list%kick(system, dth, irecp, -1) call pltpenc_list%kick(system, dth, irecp, -1) end if + if (param%lgr) then + call pl%gr_pos_kick(system, param, dth) + call tp%gr_pos_kick(system, param, dth) + end if if (param%lclose) then lplpl_collision = plplenc_list%collision_check(system, param, t+dtl, dtl, ireci) diff --git a/src/whm/whm_gr.f90 b/src/whm/whm_gr.f90 index bbe6ca30b..12ae82a35 100644 --- a/src/whm/whm_gr.f90 +++ b/src/whm/whm_gr.f90 @@ -61,7 +61,7 @@ module pure subroutine whm_gr_kick_getacch_tp(self, param) end subroutine whm_gr_kick_getacch_tp - module pure subroutine whm_gr_p4_pl(self, param, dt) + module pure subroutine whm_gr_p4_pl(self, system, param, dt) !! author: David A. Minton !! !! Position kick to massive bodies due to p**4 term in the post-Newtonian correction @@ -70,11 +70,12 @@ module pure subroutine whm_gr_p4_pl(self, param, dt) !! Adapted from David A. Minton's Swifter routine routine gr_whm_p4.f90 implicit none ! Arguments - class(whm_pl), intent(inout) :: self !! Swiftest particle object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: dt !! Step size + class(whm_pl), intent(inout) :: self !! Swiftest particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size ! Internals - integer(I4B) :: i + integer(I4B) :: i associate(pl => self, npl => self%nbody) if (npl == 0) return @@ -87,7 +88,7 @@ module pure subroutine whm_gr_p4_pl(self, param, dt) end subroutine whm_gr_p4_pl - module pure subroutine whm_gr_p4_tp(self, param, dt) + module pure subroutine whm_gr_p4_tp(self, system, param, dt) !! author: David A. Minton !! !! Position kick to test particles due to p**4 term in the post-Newtonian correction @@ -96,11 +97,12 @@ module pure subroutine whm_gr_p4_tp(self, param, dt) !! Adapted from David A. Minton's Swifter routine routine gr_whm_p4.f90 implicit none ! Arguments - class(whm_tp), intent(inout) :: self !! Swiftest particle object - class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters - real(DP), intent(in) :: dt !! Step size + class(whm_tp), intent(inout) :: self !! Swiftest particle object + class(swiftest_nbody_system), intent(inout) :: system !! Swiftest nbody system + class(swiftest_parameters), intent(in) :: param !! Current run configuration parameters + real(DP), intent(in) :: dt !! Step size ! Internals - integer(I4B) :: i + integer(I4B) :: i associate(tp => self, ntp => self%nbody) if (ntp == 0) return diff --git a/src/whm/whm_step.f90 b/src/whm/whm_step.f90 index d194e2c02..244ebbbb1 100644 --- a/src/whm/whm_step.f90 +++ b/src/whm/whm_step.f90 @@ -50,9 +50,9 @@ module subroutine whm_step_pl(self, system, param, t, dt) dth = 0.5_DP * dt call pl%kick(system, param, t, dth,lbeg=.true.) call pl%vh2vj(cb) - if (param%lgr) call pl%gr_pos_kick(param, dth) + if (param%lgr) call pl%gr_pos_kick(system, param, dth) call pl%drift(system, param, dt) - if (param%lgr) call pl%gr_pos_kick(param, dth) + if (param%lgr) call pl%gr_pos_kick(system, param, dth) call pl%j2h(cb) call pl%kick(system, param, t + dt, dth, lbeg=.false.) end associate @@ -85,9 +85,9 @@ module subroutine whm_step_tp(self, system, param, t, dt) associate(tp => self, cb => system%cb, pl => system%pl) dth = 0.5_DP * dt call tp%kick(system, param, t, dth, lbeg=.true.) - if (param%lgr) call tp%gr_pos_kick(param, dth) + if (param%lgr) call tp%gr_pos_kick(system, param, dth) call tp%drift(system, param, dt) - if (param%lgr) call tp%gr_pos_kick(param, dth) + if (param%lgr) call tp%gr_pos_kick(system, param, dth) call tp%kick(system, param, t + dt, dth, lbeg=.false.) end associate end select