Moved time variable from param to nbody_system so that the system his…

…tory dumps have the correct time associated with them

src/discard/discard.f90

@@ -135,22 +135,22 @@ subroutine discard_cb_tp(tp, system, param)
                if ((param%rmax >= 0.0_DP) .and. (rh2 > rmax2)) then
                   tp%status(i) = DISCARDED_RMAX
                   write(idstr, *) tp%id(i)
-                  write(timestr, *) param%t
+                  write(timestr, *) system%t
                   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), &
+                  call tp%info(i)%set_value(status="DISCARDED_RMAX", discard_time=system%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(timestr, *) system%t
                   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), &
+                  call tp%info(i)%set_value(status="DISCARDED_RMIN", discard_time=system%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))
@@ -159,12 +159,12 @@ subroutine discard_cb_tp(tp, system, param)
                   if ((energy > 0.0_DP) .and. (rb2 > rmaxu2)) then
                      tp%status(i) = DISCARDED_RMAXU
                      write(idstr, *) tp%id(i)
-                     write(timestr, *) param%t
+                     write(timestr, *) system%t
                      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), &
+                     call tp%info(i)%set_value(status="DISCARDED_RMAXU", discard_time=system%t, discard_xh=tp%xh(:,i), &
                                                discard_vh=tp%vh(:,i))
                   end if
                end if
@@ -194,7 +194,7 @@ subroutine discard_peri_tp(tp, system, param)
       real(DP), dimension(NDIM) :: dx
       character(len=STRMAX) :: idstr, timestr
 
-      associate(cb => system%cb, ntp => tp%nbody, pl => system%pl, npl => system%pl%nbody, t => param%t)
+      associate(cb => system%cb, ntp => tp%nbody, pl => system%pl, npl => system%pl%nbody, t => system%t)
          call tp%get_peri(system, param)
          do i = 1, ntp
             if (tp%status(i) == ACTIVE) then
@@ -211,11 +211,11 @@ subroutine discard_peri_tp(tp, system, param)
                         (tp%peri(i) <= param%qmin)) then
                         tp%status(i) = DISCARDED_PERI
                         write(idstr, *) tp%id(i)
-                        write(timestr, *) param%t
+                        write(timestr, *) system%t
                         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), &
+                        call tp%info(i)%set_value(status="DISCARDED_PERI", discard_time=system%t, discard_xh=tp%xh(:,i), &
                                                   discard_vh=tp%vh(:,i), discard_body_id=pl%id(j))
                      end if
                   end if
@@ -246,7 +246,7 @@ subroutine discard_pl_tp(tp, system, param)
       real(DP), dimension(NDIM) :: dx, dv
       character(len=STRMAX) :: idstri, idstrj, timestr
 
-      associate(ntp => tp%nbody, pl => system%pl, npl => system%pl%nbody, t => param%t, dt => param%dt)
+      associate(ntp => tp%nbody, pl => system%pl, npl => system%pl%nbody, t => system%t, dt => param%dt)
          do i = 1, ntp
             if (tp%status(i) == ACTIVE) then
                do j = 1, npl
@@ -260,12 +260,12 @@ subroutine discard_pl_tp(tp, system, param)
                      pl%ldiscard(j) = .true.
                      write(idstri, *) tp%id(i)
                      write(idstrj, *) pl%id(j)
-                     write(timestr, *) param%t
+                     write(timestr, *) system%t
                      write(*, *) "Test particle " // trim(adjustl(tp%info(i)%name)) // " ("  // trim(adjustl(idstri)) // ")" &
                                                   // "  too close to massive body " // trim(adjustl(pl%info(j)%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), &
+                     call tp%info(i)%set_value(status="DISCARDED_PLR", discard_time=system%t, discard_xh=tp%xh(:,i), &
                                                discard_vh=tp%vh(:,i), discard_body_id=pl%id(j))
                      exit
                   end if

src/io/io.f90

@@ -35,7 +35,7 @@ module subroutine io_compact_output(self, param, timer)
 
       select type(timer)
       class is (walltimer)
-         formatted_output = fmt("ILOOP",param%iloop) // fmt("T",param%t) // fmt("NPL",self%pl%nbody) // fmt("NTP",self%tp%nbody) 
+         formatted_output = fmt("ILOOP",param%iloop) // fmt("T",self%t) // fmt("NPL",self%pl%nbody) // fmt("NTP",self%tp%nbody) 
          select type(pl => self%pl)
          class is (symba_pl)
             formatted_output = formatted_output // fmt("NPLM",pl%nplm)
@@ -250,7 +250,7 @@ module subroutine io_dump_system(self, param)
       dump_param%in_netcdf = trim(adjustl(DUMP_NC_FILE(idx)))
       dump_param%nciu%id_chunk = self%pl%nbody + self%tp%nbody
       dump_param%nciu%time_chunk = 1
-      dump_param%tstart = param%t
+      dump_param%tstart = self%t
 
       call dump_param%dump(param_file_name)
 

src/main/swiftest_driver.f90

@@ -66,7 +66,7 @@ program swiftest_driver
    call setup_construct_system(nbody_system, param)
    call param%read_in(param_file_name)
 
-   associate(t               => param%t, &
+   associate(t               => nbody_system%t, &
              t0              => param%t0, &
              tstart          => param%tstart, &
              dt              => param%dt, &
@@ -137,13 +137,13 @@ program swiftest_driver
                   idump = dump_cadence
                end if
 
-               tfrac = (param%t - param%t0) / (param%tstop - param%t0)
+               tfrac = (t - t0) / (tstop - t0)
 
                select type(pl => nbody_system%pl)
                class is (symba_pl)
-                  write(display_unit, symbastatfmt) param%t, tfrac, pl%nplm, pl%nbody, nbody_system%tp%nbody
+                  write(display_unit, symbastatfmt) t, tfrac, pl%nplm, pl%nbody, nbody_system%tp%nbody
                class default
-                  write(display_unit, statusfmt) param%t, tfrac, pl%nbody, nbody_system%tp%nbody
+                  write(display_unit, statusfmt) t, tfrac, pl%nbody, nbody_system%tp%nbody
                end select
                if (param%lenergy) call nbody_system%conservation_report(param, lterminal=.true.)
                call integration_timer%report(message="Integration steps:", unit=display_unit, nsubsteps=istep_out)

src/modules/swiftest_classes.f90

@@ -37,7 +37,6 @@ module swiftest_classes
       integer(I4B)       :: maxid                = -1                 !! The current maximum particle id number 
       integer(I4B)       :: maxid_collision      = 0                  !! The current maximum collision id number
       real(DP)           :: t0                   =  0.0_DP            !! Integration reference time
-      real(DP)           :: t                    = -1.0_DP            !! Integration current time
       real(DP)           :: tstart               = -1.0_DP            !! Integration start time
       real(DP)           :: tstop                = -1.0_DP            !! Integration stop time
       real(DP)           :: dt                   = -1.0_DP            !! Time step
@@ -347,6 +346,7 @@ module swiftest_classes
       class(swiftest_tp), allocatable :: tp                   !! Test particle data structure
       class(swiftest_tp), allocatable :: tp_discards          !! Discarded test particle data structure
       class(swiftest_pl), allocatable :: pl_discards          !! Discarded massive body particle data structure
+      real(DP)                        :: t = -1.0_DP          !! Integration current time
       real(DP)                        :: GMtot = 0.0_DP       !! Total system mass - used for barycentric coordinate conversion
       real(DP)                        :: ke_orbit = 0.0_DP    !! System orbital kinetic energy
       real(DP)                        :: ke_spin = 0.0_DP     !! System spin kinetic energy

src/netcdf/netcdf.f90

@@ -793,7 +793,7 @@ module subroutine netcdf_read_hdr_system(self, iu, param)
 
       tslot = int(param%ioutput, kind=I4B) + 1
       call check( nf90_inquire_dimension(iu%ncid, iu%id_dimid, len=idmax), "netcdf_read_frame_system nf90_inquire_dimension id_dimid"  )
-      call check( nf90_get_var(iu%ncid, iu%time_varid, param%t,       start=[tslot]), "netcdf_read_hdr_system nf90_getvar time_varid"  )
+      call check( nf90_get_var(iu%ncid, iu%time_varid, self%t, start=[tslot]), "netcdf_read_hdr_system nf90_getvar time_varid"  )
 
       allocate(gmtemp(idmax))
       allocate(tpmask(idmax))
@@ -1428,7 +1428,7 @@ module subroutine netcdf_write_hdr_system(self, iu, param)
 
       tslot = int(param%ioutput, kind=I4B) + 1
 
-      call check( nf90_put_var(iu%ncid, iu%time_varid, param%t, start=[tslot]), "netcdf_write_hdr_system nf90_put_var time_varid"  )
+      call check( nf90_put_var(iu%ncid, iu%time_varid, self%t, start=[tslot]), "netcdf_write_hdr_system nf90_put_var time_varid"  )
       call check( nf90_put_var(iu%ncid, iu%npl_varid, self%pl%nbody, start=[tslot]), "netcdf_write_hdr_system nf90_put_var npl_varid"  )
       call check( nf90_put_var(iu%ncid, iu%ntp_varid, self%tp%nbody, start=[tslot]), "netcdf_write_hdr_system nf90_put_var ntp_varid"  )
       select type(pl => self%pl)

src/rmvs/rmvs_discard.f90

@@ -29,7 +29,7 @@ module subroutine rmvs_discard_tp(self, system, param)
 
       if (self%nbody == 0) return
 
-      associate(tp => self, ntp => self%nbody, pl => system%pl, t => param%t)
+      associate(tp => self, ntp => self%nbody, pl => system%pl, t => system%t)
          do i = 1, ntp
             associate(iplperP => tp%plperP(i))
                if ((tp%status(i) == ACTIVE) .and. (tp%lperi(i))) then