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/test_prt_dry.py b/autotest/test_prt_dry.py
index 72fb1b7ee0d..edadf697fcc 100644
--- a/autotest/test_prt_dry.py
+++ b/autotest/test_prt_dry.py
@@ -72,6 +72,15 @@
 
 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(
@@ -269,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]
@@ -372,13 +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)
-    if "stay" in name:
-        assert len(actual[actual.t == user_time]) == 5
+        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/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/MethodCellPassToBot.f90 b/src/Solution/ParticleTracker/MethodCellPassToBot.f90
index c56dab1b104..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,6 +47,8 @@ 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, tmax)
@@ -54,7 +58,20 @@ subroutine apply_ptb(this, particle, tmax)
     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/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.
   !!