report prev cell if we detect backtracking and terminate

src/Solution/ParticleTracker/MethodCellPassToBot.f90

@@ -9,8 +9,6 @@ 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
@@ -47,8 +45,6 @@ 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
 
     ! Check termination/reporting conditions
     call this%check(particle, this%cell%defn)
@@ -58,20 +54,6 @@ subroutine apply_ptb(this, particle, tmax)
     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

src/Solution/ParticleTracker/MethodDis.f90

@@ -195,12 +195,19 @@ subroutine load_particle(this, cell, particle)
       ic = dis%con%ja(ipos)
 
       if (ic == particle%icp) then
-        particle%istatus = 2
+        ! terminate in the previous cell.
+        ! got here by falling through the
+        ! bottom of a well.. TODO: check
+        ! if entry/exit face is top/bot?
         particle%advancing = .false.
+        particle%idomain(2) = particle%icp
+        particle%istatus = 2
+        particle%izone = particle%izp
         call this%save(particle, reason=3)
         return
       else
         particle%icp = particle%idomain(2)
+        particle%izp = particle%izone
       end if
 
       icu = dis%get_nodeuser(ic)

src/Solution/ParticleTracker/MethodDisv.f90

@@ -157,12 +157,19 @@ subroutine load_particle(this, cell, particle)
       ic = dis%con%ja(ipos)
 
       if (ic == particle%icp) then
-        particle%istatus = 2
+        ! terminate in the previous cell.
+        ! got here by falling through the
+        ! bottom of a well.. TODO: check
+        ! if entry/exit face is top/bot?
         particle%advancing = .false.
+        particle%idomain(2) = particle%icp
+        particle%istatus = 2
+        particle%izone = particle%izp
         call this%save(particle, reason=3)
         return
       else
         particle%icp = particle%idomain(2)
+        particle%izp = particle%izone
       end if
 
       icu = dis%get_nodeuser(ic)

src/Solution/ParticleTracker/Particle.f90

@@ -47,7 +47,8 @@ module ParticleModule
     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 :: izone !< current zone number
+    integer(I4B), public :: izp !< previous zone number
     integer(I4B), public :: istatus !< tracking status
     real(DP), public :: x !< x coordinate
     real(DP), public :: y !< y coordinate
@@ -92,6 +93,7 @@ module ParticleModule
     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 :: izp !< previous zone number
     integer(I4B), dimension(:), pointer, public, contiguous :: istatus !< particle status
     real(DP), dimension(:), pointer, public, contiguous :: x !< model x coord of particle
     real(DP), dimension(:), pointer, public, contiguous :: y !< model y coord of particle
@@ -135,6 +137,7 @@ subroutine allocate_particle_store(this, np, 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)
+    call mem_allocate(this%izp, np, 'PLIZP', mempath)
     call mem_allocate(this%istatus, np, 'PLISTATUS', mempath)
     call mem_allocate(this%x, np, 'PLX', mempath)
     call mem_allocate(this%y, np, 'PLY', mempath)
@@ -166,6 +169,7 @@ subroutine deallocate (this, mempath)
     call mem_deallocate(this%icu, 'PLICU', mempath)
     call mem_deallocate(this%ilay, 'PLILAY', mempath)
     call mem_deallocate(this%izone, 'PLIZONE', mempath)
+    call mem_deallocate(this%izp, 'PLIZP', mempath)
     call mem_deallocate(this%istatus, 'PLISTATUS', mempath)
     call mem_deallocate(this%x, 'PLX', mempath)
     call mem_deallocate(this%y, 'PLY', mempath)
@@ -200,6 +204,7 @@ subroutine resize(this, np, 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)
+    call mem_reallocate(this%izp, np, 'PLIZP', mempath)
     call mem_reallocate(this%istatus, np, 'PLISTATUS', mempath)
     call mem_reallocate(this%x, np, 'PLX', mempath)
     call mem_reallocate(this%y, np, 'PLY', mempath)
@@ -243,6 +248,7 @@ subroutine load_particle(this, store, imdl, iprp, ip)
     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)
@@ -279,6 +285,7 @@ subroutine save_particle(this, particle, ip)
     this%icu(ip) = particle%icu
     this%ilay(ip) = particle%ilay
     this%izone(ip) = particle%izone
+    this%izp(ip) = particle%izp
     this%istatus(ip) = particle%istatus
     this%x(ip) = particle%x
     this%y(ip) = particle%y