diff --git a/autotest/__snapshots__/test_prt_dry/test_mf6model[0-prtdry].npy b/autotest/__snapshots__/test_prt_dry/test_mf6model[0-prtdry].npy
index 5e8e184fe31..7a4f24e0652 100644
Binary files a/autotest/__snapshots__/test_prt_dry/test_mf6model[0-prtdry].npy and b/autotest/__snapshots__/test_prt_dry/test_mf6model[0-prtdry].npy differ
diff --git a/autotest/__snapshots__/test_prt_dry/test_mf6model[2-prtdry_drop].npy b/autotest/__snapshots__/test_prt_dry/test_mf6model[2-prtdry_drop].npy
index de1b2802bb4..b2dbc3f7e0f 100644
Binary files a/autotest/__snapshots__/test_prt_dry/test_mf6model[2-prtdry_drop].npy and b/autotest/__snapshots__/test_prt_dry/test_mf6model[2-prtdry_drop].npy differ
diff --git a/autotest/__snapshots__/test_prt_dry/test_mf6model[4-prtdry_stay].npy b/autotest/__snapshots__/test_prt_dry/test_mf6model[4-prtdry_stay].npy
index 5a921e5f6f3..bb0c0fdcfcc 100644
Binary files a/autotest/__snapshots__/test_prt_dry/test_mf6model[4-prtdry_stay].npy and b/autotest/__snapshots__/test_prt_dry/test_mf6model[4-prtdry_stay].npy differ
diff --git a/autotest/test_prt_dry.py b/autotest/test_prt_dry.py
index d68af937af1..edadf697fcc 100644
--- a/autotest/test_prt_dry.py
+++ b/autotest/test_prt_dry.py
@@ -70,6 +70,17 @@
 for i in range(nper):
     period_data.append((perlen[i], nstp[i], tsmult[i]))
 
+user_time = 100.0
+
+# particles are released in layer 0
+offsets = [
+    (-1, 0, 0),
+    (-1, -1, 0),
+    (-1, 1, 0),
+    (-1, 0, -1),
+    (-1, 0, 1),
+]
+
 
 def build_gwf_sim(name, gwf_ws, mf6, newton=False):
     sim = flopy.mf6.MFSimulation(
@@ -267,15 +278,6 @@ def build_prt_sim(name, gwf, prt_ws, mf6, drape=False, dry_tracking_method=False
 
     flopy.mf6.ModflowPrtmip(prt, porosity=porosity, pname="mip")
 
-    # particles are released in layer 0
-    offsets = [
-        (-1, 0, 0),
-        (-1, -1, 0),
-        (-1, 1, 0),
-        (-1, 0, -1),
-        (-1, 0, 1),
-    ]
-
     lay = 1
     row, col = (int(nrow / 4), int(ncol / 4))
     prp_cells = [(lay + k, row + i, col + j) for k, i, j in offsets]
@@ -319,6 +321,8 @@ def build_prt_sim(name, gwf, prt_ws, mf6, drape=False, dry_tracking_method=False
         trackcsv_filerecord=[trackcsvfile],
         saverecord=[("BUDGET", "ALL")],
         pname="oc",
+        ntracktimes=1 if "stay" in name else 0,
+        tracktimes=[(user_time,)] if "stay" in name else None,
     )
 
     rel_prt_folder = os.path.relpath(gwf_ws, start=prt_ws)
@@ -368,11 +372,31 @@ def check_output(idx, test, snapshot):
     strtpts = pls[pls.ireason == 0]
 
     # compare to expected results
-    decimals = 1 if "drop" in name else 2
-    actual = pls.drop(["name", "icell"], axis=1).round(decimals).reset_index(drop=True)
-    # ignore particle 4, it terminates early with optimization=2 when built with ifort
-    if "drop" in name:
+    places = 1 if "drop" in name else 2
+    actual = pls.drop(["name", "icell"], axis=1).round(places).reset_index(drop=True)
+    nparts = len(offsets)  # number of particles
+
+    if "drape" in name:
+        assert len(actual[actual.ireason == 0]) == nparts  # release
+    elif "drop" in name:
+        # ignore particle 4, it terminates early when
+        # mf6 is built with optimization=2 with ifort
         actual = actual.drop(actual[actual.irpt == 4].index)
+        nparts -= 1
+        assert len(actual[actual.ireason == 0]) == nparts  # release
+    elif "stop" in name:
+        assert len(actual[actual.ireason == 0]) == nparts  # release
+    elif "stay" in name:
+        assert len(actual[actual.ireason == 0]) == nparts  # release
+        assert len(actual[actual.t == user_time]) == nparts  # user time
+    else:
+        # immediate termination, permanently unreleased
+        assert len(actual) == nparts
+
+    # in all cases, all particles should have a termination event
+    assert len(actual[actual.ireason == 3]) == nparts
+
+    # snapshot comparison
     assert snapshot == actual.to_records(index=False)
 
     plot_pathlines = False
diff --git a/doc/ReleaseNotes/develop.tex b/doc/ReleaseNotes/develop.tex
index ca398879d04..c268269bc48 100644
--- a/doc/ReleaseNotes/develop.tex
+++ b/doc/ReleaseNotes/develop.tex
@@ -24,6 +24,7 @@
 	\item GWT, GWE and PRT FMI Packages can now read a GWF Model's binary grid file via a new GWFGRID entry. This allows FMI Packages to perform the same grid equivalence checks as exchanges, which guarantees identical error-checking behavior whether a GWT, GWE or PRT Model is coupled to a GWF Model or running as a post-processor. The GWFGRID file entry is optional but recommended. A future version may make the GWFGRID entry mandatory.
 	\item A regression was recently introduced into the PRT model's generalized tracking method, in which a coordinate transformation was carried out prematurely. This could result in incorrect particle positions in and near quad-refined cells. This bug has been fixed.
 	\item The PRT model previously allowed particles to be released at any time. Release times falling outside the bounds of the simulation's time discretization could produce undefined behavior. Any release times occurring before the simulation begins (i.e. negative times) will now be skipped with a warning message. If EXTEND\_TRACKING is not enabled, release times occurring after the end of the simulation will now be skipped with a warning message as well. If EXTEND\_TRACKING is enabled, release times after the end of the simulation are allowed.
+	\item The PRT Model did not previously report all expected tracking events. In particular, time step end and termination events could go unreported with the DRY\_TRACKING\_METHOD options DROP and STAY (only relevant for Newton formulation models), and termination events were not always reported at the end of the simulation. Reporting has been corrected for the cases identified.
 	\item A profiling module is added for more enhanced performance diagnostics of the program. It can be activated through the PROFILE\_OPTION in the simulation name file.
 \end{itemize}
 
diff --git a/src/Model/ModelUtilities/TrackFile.f90 b/src/Model/ModelUtilities/TrackFile.f90
index a09be74b7a8..e6bbcda56a9 100644
--- a/src/Model/ModelUtilities/TrackFile.f90
+++ b/src/Model/ModelUtilities/TrackFile.f90
@@ -47,12 +47,13 @@ module TrackFileModule
   !!     1: active
   !!     2: terminated at boundary face
   !!     3: terminated in weak sink cell
-  !!     4: terminated in weak source cell**
+  !!     4: terminated in weak source cell
   !!     5: terminated in cell with no exit face
   !!     6: terminated in cell with specified zone number
   !!     7: terminated in inactive cell
   !!     8: permanently unreleased***
   !!     9: terminated in subcell with no exit face*****
+  !!     10: terminated upon simulation's end
   !!
   !! PRT shares the same status enumeration as MODPATH 7. However, some
   !! don't apply to PRT; for instance, MODPATH 7 distinguishes forwards
diff --git a/src/Model/ParticleTracking/prt-prp.f90 b/src/Model/ParticleTracking/prt-prp.f90
index a3e657bb02a..ae3c154a83e 100644
--- a/src/Model/ParticleTracking/prt-prp.f90
+++ b/src/Model/ParticleTracking/prt-prp.f90
@@ -436,17 +436,11 @@ subroutine release(this, ip, trelease)
     type(ParticleType), pointer :: particle
 
     call this%initialize_particle(particle, ip, trelease)
-
-    ! Increment cumulative particle count
     np = this%nparticles + 1
     this%nparticles = np
-
-    ! Save the particle to the store
-    call this%particles%save_particle(particle, np)
+    call this%particles%put(particle, np)
     deallocate (particle)
-
-    ! Accumulate mass for release point
-    this%rptm(ip) = this%rptm(ip) + DONE
+    this%rptm(ip) = this%rptm(ip) + DONE ! TODO configurable mass
 
   end subroutine release
 
@@ -486,7 +480,6 @@ subroutine initialize_particle(this, particle, ip, trelease)
     end select
     particle%ilay = ilay
     particle%izone = this%rptzone(ic)
-
     particle%istatus = 0
     ! If the cell is inactive, either drape the particle
     ! to the top-most active cell beneath it if drape is
@@ -495,8 +488,9 @@ subroutine initialize_particle(this, particle, ip, trelease)
       ic_old = ic
       if (this%idrape > 0) then
         call this%dis%highest_active(ic, this%ibound)
-        if (ic == ic_old .or. this%ibound(ic) == 0) &
+        if (ic == ic_old .or. this%ibound(ic) == 0) then
           particle%istatus = 8
+        end if
       else
         particle%istatus = 8
       end if
diff --git a/src/Model/ParticleTracking/prt.f90 b/src/Model/ParticleTracking/prt.f90
index a2c1fc3e927..43a19937af5 100644
--- a/src/Model/ParticleTracking/prt.f90
+++ b/src/Model/ParticleTracking/prt.f90
@@ -58,7 +58,7 @@ module PrtModule
     real(DP), dimension(:), pointer, contiguous :: massstoold => null() !< particle mass storage in cells, old value
     real(DP), dimension(:), pointer, contiguous :: ratesto => null() !< particle mass storage rate in cells
   contains
-    ! -- Override BaseModelType procs
+    ! Override BaseModelType procs
     procedure :: model_df => prt_df
     procedure :: model_ar => prt_ar
     procedure :: model_rp => prt_rp
@@ -69,7 +69,7 @@ module PrtModule
     procedure :: model_da => prt_da
     procedure :: model_solve => prt_solve
 
-    ! -- Private utilities
+    ! Private utilities
     procedure :: allocate_scalars
     procedure :: allocate_arrays
     procedure, private :: package_create
@@ -110,14 +110,14 @@ module PrtModule
   data PRT_MULTIPKG/'PRP6 ', '     ', '     ', '     ', '     ', & !  5
                    &45*'     '/ ! 50
 
-  ! -- size of supported model package arrays
+  ! size of supported model package arrays
   integer(I4B), parameter :: NIUNIT_PRT = PRT_NBASEPKG + PRT_NMULTIPKG
 
 contains
 
   !> @brief Create a new particle tracking model object
   subroutine prt_cr(filename, id, modelname)
-    ! -- modules
+    ! modules
     use ListsModule, only: basemodellist
     use BaseModelModule, only: AddBaseModelToList
     use ConstantsModule, only: LINELENGTH, LENPACKAGENAME
@@ -126,41 +126,41 @@ subroutine prt_cr(filename, id, modelname)
     use MemoryManagerExtModule, only: mem_set_value
     use SimVariablesModule, only: idm_context
     use GwfNamInputModule, only: GwfNamParamFoundType
-    ! -- dummy
+    ! dummy
     character(len=*), intent(in) :: filename
     integer(I4B), intent(in) :: id
     character(len=*), intent(in) :: modelname
-    ! -- local
+    ! local
     type(PrtModelType), pointer :: this
     class(BaseModelType), pointer :: model
     character(len=LENMEMPATH) :: input_mempath
     character(len=LINELENGTH) :: lst_fname
     type(GwfNamParamFoundType) :: found
 
-    ! -- Allocate a new PRT Model (this)
+    ! Allocate a new PRT Model (this)
     allocate (this)
 
-    ! -- Set this before any allocs in the memory manager can be done
+    ! Set this before any allocs in the memory manager can be done
     this%memoryPath = create_mem_path(modelname)
 
-    ! -- Allocate track control object
+    ! Allocate track control object
     allocate (this%trackctl)
 
-    ! -- Allocate scalars and add model to basemodellist
+    ! Allocate scalars and add model to basemodellist
     call this%allocate_scalars(modelname)
     model => this
     call AddBaseModelToList(basemodellist, model)
 
-    ! -- Assign variables
+    ! Assign variables
     this%filename = filename
     this%name = modelname
     this%macronym = 'PRT'
     this%id = id
 
-    ! -- Set input model namfile memory path
+    ! Set input model namfile memory path
     input_mempath = create_mem_path(modelname, 'NAM', idm_context)
 
-    ! -- Copy options from input context
+    ! Copy options from input context
     call mem_set_value(this%iprpak, 'PRINT_INPUT', input_mempath, &
                        found%print_input)
     call mem_set_value(this%iprflow, 'PRINT_FLOWS', input_mempath, &
@@ -168,21 +168,21 @@ subroutine prt_cr(filename, id, modelname)
     call mem_set_value(this%ipakcb, 'SAVE_FLOWS', input_mempath, &
                        found%save_flows)
 
-    ! -- Create the list file
+    ! Create the list file
     call this%create_lstfile(lst_fname, filename, found%list, &
                              'PARTICLE TRACKING MODEL (PRT)')
 
-    ! -- Activate save_flows if found
+    ! Activate save_flows if found
     if (found%save_flows) then
       this%ipakcb = -1
     end if
 
-    ! -- Log options
+    ! Log options
     if (this%iout > 0) then
       call this%log_namfile_options(found)
     end if
 
-    ! -- Create model packages
+    ! Create model packages
     call this%create_packages()
   end subroutine prt_cr
 
@@ -192,21 +192,21 @@ end subroutine prt_cr
   !! (2) set variables and pointers
   !<
   subroutine prt_df(this)
-    ! -- modules
+    ! modules
     use PrtPrpModule, only: PrtPrpType
-    ! -- dummy
+    ! dummy
     class(PrtModelType) :: this
-    ! -- local
+    ! local
     integer(I4B) :: ip
     class(BndType), pointer :: packobj
 
-    ! -- Define packages and utility objects
+    ! Define packages and utility objects
     call this%dis%dis_df()
     call this%fmi%fmi_df(this%dis, 1)
     call this%oc%oc_df()
     call this%budget%budget_df(NIUNIT_PRT, 'MASS', 'M')
 
-    ! -- Define packages and assign iout for time series managers
+    ! Define packages and assign iout for time series managers
     do ip = 1, this%bndlist%Count()
       packobj => GetBndFromList(this%bndlist, ip)
       call packobj%bnd_df(this%dis%nodes, this%dis)
@@ -214,7 +214,7 @@ subroutine prt_df(this)
       packobj%TasManager%iout = this%iout
     end do
 
-    ! -- Allocate model arrays
+    ! Allocate model arrays
     call this%allocate_arrays()
 
   end subroutine prt_df
@@ -225,26 +225,26 @@ end subroutine prt_df
   !! (2) allocates memory for arrays part of this model object
   !<
   subroutine prt_ar(this)
-    ! -- modules
+    ! modules
     use ConstantsModule, only: DHNOFLO
     use PrtPrpModule, only: PrtPrpType
     use PrtMipModule, only: PrtMipType
     use MethodPoolModule, only: method_dis, method_disv
-    ! -- dummy
+    ! dummy
     class(PrtModelType) :: this
-    ! -- locals
+    ! locals
     integer(I4B) :: ip
     class(BndType), pointer :: packobj
 
-    ! -- Allocate and read modules attached to model
+    ! Allocate and read modules attached to model
     call this%fmi%fmi_ar(this%ibound)
     if (this%inmip > 0) call this%mip%mip_ar()
 
-    ! -- set up output control
+    ! set up output control
     call this%oc%oc_ar(this%dis, DHNOFLO)
     call this%budget%set_ibudcsv(this%oc%ibudcsv)
 
-    ! -- Package input files now open, so allocate and read
+    ! Package input files now open, so allocate and read
     do ip = 1, this%bndlist%Count()
       packobj => GetBndFromList(this%bndlist, ip)
       select type (packobj)
@@ -252,11 +252,11 @@ subroutine prt_ar(this)
         call packobj%prp_set_pointers(this%ibound, this%mip%izone, &
                                       this%trackctl)
       end select
-      ! -- Read and allocate package
+      ! Read and allocate package
       call packobj%bnd_ar()
     end do
 
-    ! -- Initialize tracking method
+    ! Initialize tracking method
     select type (dis => this%dis)
     type is (DisType)
       call method_dis%init( &
@@ -280,7 +280,7 @@ subroutine prt_ar(this)
       this%method => method_disv
     end select
 
-    ! -- Initialize track output files and reporting options
+    ! Initialize track output files and reporting options
     if (this%oc%itrkout > 0) &
       call this%trackctl%init_track_file(this%oc%itrkout)
     if (this%oc%itrkcsv > 0) &
@@ -297,16 +297,16 @@ end subroutine prt_ar
   !> @brief Read and prepare (calls package read and prepare routines)
   subroutine prt_rp(this)
     use TdisModule, only: readnewdata
-    ! -- dummy
+    ! dummy
     class(PrtModelType) :: this
-    ! -- local
+    ! local
     class(BndType), pointer :: packobj
     integer(I4B) :: ip
 
-    ! -- Check with TDIS on whether or not it is time to RP
+    ! Check with TDIS on whether or not it is time to RP
     if (.not. readnewdata) return
 
-    ! -- Read and prepare
+    ! Read and prepare
     if (this%inoc > 0) call this%oc%oc_rp()
     do ip = 1, this%bndlist%Count()
       packobj => GetBndFromList(this%bndlist, ip)
@@ -316,28 +316,28 @@ end subroutine prt_rp
 
   !> @brief Time step advance (calls package advance subroutines)
   subroutine prt_ad(this)
-    ! -- modules
+    ! modules
     use SimVariablesModule, only: isimcheck, iFailedStepRetry
-    ! -- dummy
+    ! dummy
     class(PrtModelType) :: this
     class(BndType), pointer :: packobj
-    ! -- local
+    ! local
     integer(I4B) :: irestore
     integer(I4B) :: ip, n, i
 
-    ! -- Reset state variable
+    ! Reset state variable
     irestore = 0
     if (iFailedStepRetry > 0) irestore = 1
 
-    ! -- Copy masssto into massstoold
+    ! Copy masssto into massstoold
     do n = 1, this%dis%nodes
       this%massstoold(n) = this%masssto(n)
     end do
 
-    ! -- Advance fmi
+    ! Advance fmi
     call this%fmi%fmi_ad()
 
-    ! -- Advance
+    ! Advance
     do ip = 1, this%bndlist%Count()
       packobj => GetBndFromList(this%bndlist, ip)
       call packobj%bnd_ad()
@@ -346,7 +346,7 @@ subroutine prt_ad(this)
       end if
     end do
     !
-    ! -- Initialize the flowja array.  Flowja is calculated each time,
+    ! Initialize the flowja array.  Flowja is calculated each time,
     !    even if output is suppressed.  (Flowja represents flow of particle
     !    mass and is positive into a cell.  Currently, each particle is assigned
     !    unit mass.)  Flowja is updated continually as particles are tracked
@@ -359,28 +359,28 @@ end subroutine prt_ad
 
   !> @brief Calculate intercell flow (flowja)
   subroutine prt_cq(this, icnvg, isuppress_output)
-    ! -- modules
+    ! modules
     use SparseModule, only: csr_diagsum
     use TdisModule, only: delt
     use PrtPrpModule, only: PrtPrpType
-    ! -- dummy
+    ! dummy
     class(PrtModelType) :: this
     integer(I4B), intent(in) :: icnvg
     integer(I4B), intent(in) :: isuppress_output
-    ! -- local
+    ! local
     integer(I4B) :: i
     integer(I4B) :: ip
     class(BndType), pointer :: packobj
     real(DP) :: tled
 
-    ! -- Flowja is calculated each time, even if output is suppressed.
+    ! Flowja is calculated each time, even if output is suppressed.
     !    Flowja represents flow of particle mass and is positive into a cell.
     !    Currently, each particle is assigned unit mass.
     !
-    ! -- Reciprocal of time step size.
+    ! Reciprocal of time step size.
     tled = DONE / delt
     !
-    ! -- Flowja was updated continually as particles were tracked over the
+    ! Flowja was updated continually as particles were tracked over the
     !    time step.  At this point, flowja contains the net particle mass
     !    exchanged between cells during the time step.  To convert these to
     !    flow rates (particle mass per time), divide by the time step size.
@@ -388,16 +388,16 @@ subroutine prt_cq(this, icnvg, isuppress_output)
       this%flowja(i) = this%flowja(i) * tled
     end do
 
-    ! -- Particle mass storage
+    ! Particle mass storage
     call this%prt_cq_sto()
 
-    ! -- Go through packages and call cq routines. Just a formality.
+    ! Go through packages and call cq routines. Just a formality.
     do ip = 1, this%bndlist%Count()
       packobj => GetBndFromList(this%bndlist, ip)
       call packobj%bnd_cq(this%masssto, this%flowja)
     end do
 
-    ! -- Finalize calculation of flowja by adding face flows to the diagonal.
+    ! Finalize calculation of flowja by adding face flows to the diagonal.
     !    This results in the flow residual being stored in the diagonal
     !    position for each cell.
     call csr_diagsum(this%dis%con%ia, this%flowja)
@@ -405,12 +405,12 @@ end subroutine prt_cq
 
   !> @brief Calculate particle mass storage
   subroutine prt_cq_sto(this)
-    ! -- modules
+    ! modules
     use TdisModule, only: delt
     use PrtPrpModule, only: PrtPrpType
-    ! -- dummy
+    ! dummy
     class(PrtModelType) :: this
-    ! -- local
+    ! local
     integer(I4B) :: ip
     class(BndType), pointer :: packobj
     integer(I4B) :: n
@@ -420,10 +420,10 @@ subroutine prt_cq_sto(this)
     real(DP) :: tled
     real(DP) :: rate
 
-    ! -- Reciprocal of time step size.
+    ! Reciprocal of time step size.
     tled = DONE / delt
 
-    ! -- Particle mass storage rate
+    ! Particle mass storage rate
     do n = 1, this%dis%nodes
       this%masssto(n) = DZERO
       this%ratesto(n) = DZERO
@@ -437,7 +437,7 @@ subroutine prt_cq_sto(this)
           ! this may need to change if istatus flags change
           if ((istatus > 0) .and. (istatus /= 8)) then
             n = packobj%particles%idomain(np, 2)
-            ! -- Each particle currently assigned unit mass
+            ! Each particle currently assigned unit mass
             this%masssto(n) = this%masssto(n) + DONE
           end if
         end do
@@ -458,20 +458,20 @@ end subroutine prt_cq_sto
   !!
   !<
   subroutine prt_bd(this, icnvg, isuppress_output)
-    ! -- modules
+    ! modules
     use TdisModule, only: delt
     use BudgetModule, only: rate_accumulator
-    ! -- dummy
+    ! dummy
     class(PrtModelType) :: this
     integer(I4B), intent(in) :: icnvg
     integer(I4B), intent(in) :: isuppress_output
-    ! -- local
+    ! local
     integer(I4B) :: ip
     class(BndType), pointer :: packobj
     real(DP) :: rin
     real(DP) :: rout
 
-    ! -- Budget routines (start by resetting).  Sole purpose of this section
+    ! Budget routines (start by resetting).  Sole purpose of this section
     !    is to add in and outs to model budget.  All ins and out for a model
     !    should be added here to this%budget.  In a subsequent exchange call,
     !    exchange flows might also be added.
@@ -488,9 +488,9 @@ end subroutine prt_bd
   !> @brief Print and/or save model output
   subroutine prt_ot(this)
     use TdisModule, only: tdis_ot, endofperiod
-    ! -- dummy
+    ! dummy
     class(PrtModelType) :: this
-    ! -- local
+    ! local
     integer(I4B) :: idvsave
     integer(I4B) :: idvprint
     integer(I4B) :: icbcfl
@@ -498,9 +498,9 @@ subroutine prt_ot(this)
     integer(I4B) :: ibudfl
     integer(I4B) :: ipflag
 
-    ! -- Note: particle tracking output is handled elsewhere
+    ! Note: particle tracking output is handled elsewhere
 
-    ! -- Set write and print flags
+    ! Set write and print flags
     idvsave = 0
     idvprint = 0
     icbcfl = 0
@@ -511,21 +511,21 @@ subroutine prt_ot(this)
     if (this%oc%oc_print('BUDGET')) ibudfl = 1
     icbcun = this%oc%oc_save_unit('BUDGET')
 
-    ! -- Override ibudfl and idvprint flags for nonconvergence
+    ! Override ibudfl and idvprint flags for nonconvergence
     !    and end of period
     ibudfl = this%oc%set_print_flag('BUDGET', 1, endofperiod)
     idvprint = this%oc%set_print_flag('CONCENTRATION', 1, endofperiod)
 
-    ! -- Save and print flows
+    ! Save and print flows
     call this%prt_ot_flow(icbcfl, ibudfl, icbcun)
 
-    ! -- Save and print dependent variables
+    ! Save and print dependent variables
     call this%prt_ot_dv(idvsave, idvprint, ipflag)
 
-    ! -- Print budget summaries
+    ! Print budget summaries
     call this%prt_ot_bdsummary(ibudfl, ipflag)
 
-    ! -- Timing Output; if any dependent variables or budgets
+    ! Timing Output; if any dependent variables or budgets
     !    are printed, then ipflag is set to 1.
     if (ipflag == 1) call tdis_ot(this%iout)
   end subroutine prt_ot
@@ -540,27 +540,27 @@ subroutine prt_ot_flow(this, icbcfl, ibudfl, icbcun)
     class(BndType), pointer :: packobj
     integer(I4B) :: ip
 
-    ! -- Save PRT flows
+    ! Save PRT flows
     call this%prt_ot_saveflow(this%dis%nja, this%flowja, icbcfl, icbcun)
     do ip = 1, this%bndlist%Count()
       packobj => GetBndFromList(this%bndlist, ip)
       call packobj%bnd_ot_model_flows(icbcfl=icbcfl, ibudfl=0, icbcun=icbcun)
     end do
 
-    ! -- Save advanced package flows
+    ! Save advanced package flows
     do ip = 1, this%bndlist%Count()
       packobj => GetBndFromList(this%bndlist, ip)
       call packobj%bnd_ot_package_flows(icbcfl=icbcfl, ibudfl=0)
     end do
 
-    ! -- Print PRT flows
+    ! Print PRT flows
     call this%prt_ot_printflow(ibudfl, this%flowja)
     do ip = 1, this%bndlist%Count()
       packobj => GetBndFromList(this%bndlist, ip)
       call packobj%bnd_ot_model_flows(icbcfl=icbcfl, ibudfl=ibudfl, icbcun=0)
     end do
 
-    ! -- Print advanced package flows
+    ! Print advanced package flows
     do ip = 1, this%bndlist%Count()
       packobj => GetBndFromList(this%bndlist, ip)
       call packobj%bnd_ot_package_flows(icbcfl=0, ibudfl=ibudfl)
@@ -569,16 +569,16 @@ end subroutine prt_ot_flow
 
   !> @brief Save intercell flows
   subroutine prt_ot_saveflow(this, nja, flowja, icbcfl, icbcun)
-    ! -- dummy
+    ! dummy
     class(PrtModelType) :: this
     integer(I4B), intent(in) :: nja
     real(DP), dimension(nja), intent(in) :: flowja
     integer(I4B), intent(in) :: icbcfl
     integer(I4B), intent(in) :: icbcun
-    ! -- local
+    ! local
     integer(I4B) :: ibinun
 
-    ! -- Set unit number for binary output
+    ! Set unit number for binary output
     if (this%ipakcb < 0) then
       ibinun = icbcun
     elseif (this%ipakcb == 0) then
@@ -588,7 +588,7 @@ subroutine prt_ot_saveflow(this, nja, flowja, icbcfl, icbcun)
     end if
     if (icbcfl == 0) ibinun = 0
 
-    ! -- Write the face flows if requested
+    ! Write the face flows if requested
     if (ibinun /= 0) then
       call this%dis%record_connection_array(flowja, ibinun, this%iout)
     end if
@@ -596,24 +596,24 @@ end subroutine prt_ot_saveflow
 
   !> @brief Print intercell flows
   subroutine prt_ot_printflow(this, ibudfl, flowja)
-    ! -- modules
+    ! modules
     use TdisModule, only: kper, kstp
     use ConstantsModule, only: LENBIGLINE
-    ! -- dummy
+    ! dummy
     class(PrtModelType) :: this
     integer(I4B), intent(in) :: ibudfl
     real(DP), intent(inout), dimension(:) :: flowja
-    ! -- local
+    ! local
     character(len=LENBIGLINE) :: line
     character(len=30) :: tempstr
     integer(I4B) :: n, ipos, m
     real(DP) :: qnm
-    ! -- formats
+    ! formats
     character(len=*), parameter :: fmtiprflow = &
                 "(/,4x,'CALCULATED INTERCELL FLOW &
                 &FOR PERIOD ', i0, ' STEP ', i0)"
 
-    ! -- Write flowja to list file if requested
+    ! Write flowja to list file if requested
     if (ibudfl /= 0 .and. this%iprflow > 0) then
       write (this%iout, fmtiprflow) kper, kstp
       do n = 1, this%dis%nodes
@@ -635,75 +635,75 @@ end subroutine prt_ot_printflow
 
   !> @brief Print dependent variables
   subroutine prt_ot_dv(this, idvsave, idvprint, ipflag)
-    ! -- dummy
+    ! dummy
     class(PrtModelType) :: this
     integer(I4B), intent(in) :: idvsave
     integer(I4B), intent(in) :: idvprint
     integer(I4B), intent(inout) :: ipflag
-    ! -- local
+    ! local
     class(BndType), pointer :: packobj
     integer(I4B) :: ip
 
-    ! -- Print advanced package dependent variables
+    ! Print advanced package dependent variables
     do ip = 1, this%bndlist%Count()
       packobj => GetBndFromList(this%bndlist, ip)
       call packobj%bnd_ot_dv(idvsave, idvprint)
     end do
 
-    ! -- save head and print head
+    ! save head and print head
     call this%oc%oc_ot(ipflag)
   end subroutine prt_ot_dv
 
   !> @brief Print budget summary
   subroutine prt_ot_bdsummary(this, ibudfl, ipflag)
-    ! -- modules
+    ! modules
     use TdisModule, only: kstp, kper, totim, delt
-    ! -- dummy
+    ! dummy
     class(PrtModelType) :: this
     integer(I4B), intent(in) :: ibudfl
     integer(I4B), intent(inout) :: ipflag
-    ! -- local
+    ! local
     class(BndType), pointer :: packobj
     integer(I4B) :: ip
 
-    ! -- Package budget summary
+    ! Package budget summary
     do ip = 1, this%bndlist%Count()
       packobj => GetBndFromList(this%bndlist, ip)
       call packobj%bnd_ot_bdsummary(kstp, kper, this%iout, ibudfl)
     end do
 
-    ! -- model budget summary
+    ! model budget summary
     call this%budget%finalize_step(delt)
     if (ibudfl /= 0) then
       ipflag = 1
-      ! -- model budget summary
+      ! model budget summary
       call this%budget%budget_ot(kstp, kper, this%iout)
     end if
 
-    ! -- Write to budget csv
+    ! Write to budget csv
     call this%budget%writecsv(totim)
   end subroutine prt_ot_bdsummary
 
   !> @brief Deallocate
   subroutine prt_da(this)
-    ! -- modules
+    ! modules
     use MemoryManagerModule, only: mem_deallocate
     use MemoryManagerExtModule, only: memorystore_remove
     use SimVariablesModule, only: idm_context
     use MethodPoolModule, only: destroy_method_pool
     use MethodCellPoolModule, only: destroy_method_cell_pool
     use MethodSubcellPoolModule, only: destroy_method_subcell_pool
-    ! -- dummy
+    ! dummy
     class(PrtModelType) :: this
-    ! -- local
+    ! local
     integer(I4B) :: ip
     class(BndType), pointer :: packobj
 
-    ! -- Deallocate idm memory
+    ! Deallocate idm memory
     call memorystore_remove(this%name, 'NAM', idm_context)
     call memorystore_remove(component=this%name, context=idm_context)
 
-    ! -- Internal packages
+    ! Internal packages
     call this%dis%dis_da()
     call this%fmi%fmi_da()
     call this%mip%mip_da()
@@ -715,19 +715,19 @@ subroutine prt_da(this)
     deallocate (this%budget)
     deallocate (this%oc)
 
-    ! -- Method objects
+    ! Method objects
     call destroy_method_subcell_pool()
     call destroy_method_cell_pool()
     call destroy_method_pool()
 
-    ! -- Boundary packages
+    ! Boundary packages
     do ip = 1, this%bndlist%Count()
       packobj => GetBndFromList(this%bndlist, ip)
       call packobj%bnd_da()
       deallocate (packobj)
     end do
 
-    ! -- Scalars
+    ! Scalars
     call mem_deallocate(this%infmi)
     call mem_deallocate(this%inmip)
     call mem_deallocate(this%inadv)
@@ -737,28 +737,26 @@ subroutine prt_da(this)
     call mem_deallocate(this%inmvt)
     call mem_deallocate(this%inoc)
 
-    ! -- Arrays
+    ! Arrays
     call mem_deallocate(this%masssto)
     call mem_deallocate(this%massstoold)
     call mem_deallocate(this%ratesto)
 
-    ! -- Track file control
     deallocate (this%trackctl)
 
-    ! -- Parent type
     call this%NumericalModelType%model_da()
   end subroutine prt_da
 
-  !> @brief Allocate memory for non-allocatable members
+  !> @brief Allocate memory for scalars
   subroutine allocate_scalars(this, modelname)
-    ! -- dummy
+    ! dummy
     class(PrtModelType) :: this
     character(len=*), intent(in) :: modelname
 
-    ! -- allocate members from parent class
+    ! allocate members from parent class
     call this%NumericalModelType%allocate_scalars(modelname)
 
-    ! -- allocate members that are part of model class
+    ! allocate members that are part of model class
     call mem_allocate(this%infmi, 'INFMI', this%memoryPath)
     call mem_allocate(this%inmip, 'INMIP', this%memoryPath)
     call mem_allocate(this%inmvt, 'INMVT', this%memoryPath)
@@ -784,18 +782,18 @@ subroutine allocate_arrays(this)
     class(PrtModelType) :: this
     integer(I4B) :: n
 
-    ! -- Allocate arrays in parent type
+    ! Allocate arrays in parent type
     this%nja = this%dis%nja
     call this%NumericalModelType%allocate_arrays()
 
-    ! -- Allocate and initialize arrays
+    ! Allocate and initialize arrays
     call mem_allocate(this%masssto, this%dis%nodes, &
                       'MASSSTO', this%memoryPath)
     call mem_allocate(this%massstoold, this%dis%nodes, &
                       'MASSSTOOLD', this%memoryPath)
     call mem_allocate(this%ratesto, this%dis%nodes, &
                       'RATESTO', this%memoryPath)
-    ! -- explicit model, so these must be manually allocated
+    ! explicit model, so these must be manually allocated
     call mem_allocate(this%x, this%dis%nodes, 'X', this%memoryPath)
     call mem_allocate(this%rhs, this%dis%nodes, 'RHS', this%memoryPath)
     call mem_allocate(this%ibound, this%dis%nodes, 'IBOUND', this%memoryPath)
@@ -812,11 +810,11 @@ end subroutine allocate_arrays
   !> @brief Create boundary condition packages for this model
   subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, mempath, &
                             inunit, iout)
-    ! -- modules
+    ! modules
     use ConstantsModule, only: LINELENGTH
     use PrtPrpModule, only: prp_create
     use ApiModule, only: api_create
-    ! -- dummy
+    ! dummy
     class(PrtModelType) :: this
     character(len=*), intent(in) :: filtyp
     character(len=LINELENGTH) :: errmsg
@@ -826,12 +824,12 @@ subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, mempath, &
     character(len=*), intent(in) :: mempath
     integer(I4B), intent(in) :: inunit
     integer(I4B), intent(in) :: iout
-    ! -- local
+    ! local
     class(BndType), pointer :: packobj
     class(BndType), pointer :: packobj2
     integer(I4B) :: ip
 
-    ! -- This part creates the package object
+    ! This part creates the package object
     select case (filtyp)
     case ('PRP6')
       call prp_create(packobj, ipakid, ipaknum, inunit, iout, &
@@ -844,9 +842,9 @@ subroutine package_create(this, filtyp, ipakid, ipaknum, pakname, mempath, &
       call store_error(errmsg, terminate=.TRUE.)
     end select
 
-    ! -- Packages is the bndlist that is associated with the parent model
-    ! -- The following statement puts a pointer to this package in the ipakid
-    ! -- position of packages.
+    ! Packages is the bndlist that is associated with the parent model
+    ! The following statement puts a pointer to this package in the ipakid
+    ! position of packages.
     do ip = 1, this%bndlist%Count()
       packobj2 => GetBndFromList(this%bndlist, ip)
       if (packobj2%packName == pakname) then
@@ -860,13 +858,13 @@ end subroutine package_create
 
   !> @brief Check to make sure required input files have been specified
   subroutine ftype_check(this, indis)
-    ! -- dummy
+    ! dummy
     class(PrtModelType) :: this
     integer(I4B), intent(in) :: indis
-    ! -- local
+    ! local
     character(len=LINELENGTH) :: errmsg
 
-    ! -- Check for DIS(u) and MIP. Stop if not present.
+    ! Check for DIS(u) and MIP. Stop if not present.
     if (indis == 0) then
       write (errmsg, '(1x,a)') &
         'Discretization (DIS6, DISV6, or DISU6) package not specified.'
@@ -887,31 +885,30 @@ end subroutine ftype_check
 
   !> @brief Solve the model
   subroutine prt_solve(this)
-    ! -- modules
-    use TdisModule, only: kper, kstp, totimc, nper, nstp, delt
+    ! modules
+    use TdisModule, only: kper, kstp, totimc, delt, endofsimulation
     use PrtPrpModule, only: PrtPrpType
-    ! -- dummy variables
+    ! dummy variables
     class(PrtModelType) :: this
-    ! -- local variables
+    ! local variables
     integer(I4B) :: np, ip
     class(BndType), pointer :: packobj
     type(ParticleType), pointer :: particle
     real(DP) :: tmax
     integer(I4B) :: iprp
 
-    ! -- Initialize particle
     call create_particle(particle)
 
-    ! -- Loop over PRP packages
+    ! Apply tracking solution to PRP packages
     iprp = 0
     do ip = 1, this%bndlist%Count()
       packobj => GetBndFromList(this%bndlist, ip)
       select type (packobj)
       type is (PrtPrpType)
-        ! -- Update PRP index
+        ! Update PRP index
         iprp = iprp + 1
 
-        ! -- Initialize PRP-specific track files
+        ! Initialize PRP track files
         if (kper == 1 .and. kstp == 1) then
           if (packobj%itrkout > 0) then
             call this%trackctl%init_track_file( &
@@ -926,15 +923,18 @@ subroutine prt_solve(this)
           end if
         end if
 
-        ! -- Loop over particles in package
+        ! Track particles
         do np = 1, packobj%nparticles
-          ! -- Load particle from storage
-          call particle%load_particle(packobj%particles, &
-                                      this%id, iprp, np)
-
-          ! -- If particle is permanently unreleased, record its initial/terminal state
-          if (particle%istatus == 8) &
-            call this%method%save(particle, reason=3)
+          ! Load particle from storage
+          call packobj%particles%get(particle, this%id, iprp, np)
+
+          ! If particle is permanently unreleased, record
+          ! unreleased state if first time step and cycle
+          if (particle%istatus == 8) then
+            if (kper == 1 .and. kstp == 1) &
+              call this%method%save(particle, reason=3)
+            cycle
+          end if
 
           ! If particle is inactive or not yet to be released, cycle
           if (particle%istatus > 1) cycle
@@ -948,9 +948,7 @@ subroutine prt_solve(this)
           ! stop time, whichever comes first, unless it's the final
           ! time step and the extend option is on, in which case
           ! it's just the particle stop time.
-          if (nper == kper .and. &
-              nstp(kper) == kstp .and. &
-              particle%iextend > 0) then
+          if (endofsimulation .and. particle%iextend > 0) then
             tmax = particle%tstop
           else
             tmax = min(totimc + delt, particle%tstop)
@@ -963,23 +961,30 @@ subroutine prt_solve(this)
           particle%icp = 0
           particle%izp = 0
 
+          ! Terminate particles still active at end of simulation
+          if (endofsimulation .and. &
+              particle%iextend == 0 .and. &
+              particle%istatus < 2) then
+            particle%istatus = 10
+            call this%method%save(particle, reason=3)
+          end if
+
           ! Update particle storage
-          call packobj%particles%save_particle(particle, np)
+          call packobj%particles%put(particle, np)
         end do
       end select
     end do
 
-    ! -- Deallocate particle
     deallocate (particle)
   end subroutine prt_solve
 
   !> @brief Source package info and begin to process
   subroutine create_bndpkgs(this, bndpkgs, pkgtypes, pkgnames, &
                             mempaths, inunits)
-    ! -- modules
+    ! modules
     use ConstantsModule, only: LINELENGTH, LENPACKAGENAME
     use CharacterStringModule, only: CharacterStringType
-    ! -- dummy
+    ! dummy
     class(PrtModelType) :: this
     integer(I4B), dimension(:), allocatable, intent(inout) :: bndpkgs
     type(CharacterStringType), dimension(:), contiguous, &
@@ -990,7 +995,7 @@ subroutine create_bndpkgs(this, bndpkgs, pkgtypes, pkgnames, &
       pointer, intent(inout) :: mempaths
     integer(I4B), dimension(:), contiguous, &
       pointer, intent(inout) :: inunits
-    ! -- local
+    ! local
     integer(I4B) :: ipakid, ipaknum
     character(len=LENFTYPE) :: pkgtype, bndptype
     character(len=LENPACKAGENAME) :: pkgname
@@ -1000,7 +1005,7 @@ subroutine create_bndpkgs(this, bndpkgs, pkgtypes, pkgnames, &
 
     if (allocated(bndpkgs)) then
       !
-      ! -- create stress packages
+      ! create stress packages
       ipakid = 1
       bndptype = ''
       do n = 1, size(bndpkgs)
@@ -1021,7 +1026,7 @@ subroutine create_bndpkgs(this, bndpkgs, pkgtypes, pkgnames, &
         ipaknum = ipaknum + 1
       end do
       !
-      ! -- cleanup
+      ! cleanup
       deallocate (bndpkgs)
     end if
 
@@ -1029,7 +1034,7 @@ end subroutine create_bndpkgs
 
   !> @brief Source package info and begin to process
   subroutine create_packages(this)
-    ! -- modules
+    ! modules
     use ConstantsModule, only: LINELENGTH, LENPACKAGENAME
     use CharacterStringModule, only: CharacterStringType
     use ArrayHandlersModule, only: expandarray
@@ -1043,9 +1048,9 @@ subroutine create_packages(this)
     use PrtMipModule, only: mip_cr
     use PrtFmiModule, only: fmi_cr
     use PrtOcModule, only: oc_cr
-    ! -- dummy
+    ! dummy
     class(PrtModelType) :: this
-    ! -- local
+    ! local
     type(CharacterStringType), dimension(:), contiguous, &
       pointer :: pkgtypes => null()
     type(CharacterStringType), dimension(:), contiguous, &
@@ -1064,10 +1069,10 @@ subroutine create_packages(this)
     integer(I4B) :: indis = 0 ! DIS enabled flag
     character(len=LENMEMPATH) :: mempathmip = ''
 
-    ! -- set input memory paths, input/model and input/model/namfile
+    ! set input memory paths, input/model and input/model/namfile
     model_mempath = create_mem_path(component=this%name, context=idm_context)
 
-    ! -- set pointers to model path package info
+    ! set pointers to model path package info
     call mem_setptr(pkgtypes, 'PKGTYPES', model_mempath)
     call mem_setptr(pkgnames, 'PKGNAMES', model_mempath)
     call mem_setptr(mempaths, 'MEMPATHS', model_mempath)
@@ -1080,7 +1085,7 @@ subroutine create_packages(this)
       mempath = mempaths(n)
       inunit => inunits(n)
 
-      ! -- create dis package first as it is a prerequisite for other packages
+      ! create dis package first as it is a prerequisite for other packages
       select case (pkgtype)
       case ('DIS6')
         indis = 1
@@ -1106,23 +1111,23 @@ subroutine create_packages(this)
       end select
     end do
 
-    ! -- Create budget manager
+    ! Create budget manager
     call budget_cr(this%budget, this%name)
 
-    ! -- Create tracking method pools
+    ! Create tracking method pools
     call create_method_pool()
     call create_method_cell_pool()
     call create_method_subcell_pool()
 
-    ! -- Create packages that are tied directly to model
+    ! Create packages that are tied directly to model
     call mip_cr(this%mip, this%name, mempathmip, this%inmip, this%iout, this%dis)
     call fmi_cr(this%fmi, this%name, this%infmi, this%iout)
     call oc_cr(this%oc, this%name, this%inoc, this%iout)
 
-    ! -- Check to make sure that required ftype's have been specified
+    ! Check to make sure that required ftype's have been specified
     call this%ftype_check(indis)
 
-    ! -- Create boundary packages
+    ! Create boundary packages
     call this%create_bndpkgs(bndpkgs, pkgtypes, pkgnames, mempaths, inunits)
   end subroutine create_packages
 
diff --git a/src/Solution/ParticleTracker/Method.f90 b/src/Solution/ParticleTracker/Method.f90
index 4035a0ad8a1..797361e7132 100644
--- a/src/Solution/ParticleTracker/Method.f90
+++ b/src/Solution/ParticleTracker/Method.f90
@@ -191,15 +191,18 @@ end subroutine save
   !! tracking the particle or terminate it, as well as whether to
   !! record any output data as per selected reporting conditions.
   !<
-  subroutine check(this, particle, cell_defn)
+  subroutine check(this, particle, cell_defn, tmax)
     ! modules
-    use TdisModule, only: endofsimulation, totim
+    use TdisModule, only: endofsimulation, totimc, totim
     ! dummy
     class(MethodType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
     type(CellDefnType), pointer, intent(inout) :: cell_defn
+    real(DP), intent(in) :: tmax
     ! local
     logical(LGP) :: dry_cell, dry_particle, no_exit_face, stop_zone, weak_sink
+    integer(I4B) :: i
+    real(DP) :: t, ttrackmax
 
     dry_cell = this%fmi%ibdgwfsat0(cell_defn%icell) == 0
     dry_particle = particle%z > cell_defn%top
@@ -235,6 +238,8 @@ subroutine check(this, particle, cell_defn)
 
     if (dry_cell) then
       if (particle%idrymeth == 0) then
+        ! drop to cell bottom. handled by pass
+        ! to bottom method, nothing to do here
         no_exit_face = .false.
       else if (particle%idrymeth == 1) then
         ! stop
@@ -244,33 +249,91 @@ subroutine check(this, particle, cell_defn)
         return
       else if (particle%idrymeth == 2) then
         ! stay
-        no_exit_face = .false.
         particle%advancing = .false.
+        no_exit_face = .false.
+
+        ! we might report tracking times
+        ! out of order here, but we want
+        ! the particle termination event
+        ! (if this is the last time step)
+        ! to have the maximum tracking t,
+        ! so we need to keep tabs on it.
+        ttrackmax = totim
+
+        ! update tracking time to time
+        ! step end time and save record
         particle%ttrack = totim
+        call this%save(particle, reason=2)
+
+        ! record user tracking times
+        call this%tracktimes%advance()
+        if (this%tracktimes%any()) then
+          do i = this%tracktimes%selection(1), this%tracktimes%selection(2)
+            t = this%tracktimes%times(i)
+            if (t < totimc) cycle
+            if (t >= tmax) exit
+            particle%ttrack = t
+            call this%save(particle, reason=5)
+            if (t > ttrackmax) ttrackmax = t
+          end do
+        end if
+
         ! terminate if last period/step
         if (endofsimulation) then
           particle%istatus = 5
+          particle%ttrack = ttrackmax
           call this%save(particle, reason=3)
           return
         end if
-        call this%save(particle, reason=2)
       end if
     else if (dry_particle .and. this%name /= "passtobottom") then
-      ! dry particle
       if (particle%idrymeth == 0) then
         ! drop to water table
         particle%z = cell_defn%top
         call this%save(particle, reason=1)
       else if (particle%idrymeth == 1) then
-        ! terminate
+        ! stop
         particle%advancing = .false.
         particle%istatus = 7
         call this%save(particle, reason=3)
         return
       else if (particle%idrymeth == 2) then
         ! stay
-        no_exit_face = .false.
         particle%advancing = .false.
+        no_exit_face = .false.
+
+        ! we might report tracking times
+        ! out of order here, but we want
+        ! the particle termination event
+        ! (if this is the last time step)
+        ! to have the maximum tracking t,
+        ! so we need to keep tabs on it.
+        ttrackmax = totim
+
+        ! update tracking time to time
+        ! step end time and save record
+        particle%ttrack = totim
+        call this%save(particle, reason=2)
+
+        ! record user tracking times
+        call this%tracktimes%advance()
+        if (this%tracktimes%any()) then
+          do i = this%tracktimes%selection(1), this%tracktimes%selection(2)
+            t = this%tracktimes%times(i)
+            if (t < totimc) cycle
+            if (t >= tmax) exit
+            particle%ttrack = t
+            call this%save(particle, reason=5)
+            if (t > ttrackmax) ttrackmax = t
+          end do
+        end if
+
+        ! terminate if last period/step
+        if (endofsimulation) then
+          particle%istatus = 5
+          particle%ttrack = ttrackmax
+          call this%save(particle, reason=3)
+        end if
         return
       end if
     end if
diff --git a/src/Solution/ParticleTracker/MethodCellPassToBot.f90 b/src/Solution/ParticleTracker/MethodCellPassToBot.f90
index 433b3ca9650..bad2d8a5fa4 100644
--- a/src/Solution/ParticleTracker/MethodCellPassToBot.f90
+++ b/src/Solution/ParticleTracker/MethodCellPassToBot.f90
@@ -5,6 +5,8 @@ module MethodCellPassToBotModule
   use CellDefnModule, only: CellDefnType, create_defn
   use PrtFmiModule, only: PrtFmiType
   use BaseDisModule, only: DisBaseType
+  use DisModule, only: DisType
+  use DisvModule, only: DisvType
   use ParticleModule, only: ParticleType
   use CellModule, only: CellType
   use SubcellModule, only: SubcellType
@@ -45,16 +47,31 @@ subroutine apply_ptb(this, particle, tmax)
     class(MethodCellPassToBotType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
     real(DP), intent(in) :: tmax
+    ! local
+    integer(I4B) :: nlay
 
     ! Check termination/reporting conditions
-    call this%check(particle, this%cell%defn)
+    call this%check(particle, this%cell%defn, tmax)
     if (.not. particle%advancing) return
 
     ! Pass to bottom face
     particle%z = this%cell%defn%bot
     particle%iboundary(2) = this%cell%defn%npolyverts + 2
 
-    ! Save datum
+    ! Terminate if in bottom layer
+    select type (dis => this%fmi%dis)
+    type is (DisType)
+      nlay = dis%nlay
+    type is (DisvType)
+      nlay = dis%nlay
+    end select
+    if (particle%ilay == nlay) then
+      particle%advancing = .false.
+      particle%istatus = 5
+      call this%save(particle, reason=3)
+    end if
+
+    ! Save record
     call this%save(particle, reason=1)
   end subroutine apply_ptb
 
diff --git a/src/Solution/ParticleTracker/MethodCellPollock.f90 b/src/Solution/ParticleTracker/MethodCellPollock.f90
index 86b1d09f4ca..06790e8f7f5 100644
--- a/src/Solution/ParticleTracker/MethodCellPollock.f90
+++ b/src/Solution/ParticleTracker/MethodCellPollock.f90
@@ -130,7 +130,7 @@ subroutine apply_mcp(this, particle, tmax)
     select type (cell => this%cell)
     type is (CellRectType)
       ! Check termination/reporting conditions
-      call this%check(particle, cell%defn)
+      call this%check(particle, cell%defn, tmax)
       if (.not. particle%advancing) return
 
       ! Transform model coordinates to local cell coordinates
diff --git a/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 b/src/Solution/ParticleTracker/MethodCellPollockQuad.f90
index c5f61f2c221..5bd749e192d 100644
--- a/src/Solution/ParticleTracker/MethodCellPollockQuad.f90
+++ b/src/Solution/ParticleTracker/MethodCellPollockQuad.f90
@@ -209,7 +209,7 @@ subroutine apply_mcpq(this, particle, tmax)
     select type (cell => this%cell)
     type is (CellRectQuadType)
       ! Check termination/reporting conditions
-      call this%check(particle, cell%defn)
+      call this%check(particle, cell%defn, tmax)
       if (.not. particle%advancing) return
 
       ! Transform model coordinates to local cell coordinates
diff --git a/src/Solution/ParticleTracker/MethodCellTernary.f90 b/src/Solution/ParticleTracker/MethodCellTernary.f90
index d6447db6638..03a401173c2 100644
--- a/src/Solution/ParticleTracker/MethodCellTernary.f90
+++ b/src/Solution/ParticleTracker/MethodCellTernary.f90
@@ -165,7 +165,7 @@ subroutine apply_mct(this, particle, tmax)
     select type (cell => this%cell)
     type is (CellPolyType)
       ! Check termination/reporting conditions
-      call this%check(particle, this%cell%defn)
+      call this%check(particle, this%cell%defn, tmax)
       if (.not. particle%advancing) return
 
       ! Number of vertices
diff --git a/src/Solution/ParticleTracker/Particle.f90 b/src/Solution/ParticleTracker/Particle.f90
index 0904cd89e6a..2506bc307aa 100644
--- a/src/Solution/ParticleTracker/Particle.f90
+++ b/src/Solution/ParticleTracker/Particle.f90
@@ -69,7 +69,6 @@ module ParticleModule
     integer(I4B), public :: iextend !< whether to extend tracking beyond the end of the simulation
   contains
     procedure, public :: get_model_coords
-    procedure, public :: load_particle
     procedure, public :: transform => transform_coords
     procedure, public :: reset_transform
   end type ParticleType
@@ -108,7 +107,8 @@ module ParticleModule
     procedure, public :: deallocate
     procedure, public :: num_stored
     procedure, public :: resize
-    procedure, public :: save_particle
+    procedure, public :: get
+    procedure, public :: put
   end type ParticleStoreType
 
 contains
@@ -224,50 +224,50 @@ end subroutine resize
   !! This routine is used to initialize a particle for tracking.
   !! The advancing flag and coordinate transformation are reset.
   !<
-  subroutine load_particle(this, store, imdl, iprp, ip)
-    class(ParticleType), intent(inout) :: this !< particle
-    type(ParticleStoreType), intent(in) :: store !< particle storage
+  subroutine get(this, particle, imdl, iprp, ip)
+    class(ParticleStoreType), intent(inout) :: this !< particle store
+    class(ParticleType), intent(inout) :: particle !< particle
     integer(I4B), intent(in) :: imdl !< index of model particle originated in
     integer(I4B), intent(in) :: iprp !< index of particle release package particle originated in
     integer(I4B), intent(in) :: ip !< index into the particle list
 
-    call this%reset_transform()
-    this%imdl = imdl
-    this%iprp = iprp
-    this%irpt = store%irpt(ip)
-    this%ip = ip
-    this%name = store%name(ip)
-    this%istopweaksink = store%istopweaksink(ip)
-    this%istopzone = store%istopzone(ip)
-    this%idrymeth = store%idrymeth(ip)
-    this%icp = 0
-    this%icu = store%icu(ip)
-    this%ilay = store%ilay(ip)
-    this%izone = store%izone(ip)
-    this%izp = store%izp(ip)
-    this%istatus = store%istatus(ip)
-    this%x = store%x(ip)
-    this%y = store%y(ip)
-    this%z = store%z(ip)
-    this%trelease = store%trelease(ip)
-    this%tstop = store%tstop(ip)
-    this%ttrack = store%ttrack(ip)
-    this%advancing = .true.
-    this%idomain(1:levelmax) = &
-      store%idomain(ip, 1:levelmax)
-    this%idomain(1) = imdl
-    this%iboundary(1:levelmax) = &
-      store%iboundary(ip, 1:levelmax)
-    this%ifrctrn = store%ifrctrn(ip)
-    this%iexmeth = store%iexmeth(ip)
-    this%extol = store%extol(ip)
-    this%iextend = store%extend(ip)
-  end subroutine load_particle
+    call particle%reset_transform()
+    particle%imdl = imdl
+    particle%iprp = iprp
+    particle%irpt = this%irpt(ip)
+    particle%ip = ip
+    particle%name = this%name(ip)
+    particle%istopweaksink = this%istopweaksink(ip)
+    particle%istopzone = this%istopzone(ip)
+    particle%idrymeth = this%idrymeth(ip)
+    particle%icp = 0
+    particle%icu = this%icu(ip)
+    particle%ilay = this%ilay(ip)
+    particle%izone = this%izone(ip)
+    particle%izp = this%izp(ip)
+    particle%istatus = this%istatus(ip)
+    particle%x = this%x(ip)
+    particle%y = this%y(ip)
+    particle%z = this%z(ip)
+    particle%trelease = this%trelease(ip)
+    particle%tstop = this%tstop(ip)
+    particle%ttrack = this%ttrack(ip)
+    particle%advancing = .true.
+    particle%idomain(1:levelmax) = &
+      this%idomain(ip, 1:levelmax)
+    particle%idomain(1) = imdl
+    particle%iboundary(1:levelmax) = &
+      this%iboundary(ip, 1:levelmax)
+    particle%ifrctrn = this%ifrctrn(ip)
+    particle%iexmeth = this%iexmeth(ip)
+    particle%extol = this%extol(ip)
+    particle%iextend = this%extend(ip)
+  end subroutine get
 
   !> @brief Save a particle's state to the particle store
-  subroutine save_particle(this, particle, ip)
+  subroutine put(this, particle, ip)
     class(ParticleStoreType), intent(inout) :: this !< particle storage
-    type(ParticleType), intent(in) :: particle !< particle
+    class(ParticleType), intent(in) :: particle !< particle
     integer(I4B), intent(in) :: ip !< particle index
 
     this%imdl(ip) = particle%imdl
@@ -300,7 +300,7 @@ subroutine save_particle(this, particle, ip)
     this%iexmeth(ip) = particle%iexmeth
     this%extol(ip) = particle%extol
     this%extend(ip) = particle%iextend
-  end subroutine save_particle
+  end subroutine put
 
   !> @brief Transform particle coordinates.
   !!