refactor(prt-fmi): add boundary face awareness (#2477)

Store IFLOWFACE assigned boundary faces in PRT FMI. Revision of #2446 following review. No functional changes.
Factored out of #2475, as this will be useful to fix the top face exit issue described in #2468 too.

Author: wpbonelli

src/Model/ParticleTracking/prt-fmi.f90

@@ -1,7 +1,8 @@
 module PrtFmiModule
 
-  use KindModule, only: DP, I4B
-  use ConstantsModule, only: DZERO, LENAUXNAME, LENPACKAGENAME
+  use KindModule, only: DP, I4B, LGP
+  use ErrorUtilModule, only: pstop
+  use ConstantsModule, only: DZERO, LENAUXNAME, LENPACKAGENAME, LENVARNAME
   use SimModule, only: store_error
   use SimVariablesModule, only: errmsg
   use FlowModelInterfaceModule, only: FlowModelInterfaceType
@@ -22,12 +23,15 @@ module PrtFmiModule
     double precision, allocatable, public :: SinkFlows(:) ! cell sink flows array
     double precision, allocatable, public :: StorageFlows(:) ! cell storage flows array
     double precision, allocatable, public :: BoundaryFlows(:) ! cell boundary flows array
+    integer(I4B), allocatable, public :: BoundaryFaces(:) ! bitmask of assigned boundary faces
 
   contains
 
     procedure :: fmi_ad
     procedure :: fmi_df => prtfmi_df
     procedure, private :: accumulate_flows
+    procedure :: mark_boundary_face
+    procedure :: is_boundary_face
 
   end type PrtFmiType
 
@@ -76,8 +80,8 @@ subroutine fmi_ad(this)
      &"(/1X,'DRY CELL REACTIVATED AT ', a)"
     !
     ! Set flag to indicated that flows are being updated.  For the case where
-    !    flows may be reused (only when flows are read from a file) then set
-    !    the flag to zero to indicated that flows were not updated
+    ! flows may be reused (only when flows are read from a file) then set
+    ! the flag to zero to indicated that flows were not updated
     this%iflowsupdated = 1
     !
     ! If reading flows from a budget file, read the next set of records
@@ -102,7 +106,7 @@ subroutine fmi_ad(this)
     do n = 1, this%dis%nodes
       !
       ! Calculate the ibound-like array that has 0 if saturation
-      !    is zero and 1 otherwise
+      ! is zero and 1 otherwise
       if (this%gwfsat(n) > DZERO) then
         this%ibdgwfsat0(n) = 1
       else
@@ -150,6 +154,7 @@ subroutine prtfmi_df(this, dis, idryinactive)
     allocate (this%SourceFlows(this%dis%nodes))
     allocate (this%SinkFlows(this%dis%nodes))
     allocate (this%BoundaryFlows(this%dis%nodes * this%max_faces))
+    allocate (this%BoundaryFaces(this%dis%nodes))
 
   end subroutine prtfmi_df
 
@@ -160,7 +165,7 @@ subroutine accumulate_flows(this)
     class(PrtFmiType) :: this
     ! local
     integer(I4B) :: j, i, ip, ib
-    integer(I4B) :: ioffset, iflowface, iauxiflowface
+    integer(I4B) :: ioffset, iflowface, iauxiflowface, iface
     real(DP) :: qbnd
     character(len=LENAUXNAME) :: auxname
     integer(I4B) :: naux
@@ -176,6 +181,7 @@ subroutine accumulate_flows(this)
     this%SourceFlows = DZERO
     this%SinkFlows = DZERO
     this%BoundaryFlows = DZERO
+    this%BoundaryFaces = 0
     do ip = 1, this%nflowpack
       iauxiflowface = 0
       naux = this%gwfpackages(ip)%naux
@@ -194,16 +200,19 @@ subroutine accumulate_flows(this)
         if (this%ibound(i) <= 0) cycle
         qbnd = this%gwfpackages(ip)%get_flow(ib)
         ! todo, after initial release: default iflowface values for different packages
-        iflowface = 0
+        iflowface = 0 ! iflowface number
+        iface = 0 ! internal face number
         if (iauxiflowface > 0) then
           iflowface = NINT(this%gwfpackages(ip)%auxvar(iauxiflowface, ib))
+          iface = iflowface
           ! maps bot -2 -> max_faces - 1, top -1 -> max_faces
-          if (iflowface < 0) iflowface = iflowface + this%max_faces + 1
+          if (iface < 0) iface = iface + this%max_faces + 1
         end if
-        if (iflowface .gt. 0) then
+        if (iface > 0) then
+          call this%mark_boundary_face(i, iface)
           ioffset = (i - 1) * this%max_faces
-          this%BoundaryFlows(ioffset + iflowface) = &
-            this%BoundaryFlows(ioffset + iflowface) + qbnd
+          this%BoundaryFlows(ioffset + iface) = &
+            this%BoundaryFlows(ioffset + iface) + qbnd
         else if (qbnd .gt. DZERO) then
           this%SourceFlows(i) = this%SourceFlows(i) + qbnd
         else if (qbnd .lt. DZERO) then
@@ -214,4 +223,44 @@ subroutine accumulate_flows(this)
 
   end subroutine accumulate_flows
 
+  !> @brief Mark a face as a boundary face.
+  subroutine mark_boundary_face(this, ic, iface)
+    class(PrtFmiType) :: this
+    integer(I4B), intent(in) :: ic !< node number (reduced)
+    integer(I4B), intent(in) :: iface !< face number
+    ! local
+    integer(I4B) :: bit_pos
+
+    if (ic <= 0 .or. ic > size(this%BoundaryFaces)) return
+    if (iface == 0) return
+    bit_pos = iface - 1 ! bit position 0-based
+    if (bit_pos < 0 .or. bit_pos > 31) then
+      print *, 'Invalid bitmask position: ', iface
+      print *, 'Expected a value in range [0, 31]'
+      call pstop(1)
+    end if
+    this%BoundaryFaces(ic) = ibset(this%BoundaryFaces(ic), bit_pos)
+  end subroutine mark_boundary_face
+
+  !> @brief Check if a face is assigned to a boundary package.
+  function is_boundary_face(this, ic, iface) result(is_boundary)
+    class(PrtFmiType) :: this
+    integer(I4B), intent(in) :: ic !< node number (reduced)
+    integer(I4B), intent(in) :: iface !< face number
+    logical(LGP) :: is_boundary
+    ! local
+    integer(I4B) :: bit_pos
+
+    is_boundary = .false.
+    if (ic <= 0 .or. ic > size(this%BoundaryFaces)) return
+    if (iface == 0) return
+    bit_pos = iface - 1 ! bit position 0-based
+    if (bit_pos < 0 .or. bit_pos > 31) then
+      print *, 'Invalid bitmask position: ', iface
+      print *, 'Expected a value in range [0, 31]'
+      call pstop(1)
+    end if
+    is_boundary = btest(this%BoundaryFaces(ic), bit_pos)
+  end function is_boundary_face
+
 end module PrtFmiModule