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 f29364335ea..5a64e1610b8 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 e4fe767b0a4..461954bda8e 100644
--- a/autotest/test_prt_dry.py
+++ b/autotest/test_prt_dry.py
@@ -369,10 +369,10 @@ def check_output(idx, test, snapshot):
     hds_path = gwf_ws / hds_file
     trackcsv_file = prt_name + ".csv"
     trackcsv_path = prt_ws / trackcsv_file
-    mf6pathlines = pd.read_csv(trackcsv_path)
-    startpts = mf6pathlines[mf6pathlines.ireason == 0]
+    pls = pd.read_csv(trackcsv_path)
+    strtpts = pls[pls.ireason == 0]
 
-    assert snapshot == mf6pathlines.drop(["name", "icell"], axis=1).round(3).to_records(
+    assert snapshot == pls.drop(["name", "icell"], axis=1).round(3).to_records(
         index=False
     )
 
@@ -438,7 +438,7 @@ def plot_pathlines_and_timeseries(
             plt.scatter(pathlines["x"], pathlines["y"])
             mm.plot_pathline(pathlines, layer="all", colors=["blue"], lw=0.75)
             ax.set_title(plottitle, fontsize=12)
-            ax.scatter(startpts.x, startpts.y)
+            ax.scatter(strtpts.x, strtpts.y)
             from shapely.geometry import Polygon
 
             cellids = [105, 85, 125, 104, 106]
@@ -468,7 +468,7 @@ def plot_pathlines_and_timeseries(
         ibd[ilay, irow, icol] = 2
         ibd = np.ma.masked_equal(ibd, 0)
 
-        plot_pathlines_and_timeseries(ax, gwf.modelgrid, ibd, mf6pathlines, None, name)
+        plot_pathlines_and_timeseries(ax, gwf.modelgrid, ibd, pls, None, name)
 
     plot_3d = False
     if plot_3d:
@@ -482,7 +482,7 @@ def plot_pathlines_and_timeseries(
         vert_exag = 1
         vtk = Vtk(model=gwf, binary=False, vertical_exageration=vert_exag, smooth=False)
         vtk.add_model(gwf)
-        vtk.add_pathline_points(mf6pathlines.to_records(index=False))
+        vtk.add_pathline_points(pls.to_records(index=False))
         gwf_mesh, prt_mesh = vtk.to_pyvista()
         axes = pv.Axes(show_actor=False, actor_scale=2.0, line_width=5)
         p = pv.Plotter(window_size=[700, 700])
diff --git a/autotest/test_prt_disv1.py b/autotest/test_prt_newton_disv1.py
similarity index 98%
rename from autotest/test_prt_disv1.py
rename to autotest/test_prt_newton_disv1.py
index 98d8bbf65d8..3e5eeec7b30 100644
--- a/autotest/test_prt_disv1.py
+++ b/autotest/test_prt_newton_disv1.py
@@ -6,13 +6,7 @@
 refined cells, instead of the new ternary method
 which applies more generally to polygonal cells.
 
-The simulation includes a single stress period
-with multiple time steps. This serves to test
-whether PRT properly solves trajectories over
-"internal" time steps, i.e. within the step's
-slice of simulation time, as well as extending
-tracking to termination or particle stop times
-during the simulation's final time step.
+The simulation uses the Newton formulation.
 
 Several cases are provided:
     - default: No user-specified tracking times, MP7 in pathline mode.
diff --git a/src/Model/ParticleTracking/prt-prp.f90 b/src/Model/ParticleTracking/prt-prp.f90
index 936f4dfde68..27a49edeb26 100644
--- a/src/Model/ParticleTracking/prt-prp.f90
+++ b/src/Model/ParticleTracking/prt-prp.f90
@@ -443,7 +443,7 @@ subroutine initialize_particle(this, particle, ip, trelease)
     real(DP), intent(in) :: trelease !< release time
     ! local
     integer(I4B) :: irow, icol, ilay, icpl
-    integer(I4B) :: ic, icu
+    integer(I4B) :: ic, icu, ic_old
     real(DP) :: x, y, z
     real(DP) :: top, bot, hds
 
@@ -461,7 +461,7 @@ subroutine initialize_particle(this, particle, ip, trelease)
     particle%irpt = ip
     particle%istopweaksink = this%istopweaksink
     particle%istopzone = this%istopzone
-    particle%idry = this%idry
+    particle%idrymeth = this%idry
     particle%icu = icu
 
     select type (dis => this%dis)
@@ -478,9 +478,14 @@ subroutine initialize_particle(this, particle, ip, trelease)
     ! to the top-most active cell beneath it if drape is
     ! enabled, or else terminate permanently unreleased.
     if (this%ibound(ic) == 0) then
-      if (this%idrape > 0) &
+      ic_old = ic
+      if (this%idrape > 0) then
         call this%dis%highest_active(ic, this%ibound)
-      if (this%ibound(ic) == 0) particle%istatus = 8
+        if (ic == ic_old .or. this%ibound(ic) == 0) &
+          particle%istatus = 8
+      else
+        particle%istatus = 8
+      end if
     end if
 
     ! Load coordinates and transform if needed
@@ -928,7 +933,7 @@ subroutine prp_read_packagedata(this)
 
         if (n < 1 .or. n > this%nreleasepoints) then
           write (errmsg, '(a,1x,i0,a)') &
-            'Release point number must be greater than 0 and less than ', &
+            'Release point number must be greater than 0 and less than', &
             'or equal to', this%nreleasepoints, '.'
           call store_error(errmsg)
           cycle
diff --git a/src/Model/ParticleTracking/prt.f90 b/src/Model/ParticleTracking/prt.f90
index e47f8ff6a0b..2a9caf88e32 100644
--- a/src/Model/ParticleTracking/prt.f90
+++ b/src/Model/ParticleTracking/prt.f90
@@ -932,7 +932,7 @@ subroutine prt_solve(this)
 
           ! -- If particle is permanently unreleased, record its initial/terminal state
           if (particle%istatus == 8) &
-            call this%method%save(particle, reason=3) ! reason=3: termination
+            call this%method%save(particle, reason=3)
 
           ! If particle is inactive or not yet to be released, cycle
           if (particle%istatus > 1) cycle
@@ -940,7 +940,7 @@ subroutine prt_solve(this)
           ! If particle released this time step, record its initial state
           particle%istatus = 1
           if (particle%trelease >= totimc) &
-            call this%method%save(particle, reason=0) ! reason=0: release
+            call this%method%save(particle, reason=0)
 
           ! Maximum time is the end of the time step or the particle
           ! stop time, whichever comes first, unless it's the final
diff --git a/src/Solution/ParticleTracker/CellDefn.f90 b/src/Solution/ParticleTracker/CellDefn.f90
index d609cae38c6..836f0421633 100644
--- a/src/Solution/ParticleTracker/CellDefn.f90
+++ b/src/Solution/ParticleTracker/CellDefn.f90
@@ -17,6 +17,7 @@ module CellDefnModule
     integer(I4B), public :: npolyverts !< number of vertices for cell polygon
     real(DP), public :: porosity !< cell porosity
     real(DP), public :: retfactor !< cell retardation factor
+    integer(I4B), public :: ilay !< layer number
     integer(I4B), public :: izone !< cell zone number
     integer(I4B), public :: iweaksink !< weak sink indicator
     integer(I4B), public :: inoexitface !< no exit face indicator
diff --git a/src/Solution/ParticleTracker/Method.f90 b/src/Solution/ParticleTracker/Method.f90
index 4c01f63e434..330ef5e33a1 100644
--- a/src/Solution/ParticleTracker/Method.f90
+++ b/src/Solution/ParticleTracker/Method.f90
@@ -7,6 +7,8 @@ module MethodModule
   use SubcellModule, only: SubcellType
   use ParticleModule
   use BaseDisModule, only: DisBaseType
+  use DisModule, only: DisType
+  use DisvModule, only: DisvType
   use PrtFmiModule, only: PrtFmiType
   use CellModule, only: CellType
   use CellDefnModule, only: CellDefnType
@@ -29,7 +31,7 @@ module MethodModule
   !! depending on cell geometry (implementing the strategy pattern).
   !<
   type, abstract :: MethodType
-    character(len=40), pointer, public :: type !< method name
+    character(len=40), pointer, public :: name !< method name
     logical(LGP), public :: delegates !< whether the method delegates
     type(PrtFmiType), pointer, public :: fmi => null() !< ptr to fmi
     class(CellType), pointer, public :: cell => null() !< ptr to the current cell
@@ -52,7 +54,7 @@ module MethodModule
     procedure :: save
     procedure :: track
     procedure :: try_pass
-    procedure :: prepare
+    procedure :: check
   end type MethodType
 
   abstract interface
@@ -96,7 +98,8 @@ subroutine init(this, fmi, cell, subcell, trackctl, tracktimes, &
     if (present(retfactor)) this%retfactor => retfactor
   end subroutine init
 
-  !> @brief Track particle through subdomains
+  !> @brief Track the particle over domains of the given
+  ! level until the particle terminates or tmax elapses.
   recursive subroutine track(this, particle, level, tmax)
     ! dummy
     class(MethodType), intent(inout) :: this
@@ -108,6 +111,7 @@ recursive subroutine track(this, particle, level, tmax)
     integer(I4B) :: nextlevel
     class(methodType), pointer :: submethod
 
+    ! Advance the particle over subdomains
     advancing = .true.
     nextlevel = level + 1
     do while (advancing)
@@ -132,8 +136,9 @@ subroutine try_pass(this, particle, nextlevel, advancing)
       ! otherwise pass the particle to the next subdomain.
       ! if that leaves it on a boundary, stop advancing.
       call this%pass(particle)
-      if (particle%iboundary(nextlevel - 1) .ne. 0) &
+      if (particle%iboundary(nextlevel - 1) .ne. 0) then
         advancing = .false.
+      end if
     end if
   end subroutine try_pass
 
@@ -180,26 +185,24 @@ subroutine save(this, particle, reason)
       end if
     end if
 
-    ! Save the particle's state to any registered tracking output files
-    call this%trackctl%save(particle, kper=per, &
-                            kstp=stp, reason=reason)
+    call this%trackctl%save(particle, kper=per, kstp=stp, reason=reason)
   end subroutine save
 
-  !> @brief Prepare to apply the tracking method to the particle.
+  !> @brief Check reporting/terminating conditions before tracking.
   !!
   !! Check a number of conditions determining whether to continue
-  !! tracking the particle or terminate it. This includes a check
-  !! for any reporting conditions as well.
+  !! tracking the particle or terminate it, as well as whether to
+  !! record any output data as per selected reporting conditions.
   !<
-  subroutine prepare(this, particle, cell_defn)
+  subroutine check(this, particle, cell_defn)
+    ! modules
+    use TdisModule, only: endofsimulation, totim
     ! dummy
     class(MethodType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
     type(CellDefnType), pointer, intent(inout) :: cell_defn
     ! local
-    logical(LGP) :: cell_dry
-    logical(LGP) :: locn_dry
-    integer(I4B) :: ic
+    logical(LGP) :: no_exit_face
 
     ! stop zone
     particle%izone = cell_defn%izone
@@ -210,40 +213,57 @@ subroutine prepare(this, particle, cell_defn)
       return
     end if
 
-    ! dry
-    cell_dry = this%fmi%ibdgwfsat0(cell_defn%icell) == 0
-    locn_dry = particle%z > cell_defn%top
-    if (cell_dry .or. locn_dry) then
-      if (particle%idry == 0) then
-        ! drop
-        if (cell_dry) then
-          ! if no active cell underneath position, terminate
-          ic = particle%idomain(2)
-          call this%fmi%dis%highest_active(ic, this%fmi%ibound)
-          if (this%fmi%ibound(ic) == 0) then
-            particle%advancing = .false.
-            particle%istatus = 7
-            call this%save(particle, reason=3)
-            return
-          end if
-        else if (locn_dry) then
-          particle%z = cell_defn%top
-          call this%save(particle, reason=1)
-        end if
-      else if (particle%idry == 1) then
+    ! no exit face
+    no_exit_face = cell_defn%inoexitface > 0
+
+    ! dry cell
+    if (this%fmi%ibdgwfsat0(cell_defn%icell) == 0) then
+      if (particle%idrymeth == 0) then
+        no_exit_face = .false.
+      else if (particle%idrymeth == 1) then
         ! stop
         particle%advancing = .false.
         particle%istatus = 7
         call this%save(particle, reason=3)
         return
-      else if (particle%idry == 2) then
+      else if (particle%idrymeth == 2) then
         ! stay
+        no_exit_face = .false.
         particle%advancing = .false.
+        particle%ttrack = totim
+        ! terminate if last period/step
+        if (endofsimulation) then
+          particle%istatus = 5
+          call this%save(particle, reason=3)
+          return
+        end if
+        call this%save(particle, reason=2)
       end if
+    else if (this%name /= "passtobottom" .and. particle%z > cell_defn%top) then
+      ! dry particle
+      if (particle%idrymeth == 0) then
+        ! drop to water table
+        no_exit_face = .false.
+        particle%z = cell_defn%top
+        call this%save(particle, reason=1)
+      else if (particle%idrymeth == 1) then
+        ! terminate
+        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.
+        return
+      end if
+    end if
 
-      ! cell with no exit face (mutually exclusive with
-      ! dry because a dry cell will have no exit face)
-    else if (cell_defn%inoexitface > 0) then
+    ! cell with no exit face. handle after
+    ! dry conditions because dry cells will
+    ! also have no exit face.
+    if (no_exit_face) then
       particle%advancing = .false.
       particle%istatus = 5
       call this%save(particle, reason=3)
@@ -261,6 +281,6 @@ subroutine prepare(this, particle, cell_defn)
         return
       end if
     end if
-  end subroutine prepare
+  end subroutine check
 
 end module MethodModule
diff --git a/src/Solution/ParticleTracker/MethodCellPassToBot.f90 b/src/Solution/ParticleTracker/MethodCellPassToBot.f90
index d325d6e3767..eb417457377 100644
--- a/src/Solution/ParticleTracker/MethodCellPassToBot.f90
+++ b/src/Solution/ParticleTracker/MethodCellPassToBot.f90
@@ -9,6 +9,8 @@ module MethodCellPassToBotModule
   use CellModule, only: CellType
   use SubcellModule, only: SubcellType
   use TrackControlModule, only: TrackControlType
+  use DisModule, only: DisType
+  use DisvModule, only: DisvType
   implicit none
 
   private
@@ -28,15 +30,15 @@ module MethodCellPassToBotModule
   subroutine create_method_cell_ptb(method)
     type(MethodCellPassToBotType), pointer :: method
     allocate (method)
-    allocate (method%type)
-    method%type = "passtobottom"
+    allocate (method%name)
+    method%name = "passtobottom"
     method%delegates = .false.
   end subroutine create_method_cell_ptb
 
   !> @brief Deallocate the pass-to-bottom tracking method
   subroutine deallocate (this)
     class(MethodCellPassToBotType), intent(inout) :: this
-    deallocate (this%type)
+    deallocate (this%name)
   end subroutine deallocate
 
   !> @brief Pass particle vertically and instantaneously to the cell bottom
@@ -45,11 +47,32 @@ 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, reason
 
-    call this%prepare(particle, this%cell%defn)
+    ! Check termination/reporting conditions
+    call this%check(particle, this%cell%defn)
     if (.not. particle%advancing) return
+
+    ! Pass to bottom face
     particle%z = this%cell%defn%bot
     particle%iboundary(2) = this%cell%defn%npolyverts + 2
+
+    ! Terminate if bottom layer
+    select type (dis => this%fmi%dis)
+    type is (DisType)
+      nlay = dis%nlay
+    type is (DisvType)
+      nlay = dis%nlay
+    end select
+    if (this%cell%defn%ilay == nlay) then
+      particle%istatus = 5
+      reason = 3
+    else
+      reason = 1
+    end if
+
+    ! Save datum
     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 5a58614a32d..86b1d09f4ca 100644
--- a/src/Solution/ParticleTracker/MethodCellPollock.f90
+++ b/src/Solution/ParticleTracker/MethodCellPollock.f90
@@ -36,7 +36,7 @@ subroutine create_method_cell_pollock(method)
     allocate (method)
     call create_cell_rect(cell)
     method%cell => cell
-    method%type => method%cell%type
+    method%name => method%cell%type
     method%delegates = .true.
     call create_subcell_rect(subcell)
     method%subcell => subcell
@@ -45,7 +45,7 @@ end subroutine create_method_cell_pollock
   !> @brief Destroy the tracking method
   subroutine destroy_mcp(this)
     class(MethodCellPollockType), intent(inout) :: this
-    deallocate (this%type)
+    deallocate (this%name)
   end subroutine destroy_mcp
 
   !> @brief Load subcell tracking method
@@ -129,13 +129,12 @@ subroutine apply_mcp(this, particle, tmax)
 
     select type (cell => this%cell)
     type is (CellRectType)
-      ! Prepare to apply the tracking method and
-      ! check termination / reporting conditions
-      call this%prepare(particle, cell%defn)
+      ! Check termination/reporting conditions
+      call this%check(particle, cell%defn)
       if (.not. particle%advancing) return
 
-      ! Transform particle location into local cell coordinates
-      ! (translated and rotated but not scaled relative to model)
+      ! Transform model coordinates to local cell coordinates
+      ! (translated/rotated but not scaled relative to model)
       xOrigin = cell%xOrigin
       yOrigin = cell%yOrigin
       zOrigin = cell%zOrigin
@@ -144,15 +143,13 @@ subroutine apply_mcp(this, particle, tmax)
       call particle%transform(xOrigin, yOrigin, zOrigin, &
                               sinrot, cosrot)
 
-      ! Track the particle across the cell
+      ! Track the particle over the cell
       call this%track(particle, 2, tmax)
 
-      ! Transform particle location back to model coordinates
+      ! Transform cell coordinates back to model coordinates
       call particle%transform(xOrigin, yOrigin, zOrigin, &
                               sinrot, cosrot, invert=.true.)
-
-      ! Reset transformation and eliminate accumulated roundoff error
-      call particle%transform(reset=.true.)
+      call particle%reset_transform()
     end select
   end subroutine apply_mcp
 
diff --git a/src/Solution/ParticleTracker/MethodCellPollockQuad.f90 b/src/Solution/ParticleTracker/MethodCellPollockQuad.f90
index 0d30fb06d61..c5f61f2c221 100644
--- a/src/Solution/ParticleTracker/MethodCellPollockQuad.f90
+++ b/src/Solution/ParticleTracker/MethodCellPollockQuad.f90
@@ -37,7 +37,7 @@ subroutine create_method_cell_quad(method)
     allocate (method)
     call create_cell_rect_quad(cell)
     method%cell => cell
-    method%type => method%cell%type
+    method%name => method%cell%type
     method%delegates = .true.
     call create_subcell_rect(subcell)
     method%subcell => subcell
@@ -46,7 +46,7 @@ end subroutine create_method_cell_quad
   !> @brief Deallocate the Pollock quad-refined cell method
   subroutine deallocate (this)
     class(MethodCellPollockQuadType), intent(inout) :: this
-    deallocate (this%type)
+    deallocate (this%name)
   end subroutine deallocate
 
   !> @brief Load subcell into tracking method
@@ -208,12 +208,12 @@ subroutine apply_mcpq(this, particle, tmax)
 
     select type (cell => this%cell)
     type is (CellRectQuadType)
-      ! Prepare to apply method, return early if done advancing
-      call this%prepare(particle, cell%defn)
+      ! Check termination/reporting conditions
+      call this%check(particle, cell%defn)
       if (.not. particle%advancing) return
 
-      ! Transform particle location into local cell coordinates
-      ! (translated and rotated but not scaled relative to model).
+      ! Transform model coordinates to local cell coordinates
+      ! (translated/rotated but not scaled relative to model)
       xOrigin = cell%xOrigin
       yOrigin = cell%yOrigin
       zOrigin = cell%zOrigin
@@ -222,14 +222,13 @@ subroutine apply_mcpq(this, particle, tmax)
       call particle%transform(xOrigin, yOrigin, zOrigin, &
                               sinrot, cosrot)
 
-      ! Track the particle across the cell.
+      ! Track the particle over the cell
       call this%track(particle, 2, tmax)
 
-      ! Transform particle location back to model coordinates, then
-      ! reset transformation and eliminate accumulated roundoff error.
+      ! Transform cell coordinates back to model coordinates
       call particle%transform(xOrigin, yOrigin, zOrigin, &
                               sinrot, cosrot, invert=.true.)
-      call particle%transform(reset=.true.)
+      call particle%reset_transform()
     end select
   end subroutine apply_mcpq
 
@@ -292,7 +291,7 @@ subroutine load_subcell(this, particle, subcell)
       qintl2 = cell%qintl(isc + 1)
       qextl1 = cell%qextl1(isc)
       qextl2 = cell%qextl2(isc)
-      !
+
       subcell%dx = dx
       subcell%dy = dy
       subcell%dz = dz
diff --git a/src/Solution/ParticleTracker/MethodCellTernary.f90 b/src/Solution/ParticleTracker/MethodCellTernary.f90
index be2123b7237..a04a65d231f 100644
--- a/src/Solution/ParticleTracker/MethodCellTernary.f90
+++ b/src/Solution/ParticleTracker/MethodCellTernary.f90
@@ -62,7 +62,7 @@ subroutine create_method_cell_ternary(method)
     allocate (method)
     call create_cell_poly(cell)
     method%cell => cell
-    method%type => method%cell%type
+    method%name => method%cell%type
     method%delegates = .true.
     call create_subcell_tri(subcell)
     method%subcell => subcell
@@ -71,7 +71,7 @@ end subroutine create_method_cell_ternary
   !> @brief Destroy the tracking method
   subroutine destroy_mct(this)
     class(MethodCellTernaryType), intent(inout) :: this
-    deallocate (this%type)
+    deallocate (this%name)
   end subroutine destroy_mct
 
   !> @brief Load subcell into tracking method
@@ -156,15 +156,16 @@ subroutine apply_mct(this, particle, tmax)
     ! local
     integer(I4B) :: i
 
-    ! Prepare to apply method, return early if done advancing
-    call this%prepare(particle, this%cell%defn)
-    if (.not. particle%advancing) return
-
     ! (Re)allocate type-bound arrays
     select type (cell => this%cell)
     type is (CellPolyType)
+      ! Check termination/reporting conditions
+      call this%check(particle, this%cell%defn)
+      if (.not. particle%advancing) return
+
       ! Number of vertices
       this%nverts = cell%defn%npolyverts
+
       ! (Re)allocate type-bound arrays
       if (allocated(this%xvert)) then
         deallocate (this%xvert)
@@ -178,6 +179,7 @@ subroutine apply_mct(this, particle, tmax)
         deallocate (this%xvertnext)
         deallocate (this%yvertnext)
       end if
+
       allocate (this%xvert(this%nverts))
       allocate (this%yvert(this%nverts))
       allocate (this%vne(this%nverts))
@@ -188,15 +190,18 @@ subroutine apply_mct(this, particle, tmax)
       allocate (this%iprev(this%nverts))
       allocate (this%xvertnext(this%nverts))
       allocate (this%yvertnext(this%nverts))
+
       ! Cell vertices
       do i = 1, this%nverts
         this%xvert(i) = cell%defn%polyvert(1, i)
         this%yvert(i) = cell%defn%polyvert(2, i)
       end do
+
       ! Top, bottom, and thickness
       this%ztop = cell%defn%top
       this%zbot = cell%defn%bot
       this%dz = this%ztop - this%zbot
+
       ! Shifted arrays
       do i = 1, this%nverts
         this%iprev(i) = i
diff --git a/src/Solution/ParticleTracker/MethodDis.f90 b/src/Solution/ParticleTracker/MethodDis.f90
index 913ba6ac978..fada2f62644 100644
--- a/src/Solution/ParticleTracker/MethodDis.f90
+++ b/src/Solution/ParticleTracker/MethodDis.f90
@@ -47,17 +47,17 @@ subroutine create_method_dis(method)
     type(CellRectType), pointer :: cell
 
     allocate (method)
-    allocate (method%type)
+    allocate (method%name)
     call create_cell_rect(cell)
     method%cell => cell
-    method%type = "dis"
+    method%name = "dis"
     method%delegates = .true.
   end subroutine create_method_dis
 
   !> @brief Destructor the tracking method
   subroutine deallocate (this)
     class(MethodDisType), intent(inout) :: this
-    deallocate (this%type)
+    deallocate (this%name)
   end subroutine deallocate
 
   subroutine load_cell(this, ic, cell)
@@ -193,6 +193,17 @@ subroutine load_particle(this, cell, particle)
       idiag = this%fmi%dis%con%ia(cell%defn%icell)
       ipos = idiag + inbr
       ic = dis%con%ja(ipos)
+
+      if (ic == particle%icp) then
+        print *, "bouncing!"
+        particle%istatus = 2
+        particle%advancing = .false.
+        call this%save(particle, reason=3)
+        return
+      else
+        particle%icp = particle%idomain(2)
+      end if
+
       icu = dis%get_nodeuser(ic)
       call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, &
                    irow, icol, ilay)
@@ -276,15 +287,17 @@ subroutine pass_dis(this, particle)
         particle%advancing = .false.
         call this%save(particle, reason=3)
       else
-        ! Otherwise, load cell properties into the
-        ! particle and update intercell mass flows
+        ! Update old to new cell properties
         call this%load_particle(cell, particle)
+        if (.not. particle%advancing) return
+
+        ! Update intercell mass flows
         call this%update_flowja(cell, particle)
       end if
     end select
   end subroutine pass_dis
 
-  !> @brief Apply the method to a particle
+  !> @brief Apply the structured tracking method to a particle.
   subroutine apply_dis(this, particle, tmax)
     class(MethodDisType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
@@ -335,6 +348,8 @@ subroutine load_properties(this, ic, defn)
     class(MethodDisType), intent(inout) :: this
     integer(I4B), intent(in) :: ic
     type(CellDefnType), pointer, intent(inout) :: defn
+    ! local
+    integer(I4B) :: irow, icol, ilay, icu
 
     defn%icell = ic
     defn%npolyverts = 4 ! rectangular cell always has 4 vertices
@@ -347,6 +362,12 @@ subroutine load_properties(this, ic, defn)
     defn%sat = this%fmi%gwfsat(ic)
     defn%porosity = this%porosity(ic)
     defn%retfactor = this%retfactor(ic)
+    select type (dis => this%fmi%dis)
+    type is (DisType)
+      icu = dis%get_nodeuser(ic)
+      call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, irow, icol, ilay)
+      defn%ilay = ilay
+    end select
     defn%izone = this%izone(ic)
     defn%can_be_rect = .true.
     defn%can_be_quad = .false.
diff --git a/src/Solution/ParticleTracker/MethodDisv.f90 b/src/Solution/ParticleTracker/MethodDisv.f90
index 615cdc94a79..8ba4f475af0 100644
--- a/src/Solution/ParticleTracker/MethodDisv.f90
+++ b/src/Solution/ParticleTracker/MethodDisv.f90
@@ -52,10 +52,10 @@ subroutine create_method_disv(method)
     type(CellPolyType), pointer :: cell
 
     allocate (method)
-    allocate (method%type)
+    allocate (method%name)
     call create_cell_poly(cell)
     method%cell => cell
-    method%type = "disv"
+    method%name = "disv"
     method%delegates = .true.
     call create_defn(method%neighbor)
   end subroutine create_method_disv
@@ -63,7 +63,7 @@ end subroutine create_method_disv
   !> @brief Destroy the tracking method
   subroutine deallocate (this)
     class(MethodDisvType), intent(inout) :: this
-    deallocate (this%type)
+    deallocate (this%name)
   end subroutine deallocate
 
   !> @brief Load the cell and the tracking method
@@ -155,6 +155,17 @@ subroutine load_particle(this, cell, particle)
       inbr = cell%defn%facenbr(inface)
       ipos = idiag + inbr
       ic = dis%con%ja(ipos)
+
+      if (ic == particle%icp) then
+        print *, "bouncing!"
+        particle%istatus = 2
+        particle%advancing = .false.
+        call this%save(particle, reason=3)
+        return
+      else
+        particle%icp = particle%idomain(2)
+      end if
+
       icu = dis%get_nodeuser(ic)
 
       call get_jk(icu, dis%ncpl, dis%nlay, icpl, ilay)
@@ -216,9 +227,11 @@ subroutine pass_disv(this, particle)
         particle%advancing = .false.
         call this%save(particle, reason=3)
       else
-        ! Otherwise, load cell properties into the
-        ! particle and update intercell mass flows
+        ! Update old to new cell properties
         call this%load_particle(cell, particle)
+        if (.not. particle%advancing) return
+
+        ! Update intercell mass flows
         call this%update_flowja(cell, particle)
       end if
     end select
@@ -290,11 +303,12 @@ subroutine map_neighbor(this, defn, inface, z)
     end if
   end subroutine map_neighbor
 
-  !> @brief Apply the DISV-grid method
+  !> @brief Apply the DISV tracking method to a particle.
   subroutine apply_disv(this, particle, tmax)
     class(MethodDisvType), intent(inout) :: this
     type(ParticleType), pointer, intent(inout) :: particle
     real(DP), intent(in) :: tmax
+
     call this%track(particle, 1, tmax)
   end subroutine apply_disv
 
@@ -322,6 +336,7 @@ subroutine load_properties(this, ic, defn)
     real(DP) :: top
     real(DP) :: bot
     real(DP) :: sat
+    integer(I4B) :: icu, icpl, ilay
 
     defn%icell = ic
     defn%iatop = get_iatop(this%fmi%dis%get_ncpl(), &
@@ -336,6 +351,12 @@ subroutine load_properties(this, ic, defn)
     defn%porosity = this%porosity(ic)
     defn%retfactor = this%retfactor(ic)
     defn%izone = this%izone(ic)
+    select type (dis => this%fmi%dis)
+    type is (DisvType)
+      icu = dis%get_nodeuser(ic)
+      call get_jk(icu, dis%ncpl, dis%nlay, icpl, ilay)
+      defn%ilay = ilay
+    end select
   end subroutine load_properties
 
   subroutine load_polygon(this, defn)
diff --git a/src/Solution/ParticleTracker/MethodSubcellPollock.f90 b/src/Solution/ParticleTracker/MethodSubcellPollock.f90
index de9157b8b9e..1cbd0cce45b 100644
--- a/src/Solution/ParticleTracker/MethodSubcellPollock.f90
+++ b/src/Solution/ParticleTracker/MethodSubcellPollock.f90
@@ -36,14 +36,14 @@ subroutine create_method_subcell_pollock(method)
     allocate (method)
     call create_subcell_rect(subcell)
     method%subcell => subcell
-    method%type => method%subcell%type
+    method%name => method%subcell%type
     method%delegates = .false.
   end subroutine create_method_subcell_pollock
 
   !> @brief Deallocate the Pollock's subcell method
   subroutine deallocate (this)
     class(MethodSubcellPollockType), intent(inout) :: this
-    deallocate (this%type)
+    deallocate (this%name)
   end subroutine deallocate
 
   !> @brief Apply Pollock's method to a rectangular subcell
@@ -73,6 +73,7 @@ subroutine apply_msp(this, particle, tmax)
       call particle%transform(xOrigin, yOrigin)
       call this%track_subcell(subcell, particle, tmax)
       call particle%transform(xOrigin, yOrigin, invert=.true.)
+      call particle%reset_transform()
     end select
   end subroutine apply_msp
 
diff --git a/src/Solution/ParticleTracker/MethodSubcellTernary.f90 b/src/Solution/ParticleTracker/MethodSubcellTernary.f90
index 057eca00050..97637b4c59b 100644
--- a/src/Solution/ParticleTracker/MethodSubcellTernary.f90
+++ b/src/Solution/ParticleTracker/MethodSubcellTernary.f90
@@ -38,14 +38,14 @@ subroutine create_method_subcell_ternary(method)
     allocate (method)
     call create_subcell_tri(subcell)
     method%subcell => subcell
-    method%type => method%subcell%type
+    method%name => method%subcell%type
     method%delegates = .false.
   end subroutine create_method_subcell_ternary
 
   !> @brief Deallocate the ternary subcell tracking method.
   subroutine deallocate (this)
     class(MethodSubcellTernaryType), intent(inout) :: this
-    deallocate (this%type)
+    deallocate (this%name)
   end subroutine deallocate
 
   !> @brief Apply the ternary subcell tracking method.
diff --git a/src/Solution/ParticleTracker/Particle.f90 b/src/Solution/ParticleTracker/Particle.f90
index 188df9e4c08..32751d71247 100644
--- a/src/Solution/ParticleTracker/Particle.f90
+++ b/src/Solution/ParticleTracker/Particle.f90
@@ -40,11 +40,12 @@ module ParticleModule
     ! stop criteria
     integer(I4B), public :: istopweaksink !< weak sink option (0: do not stop, 1: stop)
     integer(I4B), public :: istopzone !< stop zone number
-    integer(I4B), public :: idry !< stop in dry cells
+    integer(I4B), public :: idrymeth !< dry tracking method
     ! state
-    integer(I4B), allocatable, public :: idomain(:) !< tracking domain hierarchy
+    integer(I4B), allocatable, public :: idomain(:) !< tracking domain hierarchy ! TODO: rename to itdomain?
     integer(I4B), allocatable, public :: iboundary(:) !< tracking domain boundaries
-    integer(I4B), public :: icu !< user cell (node) number
+    integer(I4B), public :: icp !< previous cell number (reduced)
+    integer(I4B), public :: icu !< user cell number
     integer(I4B), public :: ilay !< grid layer
     integer(I4B), public :: izone !< zone number
     integer(I4B), public :: istatus !< tracking status
@@ -69,6 +70,7 @@ module ParticleModule
     procedure, public :: get_model_coords
     procedure, public :: load_particle
     procedure, public :: transform => transform_coords
+    procedure, public :: reset_transform
   end type ParticleType
 
   !> @brief Structure of arrays to store particles.
@@ -82,11 +84,12 @@ module ParticleModule
     ! stopping criteria
     integer(I4B), dimension(:), pointer, public, contiguous :: istopweaksink !< weak sink option: 0 = do not stop, 1 = stop
     integer(I4B), dimension(:), pointer, public, contiguous :: istopzone !< stop zone number
-    integer(I4B), dimension(:), pointer, public, contiguous :: idry !< stop in dry cells
+    integer(I4B), dimension(:), pointer, public, contiguous :: idrymeth !< stop in dry cells
     ! state
     integer(I4B), dimension(:, :), pointer, public, contiguous :: idomain !< array of indices for domains in the tracking domain hierarchy
     integer(I4B), dimension(:, :), pointer, public, contiguous :: iboundary !< array of indices for tracking domain boundaries
-    integer(I4B), dimension(:), pointer, public, contiguous :: icu !< cell number (user, not reduced)
+    integer(I4B), dimension(:), pointer, public, contiguous :: icp !< previous cell number (reduced)
+    integer(I4B), dimension(:), pointer, public, contiguous :: icu !< cell number (user)
     integer(I4B), dimension(:), pointer, public, contiguous :: ilay !< layer
     integer(I4B), dimension(:), pointer, public, contiguous :: izone !< current zone number
     integer(I4B), dimension(:), pointer, public, contiguous :: istatus !< particle status
@@ -128,6 +131,7 @@ subroutine allocate_particle_store(this, np, mempath)
     call mem_allocate(this%irpt, np, 'PLIRPT', mempath)
     call mem_allocate(this%iprp, np, 'PLIPRP', mempath)
     call mem_allocate(this%name, LENBOUNDNAME, np, 'PLNAME', mempath)
+    call mem_allocate(this%icp, np, 'PLICP', mempath)
     call mem_allocate(this%icu, np, 'PLICU', mempath)
     call mem_allocate(this%ilay, np, 'PLILAY', mempath)
     call mem_allocate(this%izone, np, 'PLIZONE', mempath)
@@ -140,7 +144,7 @@ subroutine allocate_particle_store(this, np, mempath)
     call mem_allocate(this%ttrack, np, 'PLTTRACK', mempath)
     call mem_allocate(this%istopweaksink, np, 'PLISTOPWEAKSINK', mempath)
     call mem_allocate(this%istopzone, np, 'PLISTOPZONE', mempath)
-    call mem_allocate(this%idry, np, 'PLIDRY', mempath)
+    call mem_allocate(this%idrymeth, np, 'PLIDRYMETH', mempath)
     call mem_allocate(this%ifrctrn, np, 'PLIFRCTRN', mempath)
     call mem_allocate(this%iexmeth, np, 'PLIEXMETH', mempath)
     call mem_allocate(this%extol, np, 'PLEXTOL', mempath)
@@ -158,6 +162,7 @@ subroutine deallocate (this, mempath)
     call mem_deallocate(this%iprp, 'PLIPRP', mempath)
     call mem_deallocate(this%irpt, 'PLIRPT', mempath)
     call mem_deallocate(this%name, 'PLNAME', mempath)
+    call mem_deallocate(this%icp, 'PLICP', mempath)
     call mem_deallocate(this%icu, 'PLICU', mempath)
     call mem_deallocate(this%ilay, 'PLILAY', mempath)
     call mem_deallocate(this%izone, 'PLIZONE', mempath)
@@ -170,7 +175,7 @@ subroutine deallocate (this, mempath)
     call mem_deallocate(this%ttrack, 'PLTTRACK', mempath)
     call mem_deallocate(this%istopweaksink, 'PLISTOPWEAKSINK', mempath)
     call mem_deallocate(this%istopzone, 'PLISTOPZONE', mempath)
-    call mem_deallocate(this%idry, 'PLIDRY', mempath)
+    call mem_deallocate(this%idrymeth, 'PLIDRYMETH', mempath)
     call mem_deallocate(this%ifrctrn, 'PLIFRCTRN', mempath)
     call mem_deallocate(this%iexmeth, 'PLIEXMETH', mempath)
     call mem_deallocate(this%extol, 'PLEXTOL', mempath)
@@ -181,7 +186,7 @@ end subroutine deallocate
 
   !> @brief Reallocate particle arrays
   subroutine resize(this, np, mempath)
-    ! -- dummy
+    ! dummy
     class(ParticleStoreType), intent(inout) :: this !< particle store
     integer(I4B), intent(in) :: np !< number of particles
     character(*), intent(in) :: mempath !< path to memory
@@ -191,6 +196,7 @@ subroutine resize(this, np, mempath)
     call mem_reallocate(this%iprp, np, 'PLIPRP', mempath)
     call mem_reallocate(this%irpt, np, 'PLIRPT', mempath)
     call mem_reallocate(this%name, LENBOUNDNAME, np, 'PLNAME', mempath)
+    call mem_reallocate(this%icp, np, 'PLICP', mempath)
     call mem_reallocate(this%icu, np, 'PLICU', mempath)
     call mem_reallocate(this%ilay, np, 'PLILAY', mempath)
     call mem_reallocate(this%izone, np, 'PLIZONE', mempath)
@@ -203,7 +209,7 @@ subroutine resize(this, np, mempath)
     call mem_reallocate(this%ttrack, np, 'PLTTRACK', mempath)
     call mem_reallocate(this%istopweaksink, np, 'PLISTOPWEAKSINK', mempath)
     call mem_reallocate(this%istopzone, np, 'PLISTOPZONE', mempath)
-    call mem_reallocate(this%idry, np, 'PLIDRY', mempath)
+    call mem_reallocate(this%idrymeth, np, 'PLIDRYMETH', mempath)
     call mem_reallocate(this%ifrctrn, np, 'PLIFRCTRN', mempath)
     call mem_reallocate(this%iexmeth, np, 'PLIEXMETH', mempath)
     call mem_reallocate(this%extol, np, 'PLEXTOL', mempath)
@@ -224,7 +230,7 @@ subroutine load_particle(this, store, imdl, iprp, ip)
     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%transform(reset=.true.)
+    call this%reset_transform()
     this%imdl = imdl
     this%iprp = iprp
     this%irpt = store%irpt(ip)
@@ -232,7 +238,8 @@ subroutine load_particle(this, store, imdl, iprp, ip)
     this%name = store%name(ip)
     this%istopweaksink = store%istopweaksink(ip)
     this%istopzone = store%istopzone(ip)
-    this%idry = store%idry(ip)
+    this%idrymeth = store%idrymeth(ip)
+    this%icp = store%icp(ip)
     this%icu = store%icu(ip)
     this%ilay = store%ilay(ip)
     this%izone = store%izone(ip)
@@ -267,7 +274,8 @@ subroutine save_particle(this, particle, ip)
     this%name(ip) = particle%name
     this%istopweaksink(ip) = particle%istopweaksink
     this%istopzone(ip) = particle%istopzone
-    this%idry(ip) = particle%idry
+    this%idrymeth(ip) = particle%idrymeth
+    this%icp(ip) = particle%icp
     this%icu(ip) = particle%icu
     this%ilay(ip) = particle%ilay
     this%izone(ip) = particle%izone
@@ -292,9 +300,17 @@ subroutine save_particle(this, particle, ip)
     this%extend(ip) = particle%iextend
   end subroutine save_particle
 
-  !> @brief Apply the given global-to-local transformation to the particle.
+  !> @brief Transform particle coordinates.
+  !!
+  !! Apply a translation and/or rotation to particle coordinates.
+  !! No rescaling. It's also possible to invert a transformation.
+  !!
+  !! Upon inverting, be sure to untransform to eliminate roundoff
+  !! error accumulated in the particle origin/rotation variables,
+  !! and unsetting the transformed flag.
+  !<
   subroutine transform_coords(this, xorigin, yorigin, zorigin, &
-                              sinrot, cosrot, invert, reset)
+                              sinrot, cosrot, invert)
     use GeomUtilModule, only: transform, compose
     class(ParticleType), intent(inout) :: this !< particle
     real(DP), intent(in), optional :: xorigin !< x coordinate of origin
@@ -303,44 +319,33 @@ subroutine transform_coords(this, xorigin, yorigin, zorigin, &
     real(DP), intent(in), optional :: sinrot !< sine of rotation angle
     real(DP), intent(in), optional :: cosrot !< cosine of rotation angle
     logical(LGP), intent(in), optional :: invert !< whether to invert
-    logical(LGP), intent(in), optional :: reset !< whether to reset
-
-    ! Reset if requested
-    if (present(reset)) then
-      if (reset) then
-        this%xorigin = DZERO
-        this%yorigin = DZERO
-        this%zorigin = DZERO
-        this%sinrot = DZERO
-        this%cosrot = DONE
-        this%cosrot = DONE
-        this%transformed = .false.
-        return
-      end if
-    end if
 
-    ! Otherwise, transform coordinates
     call transform(this%x, this%y, this%z, &
                    this%x, this%y, this%z, &
                    xorigin, yorigin, zorigin, &
                    sinrot, cosrot, invert)
 
-    ! Modify transformation from model coordinates to particle's new
-    ! local coordinates by incorporating this latest transformation
     call compose(this%xorigin, this%yorigin, this%zorigin, &
                  this%sinrot, this%cosrot, &
                  xorigin, yorigin, zorigin, &
                  sinrot, cosrot, invert)
 
-    ! Set isTransformed flag to true. Note that there is no check
-    ! to see whether the modification brings the coordinates back
-    ! to model coordinates (in which case the origin would be very
-    ! close to zero and sinrot and cosrot would be very close to 0.
-    ! and 1., respectively, allowing for roundoff error).
     this%transformed = .true.
   end subroutine transform_coords
 
-  !> @brief Return the particle's model (global) coordinates.
+  subroutine reset_transform(this)
+    class(ParticleType), intent(inout) :: this !< particle
+
+    this%xorigin = DZERO
+    this%yorigin = DZERO
+    this%zorigin = DZERO
+    this%sinrot = DZERO
+    this%cosrot = DONE
+    this%cosrot = DONE
+    this%transformed = .false.
+  end subroutine reset_transform
+
+  !> @brief Return the particle's model coordinates.
   subroutine get_model_coords(this, x, y, z)
     use GeomUtilModule, only: transform, compose
     class(ParticleType), intent(inout) :: this !< particle
@@ -349,12 +354,11 @@ subroutine get_model_coords(this, x, y, z)
     real(DP), intent(out) :: z !< z coordinate
 
     if (this%transformed) then
-      ! Transform back from local to model coordinates
+      ! Untransform coordinates
       call transform(this%x, this%y, this%z, x, y, z, &
                      this%xorigin, this%yorigin, this%zorigin, &
-                     this%sinrot, this%cosrot, .true.)
+                     this%sinrot, this%cosrot, invert=.true.)
     else
-      ! Already in model coordinates
       x = this%x
       y = this%y
       z = this%z
diff --git a/src/Utilities/Sim.f90 b/src/Utilities/Sim.f90
index a594bc79c1f..75d292eda6e 100644
--- a/src/Utilities/Sim.f90
+++ b/src/Utilities/Sim.f90
@@ -310,7 +310,7 @@ end subroutine store_note
   !<
   subroutine ustop(stopmess, ioutlocal)
     ! -- dummy variables
-    character, optional, intent(in) :: stopmess * (*) !< optional message to print before
+    character(len=*), optional, intent(in) :: stopmess !< optional message to print before
                                                       !! stopping the simulation
     integer(I4B), optional, intent(in) :: ioutlocal !< optional output file to
                                                     !! final message to
@@ -331,7 +331,7 @@ end subroutine ustop
   !<
   subroutine print_final_message(stopmess, ioutlocal)
     ! -- dummy variables
-    character, optional, intent(in) :: stopmess * (*) !< optional message to print before
+    character(len=*), optional, intent(in) :: stopmess !< optional message to print before
                                                       !! stopping the simulation
     integer(I4B), optional, intent(in) :: ioutlocal !< optional output file to
                                                     !! final message to