From 098748428ef19233b21952440b647c0d24e1fa52 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Tue, 23 Nov 2021 15:46:45 -0500 Subject: [PATCH] Shortened and/or split lines to fit gfortran's default line width --- src/discard/discard.f90 | 30 ++++++--- src/encounter/encounter_check.f90 | 94 ++++++++++++++++++----------- src/encounter/encounter_io.f90 | 10 +-- src/fraggle/fraggle_generate.f90 | 24 +++++--- src/fraggle/fraggle_regime.f90 | 11 ++-- src/fraggle/fraggle_set.f90 | 3 +- src/fraggle/fraggle_util.f90 | 6 +- src/gr/gr.f90 | 3 +- src/io/io.f90 | 29 +++++---- src/main/swiftest_driver.f90 | 6 +- src/modules/encounter_classes.f90 | 6 +- src/modules/swiftest_classes.f90 | 3 +- src/netcdf/netcdf.f90 | 7 ++- src/rmvs/rmvs_discard.f90 | 7 ++- src/rmvs/rmvs_encounter_check.f90 | 3 +- src/rmvs/rmvs_step.f90 | 9 ++- src/setup/setup.f90 | 9 ++- src/symba/symba_collision.f90 | 51 +++++++++++----- src/symba/symba_discard.f90 | 45 +++++++++----- src/symba/symba_encounter_check.f90 | 6 +- src/symba/symba_io.f90 | 21 ++++--- src/symba/symba_setup.f90 | 6 +- src/symba/symba_step.f90 | 28 ++++++--- src/symba/symba_util.f90 | 6 +- src/util/util_peri.f90 | 6 +- src/util/util_set.f90 | 3 +- src/walltime/walltime.f90 | 9 ++- 27 files changed, 282 insertions(+), 159 deletions(-) diff --git a/src/discard/discard.f90 b/src/discard/discard.f90 index ad5b73426..b7527c330 100644 --- a/src/discard/discard.f90 +++ b/src/discard/discard.f90 @@ -128,18 +128,22 @@ subroutine discard_cb_tp(tp, system, param) tp%status(i) = DISCARDED_RMAX write(idstr, *) tp%id(i) write(timestr, *) param%t - write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too far from the central body at t = " // trim(adjustl(timestr)) + write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & + " too far from the central body at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. tp%lmask(i) = .false. - call tp%info(i)%set_value(status="DISCARDED_RMAX", discard_time=param%t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i)) + call tp%info(i)%set_value(status="DISCARDED_RMAX", discard_time=param%t, discard_xh=tp%xh(:,i), & + discard_vh=tp%vh(:,i)) else if ((param%rmin >= 0.0_DP) .and. (rh2 < rmin2)) then tp%status(i) = DISCARDED_RMIN write(idstr, *) tp%id(i) write(timestr, *) param%t - write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too close to the central body at t = " // trim(adjustl(timestr)) + write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & + " too close to the central body at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. tp%lmask(i) = .false. - call tp%info(i)%set_value(status="DISCARDED_RMIN", discard_time=param%t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i), discard_body_id=cb%id) + call tp%info(i)%set_value(status="DISCARDED_RMIN", discard_time=param%t, discard_xh=tp%xh(:,i), & + discard_vh=tp%vh(:,i), discard_body_id=cb%id) else if (param%rmaxu >= 0.0_DP) then rb2 = dot_product(tp%xb(:, i), tp%xb(:, i)) vb2 = dot_product(tp%vb(:, i), tp%vb(:, i)) @@ -148,10 +152,12 @@ subroutine discard_cb_tp(tp, system, param) tp%status(i) = DISCARDED_RMAXU write(idstr, *) tp%id(i) write(timestr, *) param%t - write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) + write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & + " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. tp%lmask(i) = .false. - call tp%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=param%t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i)) + call tp%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=param%t, discard_xh=tp%xh(:,i), & + discard_vh=tp%vh(:,i)) end if end if end if @@ -199,9 +205,11 @@ subroutine discard_peri_tp(tp, system, param) tp%status(i) = DISCARDED_PERI write(idstr, *) tp%id(i) write(timestr, *) param%t - write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " perihelion distance too small at t = " // trim(adjustl(timestr)) + write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & + " perihelion distance too small at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. - call tp%info(i)%set_value(status="DISCARDED_PERI", discard_time=param%t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i), discard_body_id=pl%id(j)) + call tp%info(i)%set_value(status="DISCARDED_PERI", discard_time=param%t, discard_xh=tp%xh(:,i), & + discard_vh=tp%vh(:,i), discard_body_id=pl%id(j)) end if end if end if @@ -247,10 +255,12 @@ subroutine discard_pl_tp(tp, system, param) write(idstrj, *) pl%id(j) write(timestr, *) param%t write(*, *) "Test particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & - // " too close to massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstrj)) & + // " too close to massive body " // trim(adjustl(pl%info(i)%name)) // & + " (" // trim(adjustl(idstrj)) & // " at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. - call tp%info(i)%set_value(status="DISCARDED_PLR", discard_time=param%t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i), discard_body_id=pl%id(j)) + call tp%info(i)%set_value(status="DISCARDED_PLR", discard_time=param%t, discard_xh=tp%xh(:,i), & + discard_vh=tp%vh(:,i), discard_body_id=pl%id(j)) exit end if end do diff --git a/src/encounter/encounter_check.f90 b/src/encounter/encounter_check.f90 index a669855ca..ccd3d23ce 100644 --- a/src/encounter/encounter_check.f90 +++ b/src/encounter/encounter_check.f90 @@ -2,7 +2,8 @@ use swiftest contains - module subroutine encounter_check_all(nenc, index1, index2, x1, v1, x2, v2, renc1, renc2, dt, lencounter, lvdotr) + module subroutine encounter_check_all(nenc, index1, index2, x1, v1, x2, v2, renc1, renc2, dt, & + lencounter, lvdotr) !! author: David A. Minton !! !! Check for encounters between n pairs of bodies. @@ -46,7 +47,8 @@ module subroutine encounter_check_all(nenc, index1, index2, x1, v1, x2, v2, renc end subroutine encounter_check_all - module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, & + nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies. Choose between the standard triangular or the Sort & Sweep method based on user inputs @@ -68,11 +70,14 @@ module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, nenc, ind logical, save :: lfirst = .true. logical, save :: skipit = .false. ! This will be used to ensure that the sort & sweep subroutine gets called at least once before timing it so that the extent array is nearly sorted when it is timed integer(I8B) :: nplpl = 0_I8B + integer(I8B) :: k + type(walltimer), save :: timer if (param%ladaptive_encounters_plpl .and. (.not. skipit)) then nplpl = (npl * (npl - 1) / 2) if (nplpl > 0) then if (lfirst) then + call timer%reset() write(itimer%loopname, *) "encounter_check_all_plpl" write(itimer%looptype, *) "ENCOUNTER_PLPL" lfirst = .false. @@ -84,12 +89,15 @@ module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, nenc, ind param%lencounter_sas_plpl = .false. end if end if - + + call timer%start() if (param%lencounter_sas_plpl) then call encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) else call encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1, index2, lvdotr) end if + call timer%stop() + call timer%report(nsubsteps=1, message="plpl Encounter check :") if (skipit) then skipit = .false. @@ -102,11 +110,19 @@ module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, nenc, ind end if end if + ! TEMP FOR TESTING + open(unit=22,file="enclist.csv", status="replace") + do k = 1_I8B, nenc + write(22,*) index1(k), index2(k) + end do + close(22) + return end subroutine encounter_check_all_plpl - module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, & + nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between fully interacting massive bodies partially interacting massive bodies. Choose between the standard triangular or the Sort & Sweep method based on user inputs @@ -172,12 +188,14 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt, call encounter_check_all_plpl(tmp_param, nplm, xplm, vplm, rencm, dt, nenc, index1, index2, lvdotr) if (param%lencounter_sas_plpl) then - call encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, plmplt_nenc, plmplt_index1, plmplt_index2, plmplt_lvdotr) + call encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, & + plmplt_nenc, plmplt_index1, plmplt_index2, plmplt_lvdotr) else - call encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, plmplt_nenc, plmplt_index1, plmplt_index2, plmplt_lvdotr) + call encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, & + plmplt_nenc, plmplt_index1, plmplt_index2, plmplt_lvdotr) end if call timer%stop() - call timer%report(nsubsteps=1, message="Encounter check :") + call timer%report(nsubsteps=1, message="plplm Encounter check :") if (skipit) then skipit = .false. @@ -222,7 +240,8 @@ module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt, end subroutine encounter_check_all_plplm - module subroutine encounter_check_all_pltp(param, npl, ntp, xpl, vpl, xtp, vtp, renc, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_all_pltp(param, npl, ntp, xpl, vpl, xtp, vtp, renc, dt, & + nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies and test particles. Choose between the standard triangular or the Sort & Sweep method based on user inputs @@ -308,7 +327,7 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, nenc, in type(encounter_bounding_box), save :: boundingbox logical, dimension(:), allocatable :: lencounter integer(I2B), dimension(npl) :: vshift_min, vshift_max - type(walltimer) :: timer + type(walltimer), save :: timer if (npl == 0) return @@ -341,7 +360,8 @@ subroutine encounter_check_all_sort_and_sweep_plpl(npl, x, v, renc, dt, nenc, in end subroutine encounter_check_all_sort_and_sweep_plpl - subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, nenc, index1, index2, lvdotr) + subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, & + nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies and test particles. @@ -415,7 +435,8 @@ subroutine encounter_check_all_sort_and_sweep_plplm(nplm, nplt, xplm, vplm, xplt end subroutine encounter_check_all_sort_and_sweep_plplm - subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, nenc, index1, index2, lvdotr) + subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, & + nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies and test particles. @@ -488,7 +509,8 @@ subroutine encounter_check_all_sort_and_sweep_pltp(npl, ntp, xpl, vpl, xtp, vtp, end subroutine encounter_check_all_sort_and_sweep_pltp - pure subroutine encounter_check_all_sweep_one(i, n, xi, yi, zi, vxi, vyi, vzi, x, y, z, vx, vy, vz, renci, renc, dt, ind_arr, lgood, nenci, index1, index2, lvdotr) + pure subroutine encounter_check_all_sweep_one(i, n, xi, yi, zi, vxi, vyi, vzi, x, y, z, vx, vy, vz, renci, renc, dt, & + ind_arr, lgood, nenci, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between the ith body and all other bodies. @@ -539,7 +561,8 @@ pure subroutine encounter_check_all_sweep_one(i, n, xi, yi, zi, vxi, vyi, vzi, x end subroutine encounter_check_all_sweep_one - pure subroutine encounter_check_all_triangular_one(i, n, xi, yi, zi, vxi, vyi, vzi, x, y, z, vx, vy, vz, renci, renc, dt, ind_arr, lenci) + pure subroutine encounter_check_all_triangular_one(i, n, xi, yi, zi, vxi, vyi, vzi, x, y, z, vx, vy, vz, renci, renc, & + dt, ind_arr, lenci) !! author: David A. Minton !! !! Check for encounters between the ith body and all other bodies. @@ -607,7 +630,7 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1 logical, dimension(npl) :: lencounteri, lvdotri integer(I4B), dimension(:), allocatable, save :: ind_arr type(encounter_list), dimension(npl) :: lenc - type(walltimer) :: timer + type(walltimer), save :: timer call util_index_array(ind_arr, npl) @@ -630,7 +653,8 @@ subroutine encounter_check_all_triangular_plpl(npl, x, v, renc, dt, nenc, index1 end subroutine encounter_check_all_triangular_plpl - subroutine encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, nenc, index1, index2, lvdotr) + subroutine encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, & + nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies and test particles. Split off from the main subroutine for performance @@ -655,7 +679,7 @@ subroutine encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vp logical, dimension(nplt) :: lencounteri, lvdotri integer(I4B), dimension(:), allocatable, save :: ind_arr type(encounter_list), dimension(nplm) :: lenc - type(walltimer) :: timer + type(walltimer), save :: timer call util_index_array(ind_arr, nplt) @@ -678,7 +702,8 @@ subroutine encounter_check_all_triangular_plplm(nplm, nplt, xplm, vplm, xplt, vp end subroutine encounter_check_all_triangular_plplm - subroutine encounter_check_all_triangular_pltp(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, nenc, index1, index2, lvdotr) + subroutine encounter_check_all_triangular_pltp(npl, ntp, xpl, vpl, xtp, vtp, renc, dt, & + nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Check for encounters between massive bodies and test particles. Split off from the main subroutine for performance @@ -937,7 +962,8 @@ module pure subroutine encounter_check_sort_aabb_1D(self, n, extent_arr) end subroutine encounter_check_sort_aabb_1D - module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, x1, v1, x2, v2, renc1, renc2, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, x1, v1, x2, v2, renc1, renc2, dt, & + nenc, index1, index2, lvdotr) !! author: David A. Minton !! !! Sweeps the sorted bounding box extents and returns the true encounters (combines broad and narrow phases) @@ -967,8 +993,7 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, x1, v1, x type(encounter_list), dimension(n1+n2) :: lenc !! Array of encounter lists (one encounter list per body) integer(I4B), dimension(:), allocatable, save :: ind_arr integer(I8B) :: ibegix, iendix - integer(I8B), pointer :: ibegx, iendx - type(walltimer) :: timer1, timer2, timer3, timer4, timer5 + type(walltimer), save :: timer1, timer2, timer3, timer4, timer5 logical, save :: lfirst = .true. if (lfirst) then @@ -983,7 +1008,7 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, x1, v1, x ntot = n1 + n2 call util_index_array(ind_arr, ntot) - associate(ext_ind => self%aabb(1)%ind(:), ibegx => self%aabb(1)%ibeg(:), iendx => self%aabb(1)%iend(:)) + associate(ext_ind => self%aabb(1)%ind(:)) ! Sweep the intervals for each of the massive bodies along one dimension ! This will build a ragged pair of index lists inside of the lenc data structure allocate(ext_ind_true, mold=ext_ind) @@ -995,17 +1020,17 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, x1, v1, x endwhere llist1(:) = ext_ind_true(:) <= n1 - loverlap(:) = (iendx(:) + 1_I8B) > (ibegx(:) - 1_I8B) + loverlap(:) = (self%aabb(1)%ibeg(:) + 1_I8B) > (self%aabb(1)%iend(:) - 1_I8B) where(.not.loverlap(:)) lenc(:)%nenc = 0 call timer1%start() !$omp parallel do default(private) schedule(static)& - !$omp shared(ext_ind_true, ind_arr, ibegx, iendx, lenc, loverlap, x1, v1, x2, v2, renc1, renc2, llist1) & + !$omp shared(self, ext_ind_true, ind_arr, lenc, loverlap, x1, v1, x2, v2, renc1, renc2, llist1) & !$omp firstprivate(ntot, n1, n2, dt) do i = 1, n1 if (loverlap(i)) then - ibegix = ibegx(i) + 1_I8B - iendix = iendx(i) - 1_I8B + ibegix = self%aabb(1)%ibeg(i) + 1_I8B + iendix = self%aabb(1)%iend(i) - 1_I8B lencounter2(:) = .false. where(.not.llist1(ibegix:iendix)) lencounter2(ext_ind_true(ibegix:iendix) - n1) = .true. @@ -1020,12 +1045,12 @@ module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, x1, v1, x call timer2%start() !$omp parallel do default(private) schedule(static)& - !$omp shared(ext_ind_true, ind_arr, ibegx, iendx, lenc, loverlap, x1, v1, x2, v2, renc1, renc2, llist1) & + !$omp shared(self, ext_ind_true, ind_arr, lenc, loverlap, x1, v1, x2, v2, renc1, renc2, llist1) & !$omp firstprivate(ntot, n1, dt) do i = n1+1, ntot if (loverlap(i)) then - ibegix = ibegx(i) + 1_I8B - iendix = iendx(i) - 1_I8B + ibegix = self%aabb(1)%ibeg(i) + 1_I8B + iendix = self%aabb(1)%iend(i) - 1_I8B lencounter1(:) = .false. where(llist1(ibegix:iendix)) lencounter1(ext_ind_true(ibegix:iendix)) = .true. ii = i - n1 @@ -1079,8 +1104,7 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt type(encounter_list), dimension(n) :: lenc !! Array of encounter lists (one encounter list per body) integer(I4B), dimension(:), allocatable, save :: ind_arr integer(I8B) :: ibegix, iendix - integer(I8B), pointer :: ibegx, iendx - type(walltimer) :: timer1 + type(walltimer), save :: timer1 logical, save :: lfirst = .true. if (lfirst) then @@ -1090,7 +1114,7 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt call util_index_array(ind_arr, n) - associate(ext_ind => self%aabb(1)%ind(:), ibegx => self%aabb(1)%ibeg(:), iendx => self%aabb(1)%iend(:)) + associate(ext_ind => self%aabb(1)%ind(:)) ! Sweep the intervals for each of the massive bodies along one dimension ! This will build a ragged pair of index lists inside of the lenc data structure allocate(ext_ind_true, mold=ext_ind) @@ -1100,17 +1124,17 @@ module subroutine encounter_check_sweep_aabb_single_list(self, n, x, v, renc, dt ext_ind_true(:) = ext_ind(:) endwhere - loverlap(:) = (iendx(:) + 1_I8B) > (ibegx(:) - 1_I8B) + loverlap(:) = (self%aabb(1)%ibeg(:) + 1_I8B) > (self%aabb(1)%iend(:) - 1_I8B) where(.not.loverlap(:)) lenc(:)%nenc = 0 call timer1%start() !$omp parallel do default(private) schedule(static)& - !$omp shared(ext_ind_true, ind_arr, ibegx, iendx, lenc, loverlap, x, v, renc) & + !$omp shared(self, ext_ind_true, ind_arr, lenc, loverlap, x, v, renc) & !$omp firstprivate(n, dt) do i = 1, n if (loverlap(i)) then - ibegix = ibegx(i) + 1_I8B - iendix = iendx(i) - 1_I8B + ibegix = self%aabb(1)%ibeg(i) + 1_I8B + iendix = self%aabb(1)%iend(i) - 1_I8B lencounteri(:) = .false. lencounteri(ext_ind_true(ibegix:iendix)) = .true. diff --git a/src/encounter/encounter_io.f90 b/src/encounter/encounter_io.f90 index 60fdf5304..d7785a71a 100644 --- a/src/encounter/encounter_io.f90 +++ b/src/encounter/encounter_io.f90 @@ -22,9 +22,9 @@ module subroutine encounter_io_write_frame(iu, t, id1, id2, Gmass1, Gmass2, radi ! Internals character(len=STRMAX) :: errmsg - write(iu, err = 667, iomsg = errmsg) t - write(iu, err = 667, iomsg = errmsg) id1, xh1(1), xh1(2), xh1(3), vh1(1), vh1(2), Gmass1, radius1 - write(iu, err = 667, iomsg = errmsg) id2, xh2(1), xh2(2), xh2(3), vh2(1), vh2(2), Gmass2, radius2 + write(iu, err=667, iomsg=errmsg) t + write(iu, err=667, iomsg=errmsg) id1, xh1(1), xh1(2), xh1(3), vh1(1), vh1(2), Gmass1, radius1 + write(iu, err=667, iomsg=errmsg) id2, xh2(1), xh2(2), xh2(3), vh2(1), vh2(2), Gmass2, radius2 return 667 continue @@ -46,10 +46,10 @@ module subroutine encounter_io_write_list(self, pl, encbody, param) if (param%enc_out == "" .or. self%nenc == 0) return - open(unit = LUN, file = param%enc_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', iostat = ierr, iomsg = errmsg) + open(unit=LUN, file=param%enc_out, status='OLD', position='APPEND', form='UNFORMATTED', iostat=ierr, iomsg=errmsg) if (ierr /= 0) then if (lfirst) then - open(unit = LUN, file = param%enc_out, status = 'NEW', form = 'UNFORMATTED', err = 667, iomsg = errmsg) + open(unit=LUN, file=param%enc_out, status='NEW', form='UNFORMATTED', err=667, iomsg=errmsg) else goto 667 end if diff --git a/src/fraggle/fraggle_generate.f90 b/src/fraggle/fraggle_generate.f90 index 0b0126702..2b443bebe 100644 --- a/src/fraggle/fraggle_generate.f90 +++ b/src/fraggle/fraggle_generate.f90 @@ -110,14 +110,16 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa lfailure = ((abs(dEtot + frag%Qloss) > FRAGGLE_ETOL) .or. (dEtot > 0.0_DP)) if (lfailure) then write(message, *) dEtot, abs(dEtot + frag%Qloss) / FRAGGLE_ETOL - call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed due to high energy error: " // trim(adjustl(message))) + call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed due to high energy error: " // & + trim(adjustl(message))) cycle end if lfailure = ((abs(dLmag) / (.mag.frag%Ltot_before)) > FRAGGLE_LTOL) if (lfailure) then write(message,*) dLmag / (.mag.frag%Ltot_before(:)) - call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed due to high angular momentum error: " // trim(adjustl(message))) + call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle failed due to high angular momentum error: " // & + trim(adjustl(message))) cycle end if @@ -131,9 +133,11 @@ module subroutine fraggle_generate_fragments(self, colliders, system, param, lfa write(message,*) try if (lfailure) then - call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle fragment generation failed after " // trim(adjustl(message)) // " tries") + call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle fragment generation failed after " // & + trim(adjustl(message)) // " tries") else - call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle fragment generation succeeded after " // trim(adjustl(message)) // " tries") + call io_log_one_message(FRAGGLE_LOG_OUT, "Fraggle fragment generation succeeded after " // & + trim(adjustl(message)) // " tries") call fraggle_io_log_generate(frag) end if @@ -240,14 +244,16 @@ subroutine fraggle_generate_spins(frag, colliders, f_spin, lfailure) frag%ke_spin = 0.0_DP do i = 1, nfrag ! Convert a fraction (f_spin) of either the remaining angular momentum or kinetic energy budget into spin, whichever gives the smaller rotation so as not to blow any budgets - rot_ke(:) = sqrt(2 * f_spin * frag%ke_budget / (nfrag * frag%mass(i) * frag%radius(i)**2 * frag%Ip(3, i))) * L_remainder(:) / norm2(L_remainder(:)) + rot_ke(:) = sqrt(2 * f_spin * frag%ke_budget / (nfrag * frag%mass(i) * frag%radius(i)**2 * frag%Ip(3, i))) & + * L_remainder(:) / norm2(L_remainder(:)) rot_L(:) = f_spin * L_remainder(:) / (nfrag * frag%mass(i) * frag%radius(i)**2 * frag%Ip(3, i)) if (norm2(rot_ke) < norm2(rot_L)) then frag%rot(:,i) = rot_ke(:) else frag%rot(:, i) = rot_L(:) end if - frag%ke_spin = frag%ke_spin + frag%mass(i) * frag%Ip(3, i) * frag%radius(i)**2 * dot_product(frag%rot(:, i), frag%rot(:, i)) + frag%ke_spin = frag%ke_spin + frag%mass(i) * frag%Ip(3, i) * frag%radius(i)**2 & + * dot_product(frag%rot(:, i), frag%rot(:, i)) end do frag%ke_spin = 0.5_DP * frag%ke_spin @@ -340,7 +346,8 @@ subroutine fraggle_generate_tan_vel(frag, colliders, lfailure) frag%v_t_mag(1:nfrag) = solve_fragment_tan_vel(v_t_mag_input=v_t_initial(7:nfrag), lfailure=lfailure) ! Perform one final shift of the radial velocity vectors to align with the center of mass of the collisional system (the origin) - frag%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(frag%v_r_mag(1:nfrag), frag%v_r_unit(:,1:nfrag), frag%v_t_mag(1:nfrag), frag%v_t_unit(:,1:nfrag), frag%mass(1:nfrag), frag%vbcom(:)) + frag%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(frag%v_r_mag(1:nfrag), frag%v_r_unit(:,1:nfrag), frag%v_t_mag(1:nfrag), & + frag%v_t_unit(:,1:nfrag), frag%mass(1:nfrag), frag%vbcom(:)) do concurrent (i = 1:nfrag) frag%v_coll(:,i) = frag%vb(:,i) - frag%vbcom(:) end do @@ -506,7 +513,8 @@ subroutine fraggle_generate_rad_vel(frag, colliders, lfailure) ! Shift the radial velocity vectors to align with the center of mass of the collisional system (the origin) frag%ke_orbit = 0.0_DP - frag%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(frag%v_r_mag(1:nfrag), frag%v_r_unit(:,1:nfrag), frag%v_t_mag(1:nfrag), frag%v_t_unit(:,1:nfrag), frag%mass(1:nfrag), frag%vbcom(:)) + frag%vb(:,1:nfrag) = fraggle_util_vmag_to_vb(frag%v_r_mag(1:nfrag), frag%v_r_unit(:,1:nfrag), & + frag%v_t_mag(1:nfrag), frag%v_t_unit(:,1:nfrag), frag%mass(1:nfrag), frag%vbcom(:)) do i = 1, nfrag frag%v_coll(:, i) = frag%vb(:, i) - frag%vbcom(:) frag%ke_orbit = frag%ke_orbit + frag%mass(i) * dot_product(frag%vb(:, i), frag%vb(:, i)) diff --git a/src/fraggle/fraggle_regime.f90 b/src/fraggle/fraggle_regime.f90 index f53888df5..abe5766f0 100644 --- a/src/fraggle/fraggle_regime.f90 +++ b/src/fraggle/fraggle_regime.f90 @@ -50,8 +50,9 @@ module subroutine fraggle_regime_colliders(self, frag, system, param) dentot = sum(mass_si(:) * density_si(:)) / mtot !! Use the positions and velocities of the parents from indside the step (at collision) to calculate the collisional regime - call fraggle_regime_collresolve(Mcb_si, mass_si(jtarg), mass_si(jproj), radius_si(jtarg), radius_si(jproj), x1_si(:), x2_si(:),& - v1_si(:), v2_si(:), density_si(jtarg), density_si(jproj), min_mfrag_si, frag%regime, mlr, mslr, frag%Qloss) + call fraggle_regime_collresolve(Mcb_si, mass_si(jtarg), mass_si(jproj), radius_si(jtarg), radius_si(jproj), & + x1_si(:), x2_si(:), v1_si(:), v2_si(:), density_si(jtarg), density_si(jproj), & + min_mfrag_si, frag%regime, mlr, mslr, frag%Qloss) frag%mass_dist(1) = min(max(mlr, 0.0_DP), mtot) frag%mass_dist(2) = min(max(mslr, 0.0_DP), mtot) @@ -73,7 +74,8 @@ module subroutine fraggle_regime_colliders(self, frag, system, param) end subroutine fraggle_regime_colliders - subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, min_mfrag, regime, Mlr, Mslr, Qloss) + subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb2, den1, den2, min_mfrag, & + regime, Mlr, Mslr, Qloss) !! Author: Jennifer L.L. Pouplin, Carlisle A. Wishard, and David A. Minton !! !! Determine the collisional regime of two colliding bodies. @@ -182,7 +184,8 @@ subroutine fraggle_regime_collresolve(Mcb, m1, m2, rad1, rad2, xh1, xh2, vb1, vb Mlr = Mtot Mslr = 0.0_DP Qloss = 0.0_DP - call io_log_one_message(FRAGGLE_LOG_OUT, "Fragments would have mass below the minimum. Converting this collision into a merger.") + call io_log_one_message(FRAGGLE_LOG_OUT, & + "Fragments would have mass below the minimum. Converting this collision into a merger.") else if( Vimp < Vescp) then regime = COLLRESOLVE_REGIME_MERGE !perfect merging regime diff --git a/src/fraggle/fraggle_set.f90 b/src/fraggle/fraggle_set.f90 index 4012eda34..6902bfb48 100644 --- a/src/fraggle/fraggle_set.f90 +++ b/src/fraggle/fraggle_set.f90 @@ -212,7 +212,8 @@ module subroutine fraggle_set_natural_scale_factors(self, colliders) associate(frag => self) ! Set scale factors - frag%Escale = 0.5_DP * (colliders%mass(1) * dot_product(colliders%vb(:,1), colliders%vb(:,1)) + colliders%mass(2) * dot_product(colliders%vb(:,2), colliders%vb(:,2))) + frag%Escale = 0.5_DP * (colliders%mass(1) * dot_product(colliders%vb(:,1), colliders%vb(:,1)) & + + colliders%mass(2) * dot_product(colliders%vb(:,2), colliders%vb(:,2))) frag%dscale = sum(colliders%radius(:)) frag%mscale = frag%mtot frag%vscale = sqrt(frag%Escale / frag%mscale) diff --git a/src/fraggle/fraggle_util.f90 b/src/fraggle/fraggle_util.f90 index f65bd81c5..f2c082242 100644 --- a/src/fraggle/fraggle_util.f90 +++ b/src/fraggle/fraggle_util.f90 @@ -166,7 +166,8 @@ module subroutine fraggle_util_get_energy_momentum(self, colliders, system, para lexclude(1:npl_after) = .false. lexclude(colliders%idx(1:colliders%ncoll)) = .true. if (.not.allocated(tmpsys)) then - write(*,*) "Error in fraggle_util_get_energy_momentum. This must be called with lbefore=.true. at least once before calling it with lbefore=.false." + write(*,*) "Error in fraggle_util_get_energy_momentum. " // & + " This must be called with lbefore=.true. at least once before calling it with lbefore=.false." call util_exit(FAILURE) end if call fraggle_util_add_fragments_to_system(frag, colliders, tmpsys, tmpparam) @@ -232,7 +233,8 @@ module subroutine fraggle_util_restructure(self, colliders, try, f_spin, r_max_s ! Linearly interpolate the last two failed solution ke deficits to find a new distance value to try ke_tot_deficit = ke_tot_deficit - (frag%ke_budget - frag%ke_orbit - frag%ke_spin) ke_avg_deficit = ke_tot_deficit / try - delta_r = (r_max_start - r_max_start_old) * (ke_avg_deficit_target - ke_avg_deficit_old) / (ke_avg_deficit - ke_avg_deficit_old) + delta_r = (r_max_start - r_max_start_old) * (ke_avg_deficit_target - ke_avg_deficit_old) & + / (ke_avg_deficit - ke_avg_deficit_old) if (abs(delta_r) > delta_r_max) delta_r = sign(delta_r_max, delta_r) end if r_max_start_old = r_max_start diff --git a/src/gr/gr.f90 b/src/gr/gr.f90 index 8b6bb1120..b05465e98 100644 --- a/src/gr/gr.f90 +++ b/src/gr/gr.f90 @@ -28,7 +28,8 @@ module pure subroutine gr_kick_getaccb_ns_body(self, system, param) rmag = norm2(self%xh(:,i)) vmag2 = dot_product(self%vh(:,i), self%vh(:,i)) rdotv = dot_product(self%xh(:,i), self%vh(:,i)) - self%agr(:, i) = self%mu * inv_c2 / rmag**3 * ((4 * self%mu(i) / rmag - vmag2) * self%xh(:,i) + 4 * rdotv * self%vh(:,i)) + self%agr(:, i) = self%mu * inv_c2 / rmag**3 * ((4 * self%mu(i) / rmag - vmag2) & + * self%xh(:,i) + 4 * rdotv * self%vh(:,i)) end do select type(self) diff --git a/src/io/io.f90 b/src/io/io.f90 index 7ec71d9ae..18a22ccb9 100644 --- a/src/io/io.f90 +++ b/src/io/io.f90 @@ -29,10 +29,11 @@ module subroutine io_conservation_report(self, param, lterminal) associate(system => self, pl => self%pl, cb => self%cb, npl => self%pl%nbody) if (((param%out_type == REAL4_TYPE) .or. (param%out_type == REAL8_TYPE)) .and. (param%energy_out /= "")) then if (param%lfirstenergy .and. (param%out_stat /= "OLD")) then - open(unit = EGYIU, file = param%energy_out, form = "formatted", status = "replace", action = "write", err = 667, iomsg = errmsg) - write(EGYIU,EGYHEADER, err = 667, iomsg = errmsg) + open(unit=EGYIU, file=param%energy_out, form="formatted", status="replace", action="write", err=667, iomsg=errmsg) + write(EGYIU,EGYHEADER, err=667, iomsg=errmsg) else - open(unit = EGYIU, file = param%energy_out, form = "formatted", status = "old", action = "write", position = "append", err = 667, iomsg = errmsg) + open(unit=EGYIU, file=param%energy_out, form="formatted", status="old", action="write", & + position="append", err=667, iomsg=errmsg) end if end if @@ -177,9 +178,9 @@ module subroutine io_dump_particle_info_base(self, param, idx) if (lfirst) then select case(param%out_stat) case('APPEND') - open(unit = LUN, file = param%particle_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', err = 667, iomsg = errmsg) + open(unit=LUN, file=param%particle_out, status='OLD', position='APPEND', form='UNFORMATTED', err=667, iomsg=errmsg) case('NEW', 'UNKNOWN', 'REPLACE') - open(unit = LUN, file = param%particle_out, status = param%out_stat, form = 'UNFORMATTED', err = 667, iomsg = errmsg) + open(unit=LUN, file=param%particle_out, status=param%out_stat, form='UNFORMATTED', err=667, iomsg=errmsg) case default write(*,*) 'Invalid status code',trim(adjustl(param%out_stat)) call util_exit(FAILURE) @@ -187,7 +188,7 @@ module subroutine io_dump_particle_info_base(self, param, idx) lfirst = .false. else - open(unit = LUN, file = param%particle_out, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', err = 667, iomsg = errmsg) + open(unit=LUN, file=param%particle_out, status='OLD', position= 'APPEND', form='UNFORMATTED', err=667, iomsg=errmsg) end if select type(self) @@ -760,7 +761,8 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) iostat = -1 return end if - if ((param%in_type /= REAL8_TYPE) .and. (param%in_type /= "ASCII") .and. (param%in_type /= NETCDF_FLOAT_TYPE) .and. (param%in_type /= NETCDF_DOUBLE_TYPE)) then + if ((param%in_type /= REAL8_TYPE) .and. (param%in_type /= "ASCII") & + .and. (param%in_type /= NETCDF_FLOAT_TYPE) .and. (param%in_type /= NETCDF_DOUBLE_TYPE)) then write(iomsg,*) 'Invalid input file type:',trim(adjustl(param%in_type)) iostat = -1 return @@ -788,7 +790,8 @@ module subroutine io_param_reader(self, unit, iotype, v_list, iostat, iomsg) iostat = -1 return end if - if ((param%out_stat /= "NEW") .and. (param%out_stat /= "REPLACE") .and. (param%out_stat /= "APPEND") .and. (param%out_stat /= "UNKNOWN")) then + if ((param%out_stat /= "NEW") .and. (param%out_stat /= "REPLACE") .and. (param%out_stat /= "APPEND") & + .and. (param%out_stat /= "UNKNOWN")) then write(iomsg,*) 'Invalid out_stat: ',trim(adjustl(param%out_stat)) iostat = -1 return @@ -1956,9 +1959,9 @@ module subroutine io_write_discard(self, param) end if select case(out_stat) case('APPEND') - open(unit = LUN, file = param%discard_out, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) + open(unit=LUN, file=param%discard_out, status='OLD', position='APPEND', form='FORMATTED', err=667, iomsg=errmsg) case('NEW', 'REPLACE', 'UNKNOWN') - open(unit = LUN, file = param%discard_out, status = param%out_stat, form = 'FORMATTED', err = 667, iomsg = errmsg) + open(unit=LUN, file=param%discard_out, status=param%out_stat, form='FORMATTED', err=667, iomsg=errmsg) case default write(*,*) 'Invalid status code for OUT_STAT: ',trim(adjustl(param%out_stat)) call util_exit(FAILURE) @@ -2137,9 +2140,9 @@ module subroutine io_write_frame_system(self, param) if (lfirst) then select case(param%out_stat) case('APPEND') - open(unit = iu, file = param%outfile, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', err = 667, iomsg = errmsg) + open(unit=iu, file=param%outfile, status='OLD', position='APPEND', form='UNFORMATTED', err=667, iomsg=errmsg) case('NEW', 'REPLACE', 'UNKNOWN') - open(unit = iu, file = param%outfile, status = param%out_stat, form = 'UNFORMATTED', err = 667, iomsg = errmsg) + open(unit=iu, file=param%outfile, status=param%out_stat, form='UNFORMATTED', err=667, iomsg=errmsg) case default write(*,*) 'Invalid status code for OUT_STAT: ',trim(adjustl(param%out_stat)) call util_exit(FAILURE) @@ -2147,7 +2150,7 @@ module subroutine io_write_frame_system(self, param) lfirst = .false. else - open(unit = iu, file = param%outfile, status = 'OLD', position = 'APPEND', form = 'UNFORMATTED', err = 667, iomsg = errmsg) + open(unit=iu, file=param%outfile, status='OLD', position= 'APPEND', form='UNFORMATTED', err=667, iomsg=errmsg) end if else if ((param%out_type == NETCDF_FLOAT_TYPE) .or. (param%out_type == NETCDF_DOUBLE_TYPE)) then diff --git a/src/main/swiftest_driver.f90 b/src/main/swiftest_driver.f90 index d4e1f2dc8..1bf4dba8b 100644 --- a/src/main/swiftest_driver.f90 +++ b/src/main/swiftest_driver.f90 @@ -22,8 +22,10 @@ program swiftest_driver real(DP) :: old_t_final = 0.0_DP !! Output time at which writing should start, in order to prevent duplicate lines being written for restarts type(walltimer) :: integration_timer !! Object used for computing elapsed wall time real(DP) :: tfrac - character(*), parameter :: statusfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, "; Number of active pl, tp = ", I5, ", ", I5)' - character(*), parameter :: symbastatfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, "; Number of active plm, pl, tp = ", I5, ", ", I5, ", ", I5)' + character(*), parameter :: statusfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // & + '"; Number of active pl, tp = ", I5, ", ", I5)' + character(*), parameter :: symbastatfmt = '("Time = ", ES12.5, "; fraction done = ", F6.3, ' // & + '"; Number of active plm, pl, tp = ", I5, ", ", I5, ", ", I5)' ierr = io_get_args(integrator, param_file_name) if (ierr /= 0) then diff --git a/src/modules/encounter_classes.f90 b/src/modules/encounter_classes.f90 index 5c4e963d6..841037bcb 100644 --- a/src/modules/encounter_classes.f90 +++ b/src/modules/encounter_classes.f90 @@ -83,7 +83,8 @@ module subroutine encounter_check_all_plpl(param, npl, x, v, renc, dt, nenc, ind integer(I8B), intent(out) :: nenc !! Total number of encounters end subroutine encounter_check_all_plpl - module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_all_plplm(param, nplm, nplt, xplm, vplm, xplt, vplt, rencm, renct, dt, & + nenc, index1, index2, lvdotr) use swiftest_classes, only: swiftest_parameters implicit none class(swiftest_parameters), intent(inout) :: param !! Current Swiftest run configuration parameter5s @@ -148,7 +149,8 @@ module pure subroutine encounter_check_sort_aabb_1D(self, n, extent_arr) real(DP), dimension(:), intent(in) :: extent_arr !! Array of extents of size 2*n end subroutine encounter_check_sort_aabb_1D - module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, x1, v1, x2, v2, renc1, renc2, dt, nenc, index1, index2, lvdotr) + module subroutine encounter_check_sweep_aabb_double_list(self, n1, n2, x1, v1, x2, v2, renc1, renc2, dt, & + nenc, index1, index2, lvdotr) implicit none class(encounter_bounding_box), intent(inout) :: self !! Multi-dimensional bounding box structure integer(I4B), intent(in) :: n1 !! Number of bodies 1 diff --git a/src/modules/swiftest_classes.f90 b/src/modules/swiftest_classes.f90 index f93fcf517..987375461 100644 --- a/src/modules/swiftest_classes.f90 +++ b/src/modules/swiftest_classes.f90 @@ -1574,7 +1574,8 @@ module subroutine util_set_mu_tp(self, cb) class(swiftest_cb), intent(inout) :: cb !! Swiftest central body object end subroutine util_set_mu_tp - module subroutine util_set_particle_info(self, name, particle_type, status, origin_type, origin_time, collision_id, origin_xh, origin_vh, discard_time, discard_xh, discard_vh, discard_body_id) + module subroutine util_set_particle_info(self, name, particle_type, status, origin_type, origin_time, collision_id, & + origin_xh, origin_vh, discard_time, discard_xh, discard_vh, discard_body_id) implicit none class(swiftest_particle_info), intent(inout) :: self character(len=*), intent(in), optional :: name !! Non-unique name diff --git a/src/netcdf/netcdf.f90 b/src/netcdf/netcdf.f90 index d28fcc788..7acc79b02 100644 --- a/src/netcdf/netcdf.f90 +++ b/src/netcdf/netcdf.f90 @@ -250,7 +250,8 @@ module subroutine netcdf_initialize_output(self, param) call check( nf90_def_var(self%ncid, ORIGIN_TIME_VARNAME, self%out_type, self%id_dimid, self%origin_time_varid) ) !call check( nf90_def_var_chunking(self%ncid, self%origin_time_varid, NF90_CHUNKED, [self%id_chunk]) ) - call check( nf90_def_var(self%ncid, ORIGIN_TYPE_VARNAME, NF90_CHAR, [self%str_dimid, self%id_dimid], self%origin_type_varid) ) + call check( nf90_def_var(self%ncid, ORIGIN_TYPE_VARNAME, NF90_CHAR, [self%str_dimid, self%id_dimid], & + self%origin_type_varid) ) !call check( nf90_def_var_chunking(self%ncid, self%origin_type_varid, NF90_CHUNKED, [NAMELEN, self%id_chunk]) ) call check( nf90_def_var(self%ncid, ORIGIN_XHX_VARNAME, self%out_type, self%id_dimid, self%origin_xhx_varid) ) !call check( nf90_def_var_chunking(self%ncid, self%origin_xhx_varid, NF90_CHUNKED, [self%id_chunk]) ) @@ -936,7 +937,9 @@ module subroutine netcdf_write_frame_base(self, iu, param) select type(self) class is (swiftest_pl) ! Additional output if the passed polymorphic object is a massive body call check( nf90_put_var(iu%ncid, iu%Gmass_varid, self%Gmass(j), start=[idslot, tslot]) ) - if (param%lrhill_present) call check( nf90_put_var(iu%ncid, iu%rhill_varid, self%rhill(j), start=[idslot, tslot]) ) + if (param%lrhill_present) then + call check( nf90_put_var(iu%ncid, iu%rhill_varid, self%rhill(j), start=[idslot, tslot]) ) + end if if (param%lclose) call check( nf90_put_var(iu%ncid, iu%radius_varid, self%radius(j), start=[idslot, tslot]) ) if (param%lrotation) then call check( nf90_put_var(iu%ncid, iu%Ip1_varid, self%Ip(1, j), start=[idslot, tslot]) ) diff --git a/src/rmvs/rmvs_discard.f90 b/src/rmvs/rmvs_discard.f90 index b179178dc..8e3b1d7f0 100644 --- a/src/rmvs/rmvs_discard.f90 +++ b/src/rmvs/rmvs_discard.f90 @@ -30,11 +30,12 @@ module subroutine rmvs_discard_tp(self, system, param) write(idstrj, *) pl%id(iplperP) write(timestr, *) t write(*, *) "Particle " // trim(adjustl(tp%info(i)%name)) // " (" // trim(adjustl(idstri)) & - // ") q with respect to massive body " // trim(adjustl(pl%info(iplperP)%name)) // " (" // trim(adjustl(idstrj)) & - // ") is too small at t = " // trim(adjustl(timestr)) + // ") q with respect to massive body " // trim(adjustl(pl%info(iplperP)%name)) & + // " (" // trim(adjustl(idstrj)) // ") is too small at t = " // trim(adjustl(timestr)) tp%ldiscard(i) = .true. tp%lmask(i) = .false. - call tp%info(i)%set_value(status="DISCARDED_PLQ", discard_time=t, discard_xh=tp%xh(:,i), discard_vh=tp%vh(:,i), discard_body_id=pl%id(iplperP)) + call tp%info(i)%set_value(status="DISCARDED_PLQ", discard_time=t, discard_xh=tp%xh(:,i), & + discard_vh=tp%vh(:,i), discard_body_id=pl%id(iplperP)) end if end if end associate diff --git a/src/rmvs/rmvs_encounter_check.f90 b/src/rmvs/rmvs_encounter_check.f90 index f698ea4f0..bbb5a4382 100644 --- a/src/rmvs/rmvs_encounter_check.f90 +++ b/src/rmvs/rmvs_encounter_check.f90 @@ -36,7 +36,8 @@ module function rmvs_encounter_check_tp(self, param, system, dt) result(lencount class is (rmvs_pl) associate(tp => self, ntp => self%nbody, npl => pl%nbody) tp%plencP(1:ntp) = 0 - call encounter_check_all_pltp(param, npl, ntp, pl%xbeg, pl%vbeg, tp%xh, tp%vh, pl%renc, dt, nenc, index1, index2, lvdotr) + call encounter_check_all_pltp(param, npl, ntp, pl%xbeg, pl%vbeg, tp%xh, tp%vh, pl%renc, dt, & + nenc, index1, index2, lvdotr) lencounter = (nenc > 0_I8B) if (lencounter) then diff --git a/src/rmvs/rmvs_step.f90 b/src/rmvs/rmvs_step.f90 index 673188f0e..76c5a802c 100644 --- a/src/rmvs/rmvs_step.f90 +++ b/src/rmvs/rmvs_step.f90 @@ -483,8 +483,10 @@ subroutine rmvs_make_planetocentric(param, cb, pl, tp) 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) + 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 call tpenci%set_mu(cbenci) @@ -554,7 +556,8 @@ subroutine rmvs_peri_tp(tp, pl, t, dt, lfirst, inner_index, ipleP, param) id2 = tp%id(i) xh2(:) = xpc(:, i) + xh1(:) vh2(:) = xpc(:, i) + vh1(:) - call rmvs_io_write_encounter(t, id1, id2, mu, 0.0_DP, rpl, 0.0_DP, xh1(:), xh2(:), vh1(:), vh2(:), param%enc_out) + call rmvs_io_write_encounter(t, id1, id2, mu, 0.0_DP, rpl, 0.0_DP, & + xh1(:), xh2(:), vh1(:), vh2(:), param%enc_out) end if if (tp%lperi(i)) then if (peri < tp%peri(i)) then diff --git a/src/setup/setup.f90 b/src/setup/setup.f90 index 3f1932d24..e6678effa 100644 --- a/src/setup/setup.f90 +++ b/src/setup/setup.f90 @@ -106,12 +106,15 @@ module subroutine setup_initialize_particle_info_system(self, param) associate(cb => self%cb, pl => self%pl, npl => self%pl%nbody, tp => self%tp, ntp => self%tp%nbody) - call cb%info%set_value(particle_type=CB_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", origin_time=param%t0, origin_xh=[0.0_DP, 0.0_DP, 0.0_DP], origin_vh=[0.0_DP, 0.0_DP, 0.0_DP]) + call cb%info%set_value(particle_type=CB_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", & + origin_time=param%t0, origin_xh=[0.0_DP, 0.0_DP, 0.0_DP], origin_vh=[0.0_DP, 0.0_DP, 0.0_DP]) do i = 1, self%pl%nbody - call pl%info(i)%set_value(particle_type=PL_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", origin_time=param%t0, origin_xh=self%pl%xh(:,i), origin_vh=self%pl%vh(:,i)) + call pl%info(i)%set_value(particle_type=PL_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", & + origin_time=param%t0, origin_xh=self%pl%xh(:,i), origin_vh=self%pl%vh(:,i)) end do do i = 1, self%tp%nbody - call tp%info(i)%set_value(particle_type=TP_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", origin_time=param%t0, origin_xh=self%tp%xh(:,i), origin_vh=self%tp%vh(:,i)) + call tp%info(i)%set_value(particle_type=TP_TYPE_NAME, status="ACTIVE", origin_type="Initial conditions", & + origin_time=param%t0, origin_xh=self%tp%xh(:,i), origin_vh=self%tp%vh(:,i)) end do end associate diff --git a/src/symba/symba_collision.f90 b/src/symba/symba_collision.f90 index 8b429a83d..6ebd8532c 100644 --- a/src/symba/symba_collision.f90 +++ b/src/symba/symba_collision.f90 @@ -312,7 +312,8 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec vr(:) = pl%vb(:, i) - pl%vb(:, j) rlim = pl%radius(i) + pl%radius(j) Gmtot = pl%Gmass(i) + pl%Gmass(j) - lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), Gmtot, rlim, dt, self%lvdotr(k)) + lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), & + Gmtot, rlim, dt, self%lvdotr(k)) end do else do concurrent(k = 1:nenc, lmask(k)) @@ -320,7 +321,8 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec j = self%index2(k) xr(:) = pl%xh(:, i) - tp%xh(:, j) vr(:) = pl%vb(:, i) - tp%vb(:, j) - lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%Gmass(i), pl%radius(i), dt, self%lvdotr(k)) + lcollision(k) = symba_collision_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), & + pl%Gmass(i), pl%radius(i), dt, self%lvdotr(k)) end do end if @@ -356,8 +358,8 @@ module function symba_collision_check_encounter(self, system, param, t, dt, irec write(timestr, *) t call tp%info(j)%set_value(status="DISCARDED_PLR", discard_time=t, discard_xh=tp%xh(:,j), discard_vh=tp%vh(:,j)) write(message, *) "Particle " // trim(adjustl(tp%info(j)%name)) // " (" // trim(adjustl(idstrj)) // ")" & - // " collided with massive body " // trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstri)) // ")" & - // " at t = " // trim(adjustl(timestr)) + // " collided with massive body " // trim(adjustl(pl%info(i)%name)) & + // " (" // trim(adjustl(idstri)) // ")" // " at t = " // trim(adjustl(timestr)) call io_log_one_message(FRAGGLE_LOG_OUT, message) end if end if @@ -528,7 +530,9 @@ function symba_collision_consolidate_colliders(pl, cb, param, idx_parent, collid xcrossv(:) = xc(:) .cross. vc(:) colliders%L_spin(:, j) = colliders%L_spin(:, j) + mchild * xcrossv(:) - colliders%L_spin(:, j) = colliders%L_spin(:, j) + mchild * pl%Ip(3, idx_child) * pl%radius(idx_child)**2 * pl%rot(:, idx_child) + colliders%L_spin(:, j) = colliders%L_spin(:, j) + mchild * pl%Ip(3, idx_child) & + * pl%radius(idx_child)**2 & + * pl%rot(:, idx_child) colliders%Ip(:, j) = colliders%Ip(:, j) + mchild * pl%Ip(:, idx_child) end if @@ -748,7 +752,9 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) plnew%status(1:nfrag) = NEW_PARTICLE do i = 1, nfrag write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Disruption", origin_time=param%t, name=newname, origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), collision_id=param%maxid_collision) + call plnew%info(i)%set_value(origin_type="Disruption", origin_time=param%t, name=newname, & + origin_xh=plnew%xh(:,i), & + origin_vh=plnew%vh(:,i), collision_id=param%maxid_collision) end do do i = 1, ncolliders if (colliders%idx(i) == ibiggest) then @@ -756,13 +762,16 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) else iother = ibiggest end if - call pl%info(colliders%idx(i))%set_value(status="Disruption", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) + call pl%info(colliders%idx(i))%set_value(status="Disruption", discard_time=param%t, & + discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) end do case(SUPERCATASTROPHIC) plnew%status(1:nfrag) = NEW_PARTICLE do i = 1, nfrag write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=param%t, name=newname, origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), collision_id=param%maxid_collision) + call plnew%info(i)%set_value(origin_type="Supercatastrophic", origin_time=param%t, name=newname, & + origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), & + collision_id=param%maxid_collision) end do do i = 1, ncolliders if (colliders%idx(i) == ibiggest) then @@ -770,19 +779,25 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) else iother = ibiggest end if - call pl%info(colliders%idx(i))%set_value(status="Supercatastrophic", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) + call pl%info(colliders%idx(i))%set_value(status="Supercatastrophic", discard_time=param%t, & + discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), & + discard_body_id=iother) end do case(HIT_AND_RUN_DISRUPT) call plnew%info(1)%copy(pl%info(ibiggest)) plnew%status(1) = OLD_PARTICLE do i = 2, nfrag write(newname, FRAGFMT) frag%id(i) - call plnew%info(i)%set_value(origin_type="Hit and run fragment", origin_time=param%t, name=newname, origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), collision_id=param%maxid_collision) + call plnew%info(i)%set_value(origin_type="Hit and run fragment", origin_time=param%t, name=newname, & + origin_xh=plnew%xh(:,i), origin_vh=plnew%vh(:,i), & + collision_id=param%maxid_collision) end do do i = 1, ncolliders if (colliders%idx(i) == ibiggest) cycle iother = ibiggest - call pl%info(colliders%idx(i))%set_value(status="Hit and run fragmention", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) + call pl%info(colliders%idx(i))%set_value(status="Hit and run fragmention", discard_time=param%t, & + discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), & + discard_body_id=iother) end do case(MERGED) call plnew%info(1)%copy(pl%info(ibiggest)) @@ -791,7 +806,8 @@ subroutine symba_collision_mergeaddsub(system, param, colliders, frag, status) if (colliders%idx(i) == ibiggest) cycle iother = ibiggest - call pl%info(colliders%idx(i))%set_value(status="MERGED", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=iother) + call pl%info(colliders%idx(i))%set_value(status="MERGED", discard_time=param%t, discard_xh=pl%xh(:,i), & + discard_vh=pl%vh(:,i), discard_body_id=iother) end do end select @@ -991,9 +1007,12 @@ module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, ir do write(timestr,*) t call io_log_one_message(FRAGGLE_LOG_OUT, "") - call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") - call io_log_one_message(FRAGGLE_LOG_OUT, "Collision between massive bodies detected at time t = " // trim(adjustl(timestr))) - call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & + "***********************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "Collision between massive bodies detected at time t = " // & + trim(adjustl(timestr))) + call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & + "***********************************************************") allocate(tmp_param, source=param) tmp_param%t = t if (param%lfragmentation) then @@ -1003,7 +1022,7 @@ module subroutine symba_collision_resolve_plplenc(self, system, param, t, dt, ir end if ! Destroy the collision list now that the collisions are resolved - call plplcollision_list%setup(0) + call plplcollision_list%setup(0_I8B) if ((system%pl_adds%nbody == 0) .and. (system%pl_discards%nbody == 0)) exit diff --git a/src/symba/symba_discard.f90 b/src/symba/symba_discard.f90 index 0ae89587f..385f6f25b 100644 --- a/src/symba/symba_discard.f90 +++ b/src/symba/symba_discard.f90 @@ -36,26 +36,34 @@ subroutine symba_discard_cb_pl(pl, system, param) pl%status(i) = DISCARDED_RMAX write(idstr, *) pl%id(i) write(timestr, *) param%t - write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too far from the central body at t = " // trim(adjustl(timestr)) + write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & + " too far from the central body at t = " // trim(adjustl(timestr)) call io_log_one_message(FRAGGLE_LOG_OUT, "") - call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & + "***********************************************************") call io_log_one_message(FRAGGLE_LOG_OUT, message) - call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************" // & + "***********************************************************") call io_log_one_message(FRAGGLE_LOG_OUT, "") - call pl%info(i)%set_value(status="DISCARDED_RMAX", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i)) + call pl%info(i)%set_value(status="DISCARDED_RMAX", discard_time=param%t, discard_xh=pl%xh(:,i), & + discard_vh=pl%vh(:,i)) else if ((param%rmin >= 0.0_DP) .and. (rh2 < rmin2)) then pl%ldiscard(i) = .true. pl%lcollision(i) = .false. pl%status(i) = DISCARDED_RMIN write(idstr, *) pl%id(i) write(timestr, *) param%t - write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " too close to the central body at t = " // trim(adjustl(timestr)) + write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & + " too close to the central body at t = " // trim(adjustl(timestr)) call io_log_one_message(FRAGGLE_LOG_OUT, "") - call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "************************************************************" // & + "************************************************************") call io_log_one_message(FRAGGLE_LOG_OUT, message) - call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "************************************************************" // & + "************************************************************") call io_log_one_message(FRAGGLE_LOG_OUT, "") - call pl%info(i)%set_value(status="DISCARDED_RMIN", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=cb%id) + call pl%info(i)%set_value(status="DISCARDED_RMIN", discard_time=param%t, discard_xh=pl%xh(:,i), & + discard_vh=pl%vh(:,i), discard_body_id=cb%id) else if (param%rmaxu >= 0.0_DP) then rb2 = dot_product(pl%xb(:,i), pl%xb(:,i)) vb2 = dot_product(pl%vb(:,i), pl%vb(:,i)) @@ -66,13 +74,17 @@ subroutine symba_discard_cb_pl(pl, system, param) pl%status(i) = DISCARDED_RMAXU write(idstr, *) pl%id(i) write(timestr, *) param%t - write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) + write(message, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ")" // & + " is unbound and too far from barycenter at t = " // trim(adjustl(timestr)) call io_log_one_message(FRAGGLE_LOG_OUT, "") - call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "************************************************************" // & + "************************************************************") call io_log_one_message(FRAGGLE_LOG_OUT, message) - call io_log_one_message(FRAGGLE_LOG_OUT, "***********************************************************************************************************************") + call io_log_one_message(FRAGGLE_LOG_OUT, "************************************************************" // & + "************************************************************") call io_log_one_message(FRAGGLE_LOG_OUT, "") - call pl%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i)) + call pl%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=param%t, discard_xh=pl%xh(:,i), & + discard_vh=pl%vh(:,i)) end if end if end if @@ -138,7 +150,8 @@ subroutine symba_discard_conserve_mtm(pl, system, param, ipl, lescape_body) end do Ltot(:) = Ltot(:) - cb%mass * (cb%xb(:) .cross. cb%vb(:)) system%Lescape(:) = system%Lescape(:) + Ltot(:) - if (param%lrotation) system%Lescape(:) = system%Lescape + pl%mass(ipl) * pl%radius(ipl)**2 * pl%Ip(3, ipl) * pl%rot(:, ipl) + if (param%lrotation) system%Lescape(:) = system%Lescape + pl%mass(ipl) * pl%radius(ipl)**2 & + * pl%Ip(3, ipl) * pl%rot(:, ipl) else xcom(:) = (pl%mass(ipl) * pl%xb(:, ipl) + cb%mass * cb%xb(:)) / (cb%mass + pl%mass(ipl)) @@ -305,8 +318,10 @@ subroutine symba_discard_peri_pl(pl, system, param) pl%status(i) = DISCARDED_PERI write(timestr, *) param%t write(idstr, *) pl%id(i) - write(*, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // ") perihelion distance too small at t = " // trim(adjustl(timestr)) - call pl%info(i)%set_value(status="DISCARDED_PERI", discard_time=param%t, discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=system%cb%id) + write(*, *) trim(adjustl(pl%info(i)%name)) // " (" // trim(adjustl(idstr)) // & + ") perihelion distance too small at t = " // trim(adjustl(timestr)) + call pl%info(i)%set_value(status="DISCARDED_PERI", discard_time=param%t, & + discard_xh=pl%xh(:,i), discard_vh=pl%vh(:,i), discard_body_id=system%cb%id) end if end if end if diff --git a/src/symba/symba_encounter_check.f90 b/src/symba/symba_encounter_check.f90 index 4f8495440..7281916b2 100644 --- a/src/symba/symba_encounter_check.f90 +++ b/src/symba/symba_encounter_check.f90 @@ -40,7 +40,8 @@ module function symba_encounter_check_pl(self, param, system, dt, irec) result(l if (nplt == 0) then call encounter_check_all_plpl(param, npl, pl%xh, pl%vh, pl%renc, dt, nenc, index1, index2, lvdotr) else - call encounter_check_all_plplm(param, nplm, nplt, pl%xh(:,1:nplm), pl%vh(:,1:nplm), pl%xh(:,nplm+1:npl), pl%vh(:,nplm+1:npl), pl%renc(1:nplm), pl%renc(nplm+1:npl), dt, nenc, index1, index2, lvdotr) + call encounter_check_all_plplm(param, nplm, nplt, pl%xh(:,1:nplm), pl%vh(:,1:nplm), pl%xh(:,nplm+1:npl), & + pl%vh(:,nplm+1:npl), pl%renc(1:nplm), pl%renc(nplm+1:npl), dt, nenc, index1, index2, lvdotr) end if lany_encounter = nenc > 0_I8B @@ -146,7 +147,8 @@ module function symba_encounter_check(self, param, system, dt, irec) result(lany j = self%index2(k) xr(:) = tp%xh(:,j) - pl%xh(:,i) vr(:) = tp%vb(:,j) - pl%vb(:,i) - call encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%renc(i), dt, lencounter(lidx), self%lvdotr(k)) + call encounter_check_one(xr(1), xr(2), xr(3), vr(1), vr(2), vr(3), pl%renc(i), dt, & + lencounter(lidx), self%lvdotr(k)) if (lencounter(lidx)) then rlim2 = (pl%radius(i))**2 rji2 = dot_product(xr(:), xr(:))! Check to see if these are physically overlapping bodies first, which we should ignore diff --git a/src/symba/symba_io.f90 b/src/symba/symba_io.f90 index 3e217bc6a..19e8635db 100644 --- a/src/symba/symba_io.f90 +++ b/src/symba/symba_io.f90 @@ -74,7 +74,8 @@ module subroutine symba_io_param_reader(self, unit, iotype, v_list, iostat, ioms param_value = io_get_token(line, ifirst, ilast, iostat) read(param_value, *) param%seed(i) end do - param%seed(nseeds_from_file+1:nseeds) = [(param%seed(1) - param%seed(nseeds_from_file) + i, i=nseeds_from_file+1, nseeds)] + param%seed(nseeds_from_file+1:nseeds) = [(param%seed(1) - param%seed(nseeds_from_file) + i, & + i=nseeds_from_file+1, nseeds)] end if seed_set = .true. end select @@ -192,9 +193,9 @@ module subroutine symba_io_write_discard(self, param) end if select case(out_stat) case('APPEND') - open(unit = LUN, file = param%discard_out, status = 'OLD', position = 'APPEND', form = 'FORMATTED', err = 667, iomsg = errmsg) + open(unit=LUN, file=param%discard_out, status='OLD', position='APPEND', form='FORMATTED', err=667, iomsg=errmsg) case('NEW', 'REPLACE', 'UNKNOWN') - open(unit = LUN, file = param%discard_out, status = param%out_stat, form = 'FORMATTED', err = 667, iomsg = errmsg) + open(unit=LUN, file=param%discard_out, status=param%out_stat, form='FORMATTED', err=667, iomsg=errmsg) case default write(*,*) 'Invalid status code for OUT_STAT: ',trim(adjustl(param%out_stat)) call util_exit(FAILURE) @@ -205,7 +206,7 @@ module subroutine symba_io_write_discard(self, param) call pl_adds%pv2v(param) end if - write(LUN, HDRFMT, err = 667, iomsg = errmsg) param%t, pl_discards%nbody, param%lbig_discard + write(LUN, HDRFMT, err=667, iomsg=errmsg) param%t, pl_discards%nbody, param%lbig_discard iadd = 1 isub = 1 do while (iadd <= pl_adds%nbody) @@ -213,9 +214,9 @@ module subroutine symba_io_write_discard(self, param) nsub = pl_discards%ncomp(isub) do j = 1, nadd if (iadd <= pl_adds%nbody) then - write(LUN, NAMEFMT, err = 667, iomsg = errmsg) ADD, pl_adds%id(iadd), pl_adds%status(iadd) - write(LUN, VECFMT, err = 667, iomsg = errmsg) pl_adds%xh(1, iadd), pl_adds%xh(2, iadd), pl_adds%xh(3, iadd) - write(LUN, VECFMT, err = 667, iomsg = errmsg) pl_adds%vh(1, iadd), pl_adds%vh(2, iadd), pl_adds%vh(3, iadd) + write(LUN, NAMEFMT, err=667, iomsg=errmsg) ADD, pl_adds%id(iadd), pl_adds%status(iadd) + write(LUN, VECFMT, err=667, iomsg=errmsg) pl_adds%xh(1, iadd), pl_adds%xh(2, iadd), pl_adds%xh(3, iadd) + write(LUN, VECFMT, err=667, iomsg=errmsg) pl_adds%vh(1, iadd), pl_adds%vh(2, iadd), pl_adds%vh(3, iadd) else exit end if @@ -223,9 +224,9 @@ module subroutine symba_io_write_discard(self, param) end do do j = 1, nsub if (isub <= pl_discards%nbody) then - write(LUN, NAMEFMT, err = 667, iomsg = errmsg) SUB, pl_discards%id(isub), pl_discards%status(isub) - write(LUN, VECFMT, err = 667, iomsg = errmsg) pl_discards%xh(1, isub), pl_discards%xh(2, isub), pl_discards%xh(3, isub) - write(LUN, VECFMT, err = 667, iomsg = errmsg) pl_discards%vh(1, isub), pl_discards%vh(2, isub), pl_discards%vh(3, isub) + write(LUN,NAMEFMT,err=667,iomsg=errmsg) SUB, pl_discards%id(isub), pl_discards%status(isub) + write(LUN,VECFMT,err=667,iomsg=errmsg) pl_discards%xh(1,isub), pl_discards%xh(2,isub), pl_discards%xh(3,isub) + write(LUN,VECFMT,err=667,iomsg=errmsg) pl_discards%vh(1,isub), pl_discards%vh(2,isub), pl_discards%vh(3,isub) else exit end if diff --git a/src/symba/symba_setup.f90 b/src/symba/symba_setup.f90 index cecb4d806..cda5214d0 100644 --- a/src/symba/symba_setup.f90 +++ b/src/symba/symba_setup.f90 @@ -18,9 +18,9 @@ module subroutine symba_setup_initialize_system(self, param) ! Call parent method associate(system => self) call helio_setup_initialize_system(system, param) - call system%pltpenc_list%setup(0) - call system%plplenc_list%setup(0) - call system%plplcollision_list%setup(0) + call system%pltpenc_list%setup(0_I8B) + call system%plplenc_list%setup(0_I8B) + call system%plplcollision_list%setup(0_I8B) end associate return diff --git a/src/symba/symba_step.f90 b/src/symba/symba_step.f90 index 7bd2ff562..f2daa8846 100644 --- a/src/symba/symba_step.f90 +++ b/src/symba/symba_step.f90 @@ -110,7 +110,8 @@ module subroutine symba_step_set_recur_levels_system(self, ireci) ! Internals integer(I4B) :: k, irecp - associate(system => self, plplenc_list => self%plplenc_list, pltpenc_list => self%pltpenc_list, npl => self%pl%nbody, ntp => self%tp%nbody) + associate(system => self, plplenc_list => self%plplenc_list, pltpenc_list => self%pltpenc_list, & + npl => self%pl%nbody, ntp => self%tp%nbody) select type(pl => self%pl) class is (symba_pl) select type(tp => self%tp) @@ -119,8 +120,16 @@ module subroutine symba_step_set_recur_levels_system(self, ireci) if (npl >0) where(pl%levelg(1:npl) == irecp) pl%levelg(1:npl) = ireci if (ntp > 0) where(tp%levelg(1:ntp) == irecp) tp%levelg(1:ntp) = ireci - if (plplenc_list%nenc > 0) where(plplenc_list%level(1:plplenc_list%nenc) == irecp) plplenc_list%level(1:plplenc_list%nenc) = ireci - if (pltpenc_list%nenc > 0) where(pltpenc_list%level(1:pltpenc_list%nenc) == irecp) pltpenc_list%level(1:pltpenc_list%nenc) = ireci + if (plplenc_list%nenc > 0) then + where(plplenc_list%level(1:plplenc_list%nenc) == irecp) + plplenc_list%level(1:plplenc_list%nenc) = ireci + endwhere + end if + if (pltpenc_list%nenc > 0) then + where(pltpenc_list%level(1:pltpenc_list%nenc) == irecp) + pltpenc_list%level(1:pltpenc_list%nenc) = ireci + endwhere + end if system%irec = ireci @@ -173,7 +182,8 @@ module recursive subroutine symba_step_recur_system(self, param, t, ireci) nloops = NTENC end if do j = 1, nloops - lencounter = plplenc_list%encounter_check(param, system, dtl, irecp) .or. pltpenc_list%encounter_check(param, system, dtl, irecp) + lencounter = plplenc_list%encounter_check(param, system, dtl, irecp) & + .or. pltpenc_list%encounter_check(param, system, dtl, irecp) call plplenc_list%kick(system, dth, irecp, 1) call pltpenc_list%kick(system, dth, irecp, 1) @@ -247,10 +257,10 @@ module subroutine symba_step_reset_system(self, param) pl%ldiscard(1:npl) = .false. pl%lmask(1:npl) = .true. nenc_old = system%plplenc_list%nenc - call system%plplenc_list%setup(0) + call system%plplenc_list%setup(0_I8B) call system%plplenc_list%setup(nenc_old) - system%plplenc_list%nenc = 0 - call system%plplcollision_list%setup(0) + system%plplenc_list%nenc = 0_I8B + call system%plplcollision_list%setup(0_I8B) call system%pl%set_renc(0) end if @@ -261,9 +271,9 @@ module subroutine symba_step_reset_system(self, param) tp%lmask(1:ntp) = .true. tp%ldiscard(1:npl) = .false. nenc_old = system%pltpenc_list%nenc - call system%pltpenc_list%setup(0) + call system%pltpenc_list%setup(0_I8B) call system%pltpenc_list%setup(nenc_old) - system%pltpenc_list%nenc = 0 + system%pltpenc_list%nenc = 0_I8B end if call system%pl_adds%setup(0, param) diff --git a/src/symba/symba_util.f90 b/src/symba/symba_util.f90 index 892008aea..0cbdc9a4a 100644 --- a/src/symba/symba_util.f90 +++ b/src/symba/symba_util.f90 @@ -569,7 +569,8 @@ module subroutine symba_util_peri_pl(self, system, param) if (pl%isperi(i) == -1) then if (vdotr >= 0.0_DP) then pl%isperi(i) = 0 - CALL orbel_xv2aeq(system%Gmtot, pl%xb(1,i), pl%xb(2,i), pl%xb(3,i), pl%vb(1,i), pl%vb(2,i), pl%vb(3,i), pl%atp(i), e, pl%peri(i)) + CALL orbel_xv2aeq(system%Gmtot, pl%xb(1,i), pl%xb(2,i), pl%xb(3,i), pl%vb(1,i), pl%vb(2,i), pl%vb(3,i),& + pl%atp(i), e, pl%peri(i)) end if else if (vdotr > 0.0_DP) then @@ -603,7 +604,8 @@ module subroutine symba_util_rearray_pl(self, system, param) logical, dimension(:), allocatable :: lmask, ldump_mask class(symba_plplenc), allocatable :: plplenc_old logical :: lencounter - integer(I4B), dimension(:), allocatable :: levelg_orig_pl, levelm_orig_pl, levelg_orig_tp, levelm_orig_tp, nplenc_orig_pl, nplenc_orig_tp, ntpenc_orig_pl + integer(I4B), dimension(:), allocatable :: levelg_orig_pl, levelm_orig_pl, levelg_orig_tp, levelm_orig_tp + integer(I4B), dimension(:), allocatable :: nplenc_orig_pl, nplenc_orig_tp, ntpenc_orig_pl associate(pl => self, pl_adds => system%pl_adds) diff --git a/src/util/util_peri.f90 b/src/util/util_peri.f90 index a81748d74..d451e3dbd 100644 --- a/src/util/util_peri.f90 +++ b/src/util/util_peri.f90 @@ -29,7 +29,8 @@ module subroutine util_peri_tp(self, system, param) if (tp%isperi(i) == -1) then if (vdotr(i) >= 0.0_DP) then tp%isperi(i) = 0 - call orbel_xv2aeq(tp%mu(i), tp%xh(1,i), tp%xh(2,i), tp%xh(3,i), tp%vh(1,i), tp%vh(2,i), tp%vh(3,i), tp%atp(i), e, tp%peri(i)) + call orbel_xv2aeq(tp%mu(i), tp%xh(1,i), tp%xh(2,i), tp%xh(3,i), tp%vh(1,i), tp%vh(2,i), tp%vh(3,i), & + tp%atp(i), e, tp%peri(i)) end if else if (vdotr(i) > 0.0_DP) then @@ -45,7 +46,8 @@ module subroutine util_peri_tp(self, system, param) if (tp%isperi(i) == -1) then if (vdotr(i) >= 0.0_DP) then tp%isperi(i) = 0 - call orbel_xv2aeq(system%Gmtot, tp%xb(1,i), tp%xb(2,i), tp%xb(3,i), tp%vb(1,i), tp%vb(2,i), tp%vb(3,i), tp%atp(i), e, tp%peri(i)) + call orbel_xv2aeq(system%Gmtot, tp%xb(1,i), tp%xb(2,i), tp%xb(3,i), tp%vb(1,i), tp%vb(2,i), tp%vb(3,i), & + tp%atp(i), e, tp%peri(i)) end if else if (vdotr(i) > 0.0_DP) then diff --git a/src/util/util_set.f90 b/src/util/util_set.f90 index 98ee4f497..99dff1e0d 100644 --- a/src/util/util_set.f90 +++ b/src/util/util_set.f90 @@ -98,7 +98,8 @@ module subroutine util_set_mu_tp(self, cb) return end subroutine util_set_mu_tp - module subroutine util_set_particle_info(self, name, particle_type, status, origin_type, origin_time, collision_id, origin_xh, origin_vh, discard_time, discard_xh, discard_vh, discard_body_id) + module subroutine util_set_particle_info(self, name, particle_type, status, origin_type, origin_time, collision_id, origin_xh,& + origin_vh, discard_time, discard_xh, discard_vh, discard_body_id) !! author: David A. Minton !! !! Sets one or more values of the particle information metadata object diff --git a/src/walltime/walltime.f90 b/src/walltime/walltime.f90 index a8d82c2b5..87b0a33c2 100644 --- a/src/walltime/walltime.f90 +++ b/src/walltime/walltime.f90 @@ -39,7 +39,8 @@ module subroutine walltime_report(self, nsubsteps, message) integer(I4B), intent(in) :: nsubsteps !! Number of substeps used to compute the time per step character(len=*), intent(in) :: message !! Message to prepend to the wall time terminal output ! Internals - character(len=*), parameter :: walltimefmt = '" Total wall time: ", es12.5, "; Interval wall time: ", es12.5, "; Interval wall time/step: ", es12.5' + character(len=*), parameter :: walltimefmt = '" Total wall time: ", es12.5, "; Interval wall time: ", es12.5, ";' //& + 'Interval wall time/step: ", es12.5' character(len=STRMAX) :: fmt integer(I8B) :: count_delta_step, count_delta_main, count_now real(DP) :: wall_main !! Value of total elapsed time at the end of a timed step @@ -196,7 +197,8 @@ module subroutine walltime_interaction_adapt(self, param, ninteractions, pl) write(cstr,*) self%count_stop_step - self%count_start_step - call io_log_one_message(logfile, adjustl(lstyle) // " " // trim(adjustl(cstr)) // " " // trim(adjustl(nstr)) // " " // trim(adjustl(mstr))) + call io_log_one_message(logfile, adjustl(lstyle) // " " // trim(adjustl(cstr)) // " " // & + trim(adjustl(nstr)) // " " // trim(adjustl(mstr))) if (self%stage == 2) then if (ladvanced_final) then @@ -204,7 +206,8 @@ module subroutine walltime_interaction_adapt(self, param, ninteractions, pl) else lstyle = standardstyle end if - call io_log_one_message(logfile, trim(adjustl(self%loopname)) // ": the fastest loop method tested is " // trim(adjustl(lstyle))) + call io_log_one_message(logfile, trim(adjustl(self%loopname)) // & + ": the fastest loop method tested is " // trim(adjustl(lstyle))) end if return