fix 2014

src/Model/ParticleTracking/prt-fmi.f90

@@ -35,74 +35,74 @@ module PrtFmiModule
 
   !> @brief Create a new PrtFmi object
   subroutine fmi_cr(fmiobj, name_model, inunit, iout)
-    ! -- dummy
+    ! dummy
     type(PrtFmiType), pointer :: fmiobj
     character(len=*), intent(in) :: name_model
     integer(I4B), intent(inout) :: inunit
     integer(I4B), intent(in) :: iout
     !
-    ! -- Create the object
+    ! Create the object
     allocate (fmiobj)
     !
-    ! -- create name and memory path
+    ! create name and memory path
     call fmiobj%set_names(1, name_model, 'FMI', 'FMI')
     fmiobj%text = text
     !
-    ! -- Allocate scalars
+    ! Allocate scalars
     call fmiobj%allocate_scalars()
     !
-    ! -- Set variables
+    ! Set variables
     fmiobj%inunit = inunit
     fmiobj%iout = iout
     !
-    ! -- Initialize block parser
+    ! Initialize block parser
     call fmiobj%parser%Initialize(fmiobj%inunit, fmiobj%iout)
     !
-    ! -- Assign dependent variable label
+    ! Assign dependent variable label
     fmiobj%depvartype = 'TRACKS          '
 
   end subroutine fmi_cr
 
   !> @brief Time step advance
   subroutine fmi_ad(this)
-    ! -- modules
+    ! modules
     use ConstantsModule, only: DHDRY
-    ! -- dummy
+    ! dummy
     class(PrtFmiType) :: this
-    ! -- local
+    ! local
     integer(I4B) :: n
     character(len=15) :: nodestr
     character(len=*), parameter :: fmtdry = &
      &"(/1X,'WARNING: DRY CELL ENCOUNTERED AT ',a,';  RESET AS INACTIVE')"
     character(len=*), parameter :: fmtrewet = &
      &"(/1X,'DRY CELL REACTIVATED AT ', a)"
-    
+
     ! Set flag to indicated that flows are being updated
     this%iflowsupdated = 1
-    
+
     ! If reading flows from a budget file, read the next set of records
     if (this%iubud /= 0) &
       call this%advance_bfr()
-    
+
     ! If reading heads from a head file, read the next set of records
     if (this%iuhds /= 0) &
       call this%advance_hfr()
-    
+
     ! If mover flows are being read from file, read the next set of records
     if (this%iumvr /= 0) &
       call this%mvrbudobj%bfr_advance(this%dis, this%iout)
-    
+
     ! Accumulate flows
     call this%accumulate_flows()
-    
+
     ! Deactivate dry cells and reactivate rewet cells
     do n = 1, this%dis%nodes
       if (this%gwfsat(n) > DZERO) then
         this%ibdgwfsat0(n) = 1
       else
         this%ibdgwfsat0(n) = 0
       end if
-      
+
       if (this%ibound(n) > 0) then
         if (this%gwfhead(n) == DHDRY) then
           this%ibound(n) = 0
@@ -122,17 +122,17 @@ end subroutine fmi_ad
 
   !> @brief Define the flow model interface
   subroutine prtfmi_df(this, dis, idryinactive)
-    ! -- modules
+    ! modules
     use SimModule, only: store_error
-    ! -- dummy
+    ! dummy
     class(PrtFmiType) :: this
     class(DisBaseType), pointer, intent(in) :: dis
     integer(I4B), intent(in) :: idryinactive
     !
-    ! -- Call parent class define
+    ! Call parent class define
     call this%FlowModelInterfaceType%fmi_df(dis, idryinactive)
     !
-    ! -- Allocate arrays
+    ! Allocate arrays
     allocate (this%StorageFlows(this%dis%nodes))
     allocate (this%SourceFlows(this%dis%nodes))
     allocate (this%SinkFlows(this%dis%nodes))
@@ -143,9 +143,9 @@ end subroutine prtfmi_df
   !> @brief Accumulate flows
   subroutine accumulate_flows(this)
     implicit none
-    ! -- dummy
+    ! dummy
     class(PrtFmiType) :: this
-    ! -- local
+    ! local
     integer(I4B) :: j, i, ip, ib
     integer(I4B) :: ioffset, iflowface, iauxiflowface
     real(DP) :: qbnd

src/Model/ParticleTracking/prt-prp.f90

@@ -463,15 +463,11 @@ subroutine release(this, ip, trelease)
     particle%ilay = ilay
     particle%izone = this%rptzone(ic)
     particle%istatus = 0
-    if (this%ibound(ic) == 0) then
-      ! Handle inactive cells:
-      ! If drape option activated, release in highest active
-      ! cell vertically below release point. Otherwise don't
-      ! release the particle, just terminate it unreleased.
+    if (this%ibound(ic) == 0 .or. this%fmi%ibdgwfsat0(ic) == 0) then
       if (this%idrape > 0) then
         call this%dis%highest_active(ic, this%ibound)
-      else
-        particle%istatus = 8 ! permanently unreleased
+      else if (this%ibound(ic) == 0) then
+        particle%istatus = 8
       end if
     end if
     particle%x = x

src/Solution/ParticleTracker/Method.f90

@@ -2,6 +2,7 @@
 module MethodModule
 
   use KindModule, only: DP, I4B, LGP
+  use ConstantsModule, only: DZERO
   use ErrorUtilModule, only: pstop
   use SubcellModule, only: SubcellType
   use ParticleModule
@@ -11,6 +12,7 @@ module MethodModule
   use CellDefnModule, only: CellDefnType
   use TrackControlModule, only: TrackControlType
   use TimeSelectModule, only: TimeSelectType
+  use MathUtilModule, only: is_close
   implicit none
 
   private
@@ -198,28 +200,34 @@ subroutine update(this, particle, cell_defn)
     type(CellDefnType), pointer, intent(inout) :: cell_defn
 
     particle%izone = cell_defn%izone
+    if (is_close(cell_defn%top - cell_defn%bot, DZERO)) then
+      particle%advancing = .false.
+      particle%istatus = 7
+      call this%save(particle, reason=3)
+      return
+    end if
     if (cell_defn%izone .ne. 0) then
       if (particle%istopzone .eq. cell_defn%izone) then
         particle%advancing = .false.
         particle%istatus = 6
-        call this%save(particle, reason=3) ! particle terminated
+        call this%save(particle, reason=3)
         return
       end if
     end if
     if (cell_defn%inoexitface .ne. 0) then
       particle%advancing = .false.
       particle%istatus = 5
-      call this%save(particle, reason=3) ! particle terminated
+      call this%save(particle, reason=3)
       return
     end if
     if (cell_defn%iweaksink .ne. 0) then
-      if (particle%istopweaksink .ne. 0) then
-        particle%advancing = .false.
-        particle%istatus = 3
-        call this%save(particle, reason=3) ! particle terminated
+      if (particle%istopweaksink == 0) then
+        call this%save(particle, reason=4)
         return
       else
-        call this%save(particle, reason=4) ! particle exited weak sink
+        particle%advancing = .false.
+        particle%istatus = 3
+        call this%save(particle, reason=3)
         return
       end if
     end if