feat(PRT): local z option for particle release point package (#1756)

Optional support for local z coordinates. This allows positioning particles based on the water table. The motivations here are 1) modpath 7 parity and 2) to provide the same capabilities whether prt is running alongside a gwf model or as a post-processor. Previously dynamic positioning was possible only if prt runs as a post-processor.

autotest/test_prt_stop_zones.py

@@ -54,6 +54,25 @@ def create_izone(nlay, nrow, ncol):
     return izone
 
 
+def build_gwf_sim(name, ws, mf6):
+    gwf_sim = FlopyReadmeCase.get_gwf_sim(
+        name, ws, mf6
+    )
+    gwf = gwf_sim.get_model()
+    dis = gwf.get_package("DIS")
+    nlay = int(name[-1])
+    botm = [FlopyReadmeCase.top - (k + 1) for k in range(nlay)]
+    botm_data = np.array(
+        [
+            list(repeat(b, FlopyReadmeCase.nrow * FlopyReadmeCase.ncol))
+            for b in botm
+        ]
+    ).reshape((nlay, FlopyReadmeCase.nrow, FlopyReadmeCase.ncol))
+    dis.nlay = nlay
+    dis.botm.set_data(botm_data)
+    return gwf_sim
+
+
 def build_prt_sim(name, gwf_ws, prt_ws, mf6):
     # create simulation
     sim = flopy.mf6.MFSimulation(
@@ -193,21 +212,8 @@ def build_mp7_sim(name, ws, mp7, gwf):
 
 
 def build_models(idx, test):
-    gwf_sim = FlopyReadmeCase.get_gwf_sim(
-        test.name, test.workspace, test.targets["mf6"]
-    )
+    gwf_sim = build_gwf_sim(test.name, test.workspace, test.targets["mf6"])
     gwf = gwf_sim.get_model()
-    dis = gwf.get_package("DIS")
-    nlay = int(test.name[-1])
-    botm = [FlopyReadmeCase.top - (k + 1) for k in range(nlay)]
-    botm_data = np.array(
-        [
-            list(repeat(b, FlopyReadmeCase.nrow * FlopyReadmeCase.ncol))
-            for b in botm
-        ]
-    ).reshape((nlay, FlopyReadmeCase.nrow, FlopyReadmeCase.ncol))
-    dis.nlay = nlay
-    dis.botm.set_data(botm_data)
     prt_sim = build_prt_sim(
         test.name, test.workspace, test.workspace / "prt", test.targets["mf6"]
     )

doc/mf6io/mf6ivar/dfn/prt-prp.dfn

@@ -10,6 +10,14 @@ optional true
 longname
 description REPLACE boundnames {'{#1}': 'release-point'}
 
+block options
+name local_z
+type keyword
+reader urword
+optional true
+longname whether to use local z coordinates
+description indicates that rptz defines the local z coordinate of the release point within the cell, with value between 0 and 1 (inclusive).  If the cell is unsaturated at release time, the z coordinate is computed relative to the water table rather than the top of the cell.
+
 block options
 name track_filerecord
 type record track fileout trackfile

src/Model/ParticleTracking/prt-prp.f90

@@ -46,6 +46,7 @@ module PrtPrpModule
     integer(I4B), pointer :: itrkhdr => null() !< track header file
     integer(I4B), pointer :: itrkcsv => null() !< CSV track file
     integer(I4B), pointer :: irlstls => null() !< release time file
+    logical(LGP), pointer :: localz => null() !< compute z coordinates local to the cell
     logical(LGP), pointer :: rlsall => null() !< release in all time step
     logical(LGP), pointer :: rlsfirst => null() !< release in first time step
     logical(LGP), pointer :: rlstimelist => null() !< use global release time
@@ -58,6 +59,7 @@ module PrtPrpModule
     real(DP), pointer, contiguous :: rptx(:) => null() !< release point x coordinates
     real(DP), pointer, contiguous :: rpty(:) => null() !< release point y coordinates
     real(DP), pointer, contiguous :: rptz(:) => null() !< release point z coordinates
+    real(DP), pointer, contiguous :: locz(:) => null() !< release point local z coordinates
     real(DP), pointer, contiguous :: rptmass(:) => null() !< total mass released from point
     character(len=LENBOUNDNAME), pointer, contiguous :: rptname(:) => null() !< release point names
     type(TimeSelectType), pointer :: releasetimes
@@ -140,6 +142,7 @@ subroutine prp_da(this)
     call mem_deallocate(this%rlsall)
     call mem_deallocate(this%rlsfirst)
     call mem_deallocate(this%rlstimelist)
+    call mem_deallocate(this%localz)
     call mem_deallocate(this%offset)
     call mem_deallocate(this%stoptime)
     call mem_deallocate(this%stoptraveltime)
@@ -157,6 +160,7 @@ subroutine prp_da(this)
     call mem_deallocate(this%rptx)
     call mem_deallocate(this%rpty)
     call mem_deallocate(this%rptz)
+    call mem_deallocate(this%locz)
     call mem_deallocate(this%rptnode)
     call mem_deallocate(this%rptmass)
     call mem_deallocate(this%rlskstp)
@@ -201,6 +205,7 @@ subroutine prp_allocate_arrays(this, nodelist, auxvar)
     call mem_allocate(this%rptx, this%nreleasepts, 'RPTX', this%memoryPath)
     call mem_allocate(this%rpty, this%nreleasepts, 'RPTY', this%memoryPath)
     call mem_allocate(this%rptz, this%nreleasepts, 'RPTZ', this%memoryPath)
+    call mem_allocate(this%locz, this%nreleasepts, 'LOCZ', this%memoryPath)
     call mem_allocate(this%rptmass, this%nreleasepts, 'RPTMASS', this%memoryPath)
     call mem_allocate(this%rptnode, this%nreleasepts, 'RPTNODER', &
                       this%memoryPath)
@@ -230,6 +235,7 @@ subroutine prp_allocate_scalars(this)
     call mem_allocate(this%rlsall, 'RLSALL', this%memoryPath)
     call mem_allocate(this%rlsfirst, 'RLSFIRST', this%memoryPath)
     call mem_allocate(this%rlstimelist, 'RELEASETIME', this%memoryPath)
+    call mem_allocate(this%localz, 'LOCALZ', this%memoryPath)
     call mem_allocate(this%offset, 'OFFSET', this%memoryPath)
     call mem_allocate(this%stoptime, 'STOPTIME', this%memoryPath)
     call mem_allocate(this%stoptraveltime, 'STOPTRAVELTIME', this%memoryPath)
@@ -247,6 +253,7 @@ subroutine prp_allocate_scalars(this)
     this%rlsall = .false.
     this%rlsfirst = .false.
     this%rlstimelist = .false.
+    this%localz = .false.
     this%offset = DZERO
     this%stoptime = huge(1d0)
     this%stoptraveltime = huge(1d0)
@@ -303,7 +310,7 @@ subroutine prp_ad(this)
     character(len=LINELENGTH) :: errmsg
     integer(I4B) :: ic, icu, nps, nts, nrel, &
                     nreleasets, np, irow, icol, ilay, icpl
-    real(DP) :: x, y, z, trelease, tend
+    real(DP) :: x, y, z, trelease, tend, top, bot, hds
     real(DP), allocatable :: polyverts(:, :)
     type(ParticleType), pointer :: particle
 
@@ -418,7 +425,14 @@ subroutine prp_ad(this)
         end if
         particle%x = x
         particle%y = y
-        particle%z = this%rptz(nps)
+        if (this%localz) then
+          top = this%fmi%dis%top(ic)
+          bot = this%fmi%dis%bot(ic)
+          hds = this%fmi%gwfhead(ic)
+          particle%z = bot + this%rptz(nps) * (hds - bot)
+        else
+          particle%z = this%rptz(nps)
+        end if
         particle%trelease = trelease
         ! Set stopping time to earlier of times specified by STOPTIME and STOPTRAVELTIME
         if (this%stoptraveltime == huge(1d0)) then
@@ -794,6 +808,8 @@ subroutine prp_options(this, option, found)
           &FOLLOWED BY FILEOUT')
       end if
       found = .true.
+    case ('LOCAL_Z')
+      this%localz = .true.
     case default
       found = .false.
     end select
@@ -873,6 +889,11 @@ subroutine prp_read_packagedata(this)
         y(n) = this%parser%GetDouble()
         z(n) = this%parser%GetDouble()
 
+        if (this%localz .and. (z(n) < 0 .or. z(n) > 1)) then
+          call store_error('Local z coordinate must fall in the interval [0, 1]')
+          cycle
+        end if
+
         ! -- set default boundname
         write (cno, '(i9.9)') n
         bndName = 'PRP'//cno