finally

src/Solution/ParticleTracker/MethodCellTernary.f90

@@ -267,9 +267,6 @@ subroutine load_subcell(this, particle, subcell)
           y1 = this%yvertnext(iv0)
           x2 = this%xctr
           y2 = this%yctr
-
-          ! print *, "subcell in load_subcell: ", ic, isc, x0, y0, x1, y1, x2, y2
-
           x1rel = x1 - x0
           y1rel = y1 - y0
           x2rel = x2 - x0
@@ -362,7 +359,7 @@ subroutine vertvelo(this)
     real(DP), allocatable, dimension(:) :: vm0y
     real(DP), allocatable, dimension(:) :: vm1x
     real(DP), allocatable, dimension(:) :: vm1y
-    real(DP), allocatable, dimension(:,:) :: poly
+    real(DP), allocatable, dimension(:, :) :: poly
     integer(I4B) :: nvert
 
     select type (cell => this%cell)
@@ -413,18 +410,16 @@ subroutine vertvelo(this)
         this%xctr = sum((this%xvert + this%xvertnext) * wk1) / sixa
         this%yctr = sum((this%yvert + this%yvertnext) * wk1) / sixa
       end if
+
+      ! TODO: do we need to check if center is in polygon?
+      ! only if using the centroid method which misbehaves
+      !
       ! allocate(poly(2, nvert))
       ! poly(1,:) = this%xvert
       ! poly(2,:) = this%yvert
       ! if (.not. point_in_polygon(this%xctr, this%yctr, poly)) then
-      !   print *, "not in poly!!!!!! ", this%cell%defn%icell
-      !   if (nvert == 3) then
-      !     this%xctr = sum(this%xvert) / 3.0_DP
-      !     this%yctr = sum(this%yvert) / 3.0_DP
-      !   else
-      !     print *, "error -- centroid not in cell ", this%cell%defn%icell, this%xctr, this%yctr
-      !     call pstop(1)
-      !   end if
+      !   print *, "error -- centroid not in cell ", this%cell%defn%icell, this%xctr, this%yctr
+      !   call pstop(1)
       ! end if
       ! deallocate(poly)
 

src/Solution/ParticleTracker/MethodSubcellTernary.f90

@@ -31,9 +31,9 @@ module MethodSubcellTernaryModule
 
   !> @brief Create a new ternary subcell tracking method.
   subroutine create_method_subcell_ternary(method)
-    ! -- dummy
+    ! dummy
     type(MethodSubcellTernaryType), pointer :: method
-    ! -- local
+    ! local
     type(SubcellTriType), pointer :: subcell
 
     allocate (method)
@@ -167,8 +167,6 @@ subroutine track_subcell(this, subcell, particle, tmax)
     vzbot = subcell%vzbot
     vztop = subcell%vztop
 
-    ! print *, "subcell in track_subcell: ", subcell%icell, subcell%isubcell, x0, y0, x1, y1, x2, y2
-
     ! Transform coordinates to the "canonical" configuration:
     ! barycentric in two dimensions with alpha, beta & gamma
     ! such that at f2 alpha = 0, f0 beta = 0, f1 gamma = 0.
@@ -223,14 +221,19 @@ subroutine track_subcell(this, subcell, particle, tmax)
       ! Exit through lateral subcell face
       exitFace = itrifaceexit
       dtexit = dtexitxy
-    else if (dtexitz < dtexitxy) then
+    else if (dtexitz < dtexitxy .and. dtexitz >= 0.0_DP) then
       ! Exit through top or bottom
       if (itopbotexit == -1) then
         exitFace = 4
       else
         exitFace = 5
       end if
       dtexit = dtexitz
+    else
+      particle%istatus = 9
+      particle%advancing = .false.
+      call this%save(particle, reason=3)
+      return
     end if
     texit = particle%ttrack + dtexit
     t0 = particle%ttrack
@@ -262,17 +265,17 @@ subroutine track_subcell(this, subcell, particle, tmax)
     ! (local, unscaled) and other properties. The particle may at this
     ! point lie on a boundary of the subcell or may still be within it.
     if (texit .gt. tmax) then
-      ! -- The computed exit time is greater than the maximum time, so set
-      ! -- final time for particle trajectory equal to maximum time.
+      ! The computed exit time is greater than the maximum time, so set
+      ! final time for particle trajectory equal to maximum time.
       t = tmax
       dt = t - t0
       exitFace = 0
       particle%istatus = 1
       particle%advancing = .false.
       reason = 2 ! timestep end
     else
-      ! -- The computed exit time is less than or equal to the maximum time,
-      ! -- so set final time for particle trajectory equal to exit time.
+      ! The computed exit time is less than or equal to the maximum time,
+      ! so set final time for particle trajectory equal to exit time.
       t = texit
       dt = dtexit
       reason = 1 ! (sub)cell transition
@@ -426,13 +429,6 @@ subroutine calculate_dt(v1, v2, dx, xL, v, dvdx, &
       end if
     end if
 
-    if (vr < DZERO) then
-      print *, "vr v1a v2a: ", vr, v1a, v2a
-      print *, "vr1 vr2: ", vr1, vr2
-      print *, "v1 v2: ", v1, v2
-      print *, "v: ", v
-    end if
-
     ! Compute travel time to exit face. Return with status = 0
     dt = log(abs(vr)) / dvdx
     status = 0
@@ -498,13 +494,13 @@ subroutine calculate_xyz_position(dt, rxx, rxy, ryx, ryy, sxx, sxy, syy, &
     else
       ! otherwise calculate z.
       if (izstatus .eq. 2) then
-        ! -- vz uniformly zero
+        ! vz uniformly zero
         z = zi
       else if (izstatus .eq. 1) then
-        ! -- vz uniform, nonzero
+        ! vz uniform, nonzero
         z = zi + vzi * dt
       else
-        ! -- vz nonuniform
+        ! vz nonuniform
         z = zbot + (vzi * dexp(az * dt) - vzbot) / az
       end if
     end if

src/Solution/ParticleTracker/TernarySolveTrack.f90

@@ -304,14 +304,6 @@ subroutine step_analytical(t, alp, bet)
     real(DP) :: bet
 
     if (icase .eq. 1) then
-      if (kper > 86) then
-        print *, "----"
-        print *, ca1
-        print *, ca2
-        print *, ca3
-        print *, waa
-        print *, t
-      end if
       alp = ca1 + ca2 * dexp(waa * t) + ca3 * dexp(wbb * t)
       bet = cb1 + cb2 * dexp(wbb * t)
     else if (icase .eq. -1) then