support releasetime array

autotest/prt/test_prt_release_timing.py

@@ -6,7 +6,9 @@
 the first and only time step of the first and
 only stress period, with duration 1 time unit,
 so the same value of 0.5 can be used for both
-RELEASETIME and FRACTION.
+RELEASETIME and FRACTION. A third test case
+checks that multiple values can be provided
+for RELEASETIME.
 
 Period-block FRACTION should work with FIRST
 and ALL, but flopy hangs with either option.
@@ -47,10 +49,11 @@
 
 from framework import TestFramework
 
-simname = "prtfmi05"
+simname = "prtrelt"
 cases = [
     # options block options
-    f"{simname}relt",  # RELEASETIME 0.5
+    f"{simname}sgl",  # RELEASETIME 0.5
+    f"{simname}dbl",  # RELEASETIME 0.5 0.6
     # period block options
     # f"{simname}all",  # ALL FRACTION 0.5      # todo debug flopy hanging
     # f"{simname}frst", # FIRST FRACTION 0.5    # todo debug flopy hanging
@@ -59,7 +62,7 @@
 
 
 def get_perioddata(name, periods=1, fraction=None) -> Optional[dict]:
-    if "relt" in name:
+    if "sgl" in name or "dbl" in name:
         return None
     opt = [
         "FIRST"
@@ -125,7 +128,13 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6, fraction=None):
     prp_track_csv_file = f"{prt_name}.prp.trk.csv"
     pdat = get_perioddata(prt_name, fraction=fraction)
     # fraction 0.5 equiv. to release time 0.5 since 1 period 1 step with length 1
-    trelease = fraction if "relt" in prt_name else None
+    releasetime = (
+        [fraction]
+        if "sgl" in prt_name
+        else [fraction, fraction + 0.1]
+        if "dbl" in prt_name
+        else None
+    )
     flopy.mf6.ModflowPrtprp(
         prt,
         pname="prp1",
@@ -135,7 +144,7 @@ def build_prt_sim(name, gwf_ws, prt_ws, mf6, fraction=None):
         perioddata=pdat,
         track_filerecord=[prp_track_file],
         trackcsv_filerecord=[prp_track_csv_file],
-        releasetime=trelease,
+        releasetimerecord=releasetime,
     )
 
     # create output control package
@@ -364,8 +373,12 @@ def check_output(idx, test, fraction):
     del mp7_pls["zloc"]
 
     # compare mf6 / mp7 pathline data
-    assert mf6_pls.shape == mp7_pls.shape
-    assert np.allclose(mf6_pls, mp7_pls, atol=1e-3)
+    if "dbl" in name:
+        assert len(mf6_pls) == 2 * len(mp7_pls)
+        # todo check for double mass
+    else:
+        assert mf6_pls.shape == mp7_pls.shape
+        assert np.allclose(mf6_pls, mp7_pls, atol=1e-3)
 
 
 @pytest.mark.parametrize("idx, name", enumerate(cases))

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

@@ -127,13 +127,38 @@ optional true
 longname drape
 description is a text keyword to indicate that if a particle's release point is in a cell that happens to be dry at release time, the particle is to be moved to the topmost active cell below it, if any. By default, a particle is not released into the simulation if its release point's cell is dry at release time, instead the particle is terminated immediately with ireason 3 and istatus 8.
 
+block options
+name releasetimerecord
+type record releasetime times
+shape
+reader urword
+tagged true
+optional true
+longname
+description
+
 block options
 name releasetime
+type keyword
+reader urword
+in_record true
+optional false
+tagged true
+shape
+longname
+description keyword indicating release times will follow
+
+block options
+name times
 type double precision
+shape (unknown)
 reader urword
-optional true
-longname global release time
-description real value defining the time at which to release particles. This is comparable to MODPATH 7 referencetime option 1. Overrides release settings specified in period data.
+in_record true
+optional false
+tagged false
+repeating true
+longname release times
+description times to release, relative to the beginning of the simulation.
 
 # --------------------- prt prp dimensions ---------------------
 

doc/mf6io/mf6ivar/md/mf6ivar.md

@@ -1225,7 +1225,8 @@
 | PRT | PRP | OPTIONS | STOP_AT_WEAK_SINK | KEYWORD | is a text keyword to indicate that a particle is to terminate when it enters a cell that is a weak sink.  By default, particles are allowed to pass though cells that are weak sinks. |
 | PRT | PRP | OPTIONS | ISTOPZONE | INTEGER | integer value defining the stop zone number.  If cells have been assigned IZONE values in the GRIDDATA block, a particle terminates if it enters a cell whose IZONE value matches ISTOPZONE.  An ISTOPZONE value of zero indicates that there is no stop zone.  The default value is zero. |
 | PRT | PRP | OPTIONS | DRAPE | KEYWORD | is a text keyword to indicate that if a particle's release point is in a cell that happens to be dry at release time, the particle is to be moved to the topmost active cell below it, if any. By default, a particle is not released into the simulation if its release point's cell is dry at release time, instead the particle is terminated immediately with ireason 3 and istatus 8. |
-| PRT | PRP | OPTIONS | RELEASETIME | DOUBLE PRECISION | real value defining the time at which to release particles. This is comparable to MODPATH 7 referencetime option 1. Overrides release settings specified in period data. |
+| PRT | PRP | OPTIONS | RELEASETIME | KEYWORD | keyword indicating release times will follow |
+| PRT | PRP | OPTIONS | TIMES | DOUBLE PRECISION (UNKNOWN) | times to release, relative to the beginning of the simulation. |
 | PRT | PRP | DIMENSIONS | NRELEASEPTS | INTEGER | is the number of particle release points. |
 | PRT | PRP | PACKAGEDATA | IRPTNO | INTEGER | integer value that defines the PRP release point number associated with the specified PACKAGEDATA data on the line. IRPTNO must be greater than zero and less than or equal to NRELEASEPTS.  The program will terminate with an error if information for a PRP release point number is specified more than once. |
 | PRT | PRP | PACKAGEDATA | CELLID | INTEGER (NCELLDIM) | is the cell identifier, and depends on the type of grid that is used for the simulation.  For a structured grid that uses the DIS input file, CELLID is the layer, row, and column.   For a grid that uses the DISV input file, CELLID is the layer and CELL2D number.  If the model uses the unstructured discretization (DISU) input file, CELLID is the node number for the cell. |

doc/mf6io/mf6ivar/tex/gwf-maw-desc.tex

@@ -93,7 +93,7 @@
 
 \item \texttt{scrn\_top}---value that defines the top elevation of the screen for the multi-aquifer well connection. If CONDEQN is SPECIFIED, THIEM, SKIN, or COMPOSITE, SCRN\_TOP can be any value and is set to the top of the cell. If CONDEQN is MEAN, SCRN\_TOP is set to the multi-aquifer well connection cell top if the specified value is greater than the cell top. The program will terminate with an error if the screen top is less than the screen bottom.
 
-\item \texttt{scrn\_bot}---value that defines the bottom elevation of the screen for the multi-aquifer well connection. If CONDEQN is SPECIFIED, THIEM, SKIN, or COMPOSITE, SCRN\_BOT can be any value is set to the bottom of the cell. If CONDEQN is MEAN, SCRN\_BOT is set to the multi-aquifer well connection cell bottom if the specified value is less than the cell bottom. The program will terminate with an error if the screen bottom is greater than the screen top.
+\item \texttt{scrn\_bot}---value that defines the bottom elevation of the screen for the multi-aquifer well connection. If CONDEQN is SPECIFIED, THIEM, SKIN, or COMPOSITE, SCRN\_BOT can be any value and is set to the bottom of the cell. If CONDEQN is MEAN, SCRN\_BOT is set to the multi-aquifer well connection cell bottom if the specified value is less than the cell bottom. The program will terminate with an error if the screen bottom is greater than the screen top.
 
 \item \texttt{hk\_skin}---value that defines the skin (filter pack) hydraulic conductivity (if CONDEQN for the multi-aquifer well is SKIN, CUMULATIVE, or MEAN) or conductance (if CONDEQN for the multi-aquifer well is SPECIFIED) for each GWF node connected to the multi-aquifer well (NGWFNODES). If CONDEQN is SPECIFIED, HK\_SKIN must be greater than or equal to zero.  HK\_SKIN can be any value if CONDEQN is THIEM. Otherwise, HK\_SKIN must be greater than zero. If CONDEQN is SKIN, the contrast between the cell transmissivity (the product of geometric mean horizontal hydraulic conductivity and the cell thickness) and the well transmissivity (the product of HK\_SKIN and the screen thicknesses) must be greater than one in node CELLID or the program will terminate with an error condition; if an error condition occurs, it is suggested that the HK\_SKIN be reduced to a value less than K11 and K22 in node CELLID or the THEIM or MEAN conductance equations be used for these multi-aquifer wells.
 

doc/mf6io/mf6ivar/tex/gwf-rcha-desc.tex

@@ -21,7 +21,7 @@
 
 \item \texttt{FILEIN}---keyword to specify that an input filename is expected next.
 
-\item \texttt{tas6\_filename}---defines a time-array-series file defining a time-array series that can be used to assign time-varying values. See the Time-Variable Input section for instructions on using the time-array series capability.
+\item \texttt{tas6\_filename}---defines a time-array-series file defining a time-array series that can be used to assign time-varying values. See the ``Time-Variable Input'' section for instructions on using the time-array series capability.
 
 \item \texttt{OBS6}---keyword to specify that record corresponds to an observations file.
 
@@ -35,7 +35,7 @@
 
 \item \texttt{irch}---IRCH is the layer number that defines the layer in each vertical column where recharge is applied. If IRCH is omitted, recharge by default is applied to cells in layer 1.  IRCH can only be used if READASARRAYS is specified in the OPTIONS block.  If IRCH is specified, it must be specified as the first variable in the PERIOD block or MODFLOW will terminate with an error.
 
-\item \textcolor{blue}{\texttt{recharge}---is the recharge flux rate ($LT^{-1}$).  This rate is multiplied inside the program by the surface area of the cell to calculate the volumetric recharge rate. The recharge array may be defined by a time-array series (see the ``Using Time-Array Series in a Package'' section).}
+\item \textcolor{blue}{\texttt{recharge}---is the recharge flux rate ($LT^{-1}$).  This rate is multiplied inside the program by the surface area of the cell to calculate the volumetric recharge rate. The recharge array may be defined by a time-array series (see the "Using Time-Array Series in a Package" section).}
 
 \item \textcolor{blue}{\texttt{aux}---is an array of values for auxiliary variable aux(iaux), where iaux is a value from 1 to naux, and aux(iaux) must be listed as part of the auxiliary variables.  A separate array can be specified for each auxiliary variable.  If an array is not specified for an auxiliary variable, then a value of zero is assigned.  If the value specified here for the auxiliary variable is the same as auxmultname, then the recharge array will be multiplied by this array.}
 

doc/mf6io/mf6ivar/tex/prt-prp-desc.tex

@@ -25,7 +25,9 @@
 
 \item \texttt{DRAPE}---is a text keyword to indicate that if a particle's release point is in a cell that happens to be dry at release time, the particle is to be moved to the topmost active cell below it, if any. By default, a particle is not released into the simulation if its release point's cell is dry at release time, instead the particle is terminated immediately with ireason 3 and istatus 8.
 
-\item \texttt{releasetime}---real value defining the time at which to release particles. This is comparable to MODPATH 7 referencetime option 1. Overrides release settings specified in period data.
+\item \texttt{RELEASETIME}---keyword indicating release times will follow
+
+\item \texttt{times}---times to release, relative to the beginning of the simulation.
 
 \end{description}
 \item \textbf{Block: DIMENSIONS}

doc/mf6io/mf6ivar/tex/prt-prp-options.dat

@@ -7,5 +7,5 @@ BEGIN OPTIONS
   [STOP_AT_WEAK_SINK]
   [ISTOPZONE <istopzone>]
   [DRAPE]
-  [RELEASETIME <releasetime>]
+  [RELEASETIME <times(unknown)>]
 END OPTIONS

src/Model/ModelUtilities/TrackData.f90

@@ -7,9 +7,18 @@ module TrackModule
   implicit none
 
   private save_record
-  public :: TrackControlType
+  public :: TrackFileControlType
   public :: TrackFileType
 
+  !> @brief
+  ! type :: TrackTimesType
+  !   real(DP), dimension(:), pointer, contiguous, public :: times => null()
+  !   integer(I4B) :: current = 0
+  ! contains
+  !   procedure :: init
+  !   procedure :: expand
+  ! end type TrackTimesType
+
   !> @brief Output file containing all or some particle pathlines.
   !!
   !! May be associated with a particle release point (PRP) package
@@ -37,7 +46,7 @@ module TrackModule
   !! An arbitrary number of files can be managed, internal
   !! arrays are resized as needed.
   !<
-  type :: TrackControlType
+  type :: TrackFileControlType
     private
     type(TrackFileType), public, allocatable :: trackfiles(:) !< output files
     integer(I4B), public :: ntrackfiles !< number of output files
@@ -52,7 +61,7 @@ module TrackModule
     procedure, public :: init_track_file
     procedure, public :: save
     procedure, public :: set_track_events
-  end type TrackControlType
+  end type TrackFileControlType
 
   ! Track file headers
   character(len=*), parameter, public :: TRACKHEADERS = &
@@ -123,7 +132,7 @@ module TrackModule
   !> @brief Initialize a new track file
   subroutine init_track_file(this, iun, csv, iprp)
     ! -- dummy
-    class(TrackControlType) :: this
+    class(TrackFileControlType) :: this
     integer(I4B), intent(in) :: iun
     logical(LGP), intent(in), optional :: csv
     integer(I4B), intent(in), optional :: iprp
@@ -152,7 +161,7 @@ end subroutine init_track_file
   !> @brief Expand the trackfile array, internal use only
   subroutine expand(this, increment)
     ! -- dummy
-    class(TrackControlType) :: this
+    class(TrackFileControlType) :: this
     integer(I4B), optional, intent(in) :: increment
     ! -- local
     integer(I4B) :: inclocal, isize, newsize
@@ -249,7 +258,7 @@ subroutine save_record(iun, particle, kper, kstp, reason, csv)
   !<
   subroutine save(this, particle, kper, kstp, reason, level)
     ! -- dummy
-    class(TrackControlType), intent(inout) :: this
+    class(TrackFileControlType), intent(inout) :: this
     type(ParticleType), pointer, intent(in) :: particle
     integer(I4B), intent(in) :: kper, kstp
     integer(I4B), intent(in) :: reason
@@ -300,7 +309,7 @@ subroutine set_track_events(this, &
                               terminate, &
                               weaksink, &
                               tracktimes)
-    class(TrackControlType) :: this
+    class(TrackFileControlType) :: this
     logical(LGP), intent(in) :: release, transit, timestep, &
                                 terminate, weaksink, tracktimes
     this%itrackrelease = release

src/Model/ParticleTracking/prt1.f90

@@ -19,7 +19,7 @@ module PrtModule
   use BudgetModule, only: BudgetType
   use ListModule, only: ListType
   use ParticleModule, only: ParticleType, create_particle
-  use TrackModule, only: TrackControlType, TrackFileType
+  use TrackModule, only: TrackFileControlType, TrackFileType
   use SimModule, only: count_errors, store_error, store_error_filename
   use MemoryManagerModule, only: mem_allocate
   use MethodModule, only: MethodType
@@ -44,7 +44,7 @@ module PrtModule
     type(PrtObsType), pointer :: obs => null() ! observation package
     type(BudgetType), pointer :: budget => null() ! budget object
     class(MethodType), pointer :: method => null() ! tracking method
-    type(TrackControlType), pointer :: trackctl ! track file control
+    type(TrackFileControlType), pointer :: trackctl ! track file control
     integer(I4B), pointer :: infmi => null() ! unit number FMI
     integer(I4B), pointer :: inmip => null() ! unit number MIP
     integer(I4B), pointer :: inmvt => null() ! unit number MVT
@@ -289,12 +289,12 @@ subroutine prt_ar(this)
     if (this%oc%itrkcsv > 0) &
       call this%trackctl%init_track_file(this%oc%itrkcsv, csv=.true.)
     call this%trackctl%set_track_events( &
-      this%oc%itrkrls, &
-      this%oc%itrktrs, &
-      this%oc%itrktst, &
-      this%oc%itrkter, &
-      this%oc%itrkwsk, &
-      this%oc%itrktrt)
+      this%oc%trackrelease, &
+      this%oc%tracktransit, &
+      this%oc%tracktimestep, &
+      this%oc%trackterminate, &
+      this%oc%trackweaksink, &
+      this%oc%tracktime)
   end subroutine prt_ar
 
   !> @brief Read and prepare (calls package read and prepare routines)