MODFLOW 6  version 6.7.0.dev3
USGS Modular Hydrologic Model
prtprpmodule Module Reference

Data Types

type  prtprptype
 Particle release point (PRP) package. More...
 

Functions/Subroutines

subroutine, public prp_create (packobj, id, ibcnum, inunit, iout, namemodel, pakname, input_mempath, fmi)
 Create a new particle release point package. More...
 
subroutine prp_da (this)
 Deallocate memory. More...
 
subroutine prp_set_pointers (this, ibound, izone)
 @ brief Set pointers to model variables More...
 
subroutine prp_allocate_arrays (this, nodelist, auxvar)
 Allocate arrays. More...
 
subroutine prp_allocate_scalars (this)
 Allocate scalars. More...
 
subroutine prp_ar (this)
 @ brief Allocate and read period data More...
 
subroutine prp_ad (this)
 Advance a time step and release particles if scheduled. More...
 
subroutine log_release (this)
 Log the release scheduled for this time step. More...
 
subroutine validate_release_point (this, ic, x, y, z)
 Verify that the release point is in the cell. More...
 
subroutine release (this, ip, trelease)
 Release a particle at the specified time. More...
 
subroutine initialize_particle (this, particle, ip, trelease)
 
subroutine prp_rp (this)
 @ brief Read and prepare period data for particle input More...
 
subroutine prp_cq_simrate (this, hnew, flowja, imover)
 @ brief Calculate flow between package and model. More...
 
subroutine define_listlabel (this)
 
logical function prp_obs_supported (this)
 Indicates whether observations are supported. More...
 
subroutine prp_df_obs (this)
 Store supported observations. More...
 
subroutine prp_options (this)
 @ brief Set options specific to PrtPrpType More...
 
subroutine prp_log_options (this, found, trackfile, trackcsvfile)
 @ brief Log options specific to PrtPrpType More...
 
subroutine prp_dimensions (this)
 @ brief Set dimensions specific to PrtPrpType More...
 
subroutine prp_packagedata (this)
 Load package data (release points). More...
 
subroutine prp_releasetimes (this)
 Load explicitly specified release times. More...
 
subroutine prp_load_releasetimefrequency (this)
 Load regularly spaced release times if configured. More...
 

Variables

character(len=lenftype) ftype = 'PRP'
 
character(len=16) text = ' PRP'
 

Function/Subroutine Documentation

◆ define_listlabel()

subroutine prtprpmodule::define_listlabel ( class(prtprptype), intent(inout)  this)

Definition at line 643 of file prt-prp.f90.

644  class(PrtPrpType), intent(inout) :: this
645  ! not implemented, not used

◆ initialize_particle()

subroutine prtprpmodule::initialize_particle ( class(prtprptype), intent(inout)  this,
type(particletype), intent(inout), pointer  particle,
integer(i4b), intent(in)  ip,
real(dp), intent(in)  trelease 
)
private
Parameters
[in,out]thisthis instance
[in,out]particlethe particle
[in]ipparticle index
[in]treleaserelease time

Definition at line 456 of file prt-prp.f90.

458  class(PrtPrpType), intent(inout) :: this !< this instance
459  type(ParticleType), pointer, intent(inout) :: particle !< the particle
460  integer(I4B), intent(in) :: ip !< particle index
461  real(DP), intent(in) :: trelease !< release time
462  ! local
463  integer(I4B) :: irow, icol, ilay, icpl
464  integer(I4B) :: ic, icu, ic_old
465  real(DP) :: x, y, z
466  real(DP) :: top, bot, hds
467 
468  ic = this%rptnode(ip)
469  icu = this%dis%get_nodeuser(ic)
470 
471  call create_particle(particle)
472 
473  if (size(this%boundname) /= 0) then
474  particle%name = this%boundname(ip)
475  else
476  particle%name = ''
477  end if
478 
479  particle%irpt = ip
480  particle%istopweaksink = this%istopweaksink
481  particle%istopzone = this%istopzone
482  particle%idrymeth = this%idrymeth
483  particle%icu = icu
484 
485  select type (dis => this%dis)
486  type is (distype)
487  call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, irow, icol, ilay)
488  type is (disvtype)
489  call get_jk(icu, dis%ncpl, dis%nlay, icpl, ilay)
490  end select
491  particle%ilay = ilay
492  particle%izone = this%rptzone(ic)
493  particle%istatus = 0 ! status 0 until tracking starts
494  ! If the cell is inactive, either drape the particle
495  ! to the top-most active cell beneath it if drape is
496  ! enabled, or else terminate permanently unreleased.
497  if (this%ibound(ic) == 0) then
498  ic_old = ic
499  if (this%drape) then
500  call this%dis%highest_active(ic, this%ibound)
501  if (ic == ic_old .or. this%ibound(ic) == 0) then
502  ! negative unreleased status signals to the
503  ! tracking method that we haven't yet saved
504  ! a termination record, it needs to do so.
505  particle%istatus = -1 * term_unreleased
506  end if
507  else
508  particle%istatus = -1 * term_unreleased
509  end if
510  end if
511 
512  ! Load coordinates and transform if needed
513  x = this%rptx(ip)
514  y = this%rpty(ip)
515  if (this%localz) then
516  top = this%fmi%dis%top(ic)
517  bot = this%fmi%dis%bot(ic)
518  hds = this%fmi%gwfhead(ic)
519  z = bot + this%rptz(ip) * (hds - bot)
520  else
521  z = this%rptz(ip)
522  end if
523 
524  if (this%ichkmeth > 0) &
525  call this%validate_release_point(ic, x, y, z)
526 
527  particle%x = x
528  particle%y = y
529  particle%z = z
530  particle%trelease = trelease
531 
532  ! Set stop time to earlier of STOPTIME and STOPTRAVELTIME
533  if (this%stoptraveltime == huge(1d0)) then
534  particle%tstop = this%stoptime
535  else
536  particle%tstop = particle%trelease + this%stoptraveltime
537  if (this%stoptime < particle%tstop) particle%tstop = this%stoptime
538  end if
539 
540  particle%ttrack = particle%trelease
541  particle%itrdomain(level_model) = 0
542  particle%iboundary(level_model) = 0
543  particle%itrdomain(level_feature) = ic
544  particle%iboundary(level_feature) = 0
545  particle%itrdomain(level_subfeature) = 0
546  particle%iboundary(level_subfeature) = 0
547  particle%frctrn = this%frctrn
548  particle%iexmeth = this%iexmeth
549  particle%extend = this%extend
550  particle%icycwin = this%icycwin
551  particle%extol = this%extol
@ term_unreleased
terminated permanently unreleased
Definition: Particle.f90:35
Here is the call graph for this function:

◆ log_release()

subroutine prtprpmodule::log_release ( class(prtprptype), intent(inout)  this)
Parameters
[in,out]thisprp

Definition at line 380 of file prt-prp.f90.

381  class(PrtPrpType), intent(inout) :: this !< prp
382  if (this%iprpak > 0) then
383  write (this%iout, "(1x,/1x,a,1x,i0)") &
384  'PARTICLE RELEASE FOR PRP', this%ibcnum
385  call this%schedule%log(this%iout)
386  end if

◆ prp_ad()

subroutine prtprpmodule::prp_ad ( class(prtprptype this)
private

Definition at line 319 of file prt-prp.f90.

320  use tdismodule, only: totalsimtime
321  class(PrtPrpType) :: this
322  integer(I4B) :: ip, it
323  real(DP) :: t
324 
325  ! Notes
326  ! -----
327  ! Each release point can be thought of as
328  ! a gumball machine with infinite supply:
329  ! a point can release an arbitrary number
330  ! of particles, but only one at any time.
331  ! Coincident release times are merged to
332  ! a single time by the release scheduler.
333 
334  ! Reset mass accumulators for this time step.
335  do ip = 1, this%nreleasepoints
336  this%rptm(ip) = dzero
337  end do
338 
339  ! Advance the release schedule and check if
340  ! any releases will be made this time step.
341  call this%schedule%advance()
342  if (.not. this%schedule%any()) return
343 
344  ! Log the schedule to the list file.
345  call this%log_release()
346 
347  ! Expand the particle store. We know from the
348  ! schedule how many particles will be released.
349  call this%particles%resize( &
350  this%particles%num_stored() + &
351  (this%nreleasepoints * this%schedule%count()), &
352  this%memoryPath)
353 
354  ! Release a particle from each point for
355  ! each release time in the current step.
356  do ip = 1, this%nreleasepoints
357  do it = 1, this%schedule%count()
358  t = this%schedule%times(it)
359  ! Skip the release time if it's before the simulation
360  ! starts, or if no `extend_tracking`, after it ends.
361  if (t < dzero) then
362  write (warnmsg, '(a,g0,a)') &
363  'Skipping negative release time (t=', t, ').'
364  call store_warning(warnmsg)
365  cycle
366  else if (t > totalsimtime .and. .not. this%extend) then
367  write (warnmsg, '(a,g0,a)') &
368  'Skipping release time falling after the end of the &
369  &simulation (t=', t, '). Enable EXTEND_TRACKING to &
370  &release particles after the simulation end time.'
371  call store_warning(warnmsg)
372  cycle
373  end if
374  call this%release(ip, t)
375  end do
376  end do
real(dp), pointer, public totalsimtime
time at end of simulation
Definition: tdis.f90:37
Here is the call graph for this function:

◆ prp_allocate_arrays()

subroutine prtprpmodule::prp_allocate_arrays ( class(prtprptype this,
integer(i4b), dimension(:), optional, pointer, contiguous  nodelist,
real(dp), dimension(:, :), optional, pointer, contiguous  auxvar 
)
private

Definition at line 206 of file prt-prp.f90.

207  ! dummy
208  class(PrtPrpType) :: this
209  integer(I4B), dimension(:), pointer, contiguous, optional :: nodelist
210  real(DP), dimension(:, :), pointer, contiguous, optional :: auxvar
211  ! local
212  integer(I4B) :: nps
213 
214  call this%BndExtType%allocate_arrays()
215 
216  ! Allocate particle store, starting with the number
217  ! of release points (arrays resized if/when needed)
218  call create_particle_store( &
219  this%particles, &
220  this%nreleasepoints, &
221  this%memoryPath)
222 
223  ! Allocate arrays
224  call mem_allocate(this%rptx, this%nreleasepoints, 'RPTX', this%memoryPath)
225  call mem_allocate(this%rpty, this%nreleasepoints, 'RPTY', this%memoryPath)
226  call mem_allocate(this%rptz, this%nreleasepoints, 'RPTZ', this%memoryPath)
227  call mem_allocate(this%rptm, this%nreleasepoints, 'RPTMASS', &
228  this%memoryPath)
229  call mem_allocate(this%rptnode, this%nreleasepoints, 'RPTNODER', &
230  this%memoryPath)
231  call mem_allocate(this%rptname, lenboundname, this%nreleasepoints, &
232  'RPTNAME', this%memoryPath)
233 
234  ! Initialize arrays
235  do nps = 1, this%nreleasepoints
236  this%rptm(nps) = dzero
237  end do
Here is the call graph for this function:

◆ prp_allocate_scalars()

subroutine prtprpmodule::prp_allocate_scalars ( class(prtprptype this)
private

Definition at line 241 of file prt-prp.f90.

242  class(PrtPrpType) :: this
243 
244  ! Allocate parent's scalars
245  call this%BndExtType%allocate_scalars()
246 
247  ! Allocate scalars for this type
248  call mem_allocate(this%localz, 'LOCALZ', this%memoryPath)
249  call mem_allocate(this%extend, 'EXTEND', this%memoryPath)
250  call mem_allocate(this%offset, 'OFFSET', this%memoryPath)
251  call mem_allocate(this%stoptime, 'STOPTIME', this%memoryPath)
252  call mem_allocate(this%stoptraveltime, 'STOPTRAVELTIME', this%memoryPath)
253  call mem_allocate(this%istopweaksink, 'ISTOPWEAKSINK', this%memoryPath)
254  call mem_allocate(this%istopzone, 'ISTOPZONE', this%memoryPath)
255  call mem_allocate(this%drape, 'DRAPE', this%memoryPath)
256  call mem_allocate(this%idrymeth, 'IDRYMETH', this%memoryPath)
257  call mem_allocate(this%nreleasepoints, 'NRELEASEPOINTS', this%memoryPath)
258  call mem_allocate(this%nreleasetimes, 'NRELEASETIMES', this%memoryPath)
259  call mem_allocate(this%nparticles, 'NPARTICLES', this%memoryPath)
260  call mem_allocate(this%itrkout, 'ITRKOUT', this%memoryPath)
261  call mem_allocate(this%itrkhdr, 'ITRKHDR', this%memoryPath)
262  call mem_allocate(this%itrkcsv, 'ITRKCSV', this%memoryPath)
263  call mem_allocate(this%irlstls, 'IRLSTLS', this%memoryPath)
264  call mem_allocate(this%frctrn, 'FRCTRN', this%memoryPath)
265  call mem_allocate(this%iexmeth, 'IEXMETH', this%memoryPath)
266  call mem_allocate(this%ichkmeth, 'ICHKMETH', this%memoryPath)
267  call mem_allocate(this%icycwin, 'ICYCWIN', this%memoryPath)
268  call mem_allocate(this%extol, 'EXTOL', this%memoryPath)
269  call mem_allocate(this%rttol, 'RTTOL', this%memoryPath)
270  call mem_allocate(this%rtfreq, 'RTFREQ', this%memoryPath)
271 
272  ! Set values
273  this%localz = .false.
274  this%extend = .false.
275  this%offset = dzero
276  this%stoptime = huge(1d0)
277  this%stoptraveltime = huge(1d0)
278  this%istopweaksink = 0
279  this%istopzone = 0
280  this%drape = .false.
281  this%idrymeth = 0
282  this%nreleasepoints = 0
283  this%nreleasetimes = 0
284  this%nparticles = 0
285  this%itrkout = 0
286  this%itrkhdr = 0
287  this%itrkcsv = 0
288  this%irlstls = 0
289  this%frctrn = .false.
290  this%iexmeth = 0
291  this%ichkmeth = 1
292  this%icycwin = 0
293  this%extol = dem5
294  this%rttol = dsame * dep9
295  this%rtfreq = dzero
296 

◆ prp_ar()

subroutine prtprpmodule::prp_ar ( class(prtprptype), intent(inout)  this)
private

Definition at line 300 of file prt-prp.f90.

301  ! dummy variables
302  class(PrtPrpType), intent(inout) :: this
303  ! local variables
304  integer(I4B) :: n
305 
306  call this%obs%obs_ar()
307 
308  if (this%inamedbound /= 0) then
309  do n = 1, this%nreleasepoints
310  this%boundname(n) = this%rptname(n)
311  end do
312  end if
313  do n = 1, this%nreleasepoints
314  this%nodelist(n) = this%rptnode(n)
315  end do

◆ prp_cq_simrate()

subroutine prtprpmodule::prp_cq_simrate ( class(prtprptype this,
real(dp), dimension(:), intent(in)  hnew,
real(dp), dimension(:), intent(inout)  flowja,
integer(i4b), intent(in)  imover 
)
Parameters
[in,out]flowjaflow between package and model
[in]imoverflag indicating if the mover package is active

Definition at line 609 of file prt-prp.f90.

610  ! modules
611  use tdismodule, only: delt
612  ! dummy variables
613  class(PrtPrpType) :: this
614  real(DP), dimension(:), intent(in) :: hnew
615  real(DP), dimension(:), intent(inout) :: flowja !< flow between package and model
616  integer(I4B), intent(in) :: imover !< flag indicating if the mover package is active
617  ! local variables
618  integer(I4B) :: i
619  integer(I4B) :: node
620  integer(I4B) :: idiag
621  real(DP) :: rrate
622 
623  ! If no boundaries, skip flow calculations.
624  if (this%nbound <= 0) return
625 
626  ! Loop through each boundary calculating flow.
627  do i = 1, this%nbound
628  node = this%nodelist(i)
629  rrate = dzero
630  ! If cell is no-flow or constant-head, then ignore it.
631  if (node > 0) then
632  ! Calculate the flow rate into the cell.
633  idiag = this%dis%con%ia(node)
634  rrate = this%rptm(i) * (done / delt) ! reciprocal of tstp length
635  flowja(idiag) = flowja(idiag) + rrate
636  end if
637 
638  ! Save simulated value to simvals array.
639  this%simvals(i) = rrate
640  end do
real(dp), pointer, public delt
length of the current time step
Definition: tdis.f90:29

◆ prp_create()

subroutine, public prtprpmodule::prp_create ( class(bndtype), pointer  packobj,
integer(i4b), intent(in)  id,
integer(i4b), intent(in)  ibcnum,
integer(i4b), intent(in)  inunit,
integer(i4b), intent(in)  iout,
character(len=*), intent(in)  namemodel,
character(len=*), intent(in)  pakname,
character(len=*), intent(in)  input_mempath,
type(prtfmitype), pointer  fmi 
)

Definition at line 101 of file prt-prp.f90.

103  ! dummy
104  class(BndType), pointer :: packobj
105  integer(I4B), intent(in) :: id
106  integer(I4B), intent(in) :: ibcnum
107  integer(I4B), intent(in) :: inunit
108  integer(I4B), intent(in) :: iout
109  character(len=*), intent(in) :: namemodel
110  character(len=*), intent(in) :: pakname
111  character(len=*), intent(in) :: input_mempath
112  type(PrtFmiType), pointer :: fmi
113  ! local
114  type(PrtPrpType), pointer :: prpobj
115  ! formats
116  character(len=*), parameter :: fmtheader = &
117  "(1x, /1x, 'PRP PARTICLE RELEASE POINT PACKAGE', &
118  &' INPUT READ FROM MEMPATH: ', a, /)"
119 
120  ! allocate the object and assign values to object variables
121  allocate (prpobj)
122  packobj => prpobj
123 
124  ! create name and memory path
125  call packobj%set_names(ibcnum, namemodel, pakname, ftype, input_mempath)
126  prpobj%text = text
127 
128  ! allocate scalars
129  call prpobj%prp_allocate_scalars()
130 
131  ! initialize package
132  call packobj%pack_initialize()
133 
134  packobj%inunit = inunit
135  packobj%iout = iout
136  packobj%id = id
137  packobj%ibcnum = ibcnum
138  packobj%ncolbnd = 4
139  packobj%iscloc = 1
140 
141  ! store pointer to flow model interface
142  prpobj%fmi => fmi
143 
144  ! if prp is enabled, print a message identifying it
145  if (inunit > 0) write (iout, fmtheader) input_mempath
Here is the caller graph for this function:

◆ prp_da()

subroutine prtprpmodule::prp_da ( class(prtprptype this)
private

Definition at line 149 of file prt-prp.f90.

150  class(PrtPrpType) :: this
151 
152  ! Deallocate parent
153  call this%BndExtType%bnd_da()
154 
155  ! Deallocate scalars
156  call mem_deallocate(this%localz)
157  call mem_deallocate(this%extend)
158  call mem_deallocate(this%offset)
159  call mem_deallocate(this%stoptime)
160  call mem_deallocate(this%stoptraveltime)
161  call mem_deallocate(this%istopweaksink)
162  call mem_deallocate(this%istopzone)
163  call mem_deallocate(this%drape)
164  call mem_deallocate(this%idrymeth)
165  call mem_deallocate(this%nreleasepoints)
166  call mem_deallocate(this%nreleasetimes)
167  call mem_deallocate(this%nparticles)
168  call mem_deallocate(this%itrkout)
169  call mem_deallocate(this%itrkhdr)
170  call mem_deallocate(this%itrkcsv)
171  call mem_deallocate(this%irlstls)
172  call mem_deallocate(this%frctrn)
173  call mem_deallocate(this%iexmeth)
174  call mem_deallocate(this%ichkmeth)
175  call mem_deallocate(this%icycwin)
176  call mem_deallocate(this%extol)
177  call mem_deallocate(this%rttol)
178  call mem_deallocate(this%rtfreq)
179 
180  ! Deallocate arrays
181  call mem_deallocate(this%rptx)
182  call mem_deallocate(this%rpty)
183  call mem_deallocate(this%rptz)
184  call mem_deallocate(this%rptnode)
185  call mem_deallocate(this%rptm)
186  call mem_deallocate(this%rptname, 'RPTNAME', this%memoryPath)
187 
188  ! Deallocate objects
189  call this%particles%destroy(this%memoryPath)
190  call this%schedule%destroy()
191  deallocate (this%particles)
192  deallocate (this%schedule)

◆ prp_df_obs()

subroutine prtprpmodule::prp_df_obs ( class(prtprptype this)
private

Definition at line 655 of file prt-prp.f90.

656  ! dummy
657  class(PrtPrpType) :: this
658  ! local
659  integer(I4B) :: indx
660  call this%obs%StoreObsType('prp', .true., indx)
661  this%obs%obsData(indx)%ProcessIdPtr => defaultobsidprocessor
662 
663  ! Store obs type and assign procedure pointer
664  ! for to-mvr observation type.
665  call this%obs%StoreObsType('to-mvr', .true., indx)
666  this%obs%obsData(indx)%ProcessIdPtr => defaultobsidprocessor
Here is the call graph for this function:

◆ prp_dimensions()

subroutine prtprpmodule::prp_dimensions ( class(prtprptype), intent(inout)  this)

Definition at line 845 of file prt-prp.f90.

846  ! -- modules
849  ! -- dummy variables
850  class(PrtPrpType), intent(inout) :: this
851  ! -- local variables
852  type(PrtPrpParamFoundType) :: found
853 
854  call mem_set_value(this%nreleasepoints, 'NRELEASEPTS', this%input_mempath, &
855  found%nreleasepts)
856  call mem_set_value(this%nreleasetimes, 'NRELEASETIMES', this%input_mempath, &
857  found%nreleasetimes)
858 
859  write (this%iout, '(1x,a)') 'PROCESSING PARTICLE INPUT DIMENSIONS'
860  write (this%iout, '(4x,a,i0)') 'NRELEASEPTS = ', this%nreleasepoints
861  write (this%iout, '(4x,a,i0)') 'NRELEASETIMES = ', this%nreleasetimes
862  write (this%iout, '(1x,a)') 'END OF PARTICLE INPUT DIMENSIONS'
863 
864  ! set maxbound and nbound to nreleasepts
865  this%maxbound = this%nreleasepoints
866  this%nbound = this%nreleasepoints
867 
868  ! allocate arrays for prp package
869  call this%prp_allocate_arrays()
870 
871  ! read packagedata and releasetimes blocks
872  call this%prp_packagedata()
873  call this%prp_releasetimes()
874  call this%prp_load_releasetimefrequency()

◆ prp_load_releasetimefrequency()

subroutine prtprpmodule::prp_load_releasetimefrequency ( class(prtprptype), intent(inout)  this)

Definition at line 1053 of file prt-prp.f90.

1054  ! modules
1055  use tdismodule, only: totalsimtime
1056  ! dummy
1057  class(PrtPrpType), intent(inout) :: this
1058  ! local
1059  real(DP), allocatable :: times(:)
1060 
1061  ! check if a release time frequency is configured
1062  if (this%rtfreq <= dzero) return
1063 
1064  ! create array of regularly-spaced release times
1065  times = arange( &
1066  start=dzero, &
1067  stop=totalsimtime, &
1068  step=this%rtfreq)
1069 
1070  ! register times with release schedule
1071  call this%schedule%time_select%extend(times)
1072 
1073  ! make sure times strictly increase
1074  if (.not. this%schedule%time_select%increasing()) then
1075  errmsg = "Release times must strictly increase"
1076  call store_error(errmsg)
1077  call store_error_filename(this%input_fname)
1078  end if
1079 
1080  ! deallocate
1081  deallocate (times)
1082 
Here is the call graph for this function:

◆ prp_log_options()

subroutine prtprpmodule::prp_log_options ( class(prtprptype), intent(inout)  this,
type(prtprpparamfoundtype), intent(in)  found,
character(len=*), intent(in)  trackfile,
character(len=*), intent(in)  trackcsvfile 
)

Definition at line 809 of file prt-prp.f90.

810  ! -- modules
812  ! -- dummy variables
813  class(PrtPrpType), intent(inout) :: this
814  type(PrtPrpParamFoundType), intent(in) :: found
815  character(len=*), intent(in) :: trackfile
816  character(len=*), intent(in) :: trackcsvfile
817  ! -- local variables
818  ! formats
819  character(len=*), parameter :: fmttrkbin = &
820  "(4x, 'PARTICLE TRACKS WILL BE SAVED TO BINARY FILE: ', a, /4x, &
821  &'OPENED ON UNIT: ', I0)"
822  character(len=*), parameter :: fmttrkcsv = &
823  "(4x, 'PARTICLE TRACKS WILL BE SAVED TO CSV FILE: ', a, /4x, &
824  &'OPENED ON UNIT: ', I0)"
825 
826  write (this%iout, '(1x,a)') 'PROCESSING PARTICLE INPUT DIMENSIONS'
827 
828  if (found%frctrn) then
829  write (this%iout, '(4x,a)') &
830  'IF DISV, TRACKING WILL USE THE TERNARY METHOD REGARDLESS OF CELL TYPE'
831  end if
832 
833  if (found%trackfile) then
834  write (this%iout, fmttrkbin) trim(adjustl(trackfile)), this%itrkout
835  end if
836 
837  if (found%trackcsvfile) then
838  write (this%iout, fmttrkcsv) trim(adjustl(trackcsvfile)), this%itrkcsv
839  end if
840 
841  write (this%iout, '(1x,a)') 'END OF PARTICLE INPUT DIMENSIONS'

◆ prp_obs_supported()

logical function prtprpmodule::prp_obs_supported ( class(prtprptype this)
private

Definition at line 649 of file prt-prp.f90.

650  class(PrtPrpType) :: this
651  prp_obs_supported = .true.

◆ prp_options()

subroutine prtprpmodule::prp_options ( class(prtprptype), intent(inout)  this)
private

Definition at line 670 of file prt-prp.f90.

671  ! -- modules
674  use openspecmodule, only: access, form
677  ! -- dummy variables
678  class(PrtPrpType), intent(inout) :: this
679  ! -- local variables
680  character(len=LENVARNAME), dimension(3) :: drytrack_method = &
681  &[character(len=LENVARNAME) :: 'DROP', 'STOP', 'STAY']
682  character(len=lenvarname), dimension(2) :: coorcheck_method = &
683  &[character(len=LENVARNAME) :: 'NONE', 'EAGER']
684  character(len=LINELENGTH) :: trackfile, trackcsvfile, fname
685  type(PrtPrpParamFoundType) :: found
686 
687  ! -- source base class options
688  call this%BndExtType%source_options()
689 
690  ! -- update defaults from input context
691  call mem_set_value(this%stoptime, 'STOPTIME', this%input_mempath, &
692  found%stoptime)
693  call mem_set_value(this%stoptraveltime, 'STOPTRAVELTIME', &
694  this%input_mempath, found%stoptraveltime)
695  call mem_set_value(this%istopweaksink, 'ISTOPWEAKSINK', this%input_mempath, &
696  found%istopweaksink)
697  call mem_set_value(this%istopzone, 'ISTOPZONE', this%input_mempath, &
698  found%istopzone)
699  call mem_set_value(this%drape, 'DRAPE', this%input_mempath, &
700  found%drape)
701  call mem_set_value(this%idrymeth, 'IDRYMETH', this%input_mempath, &
702  drytrack_method, found%idrymeth)
703  call mem_set_value(trackfile, 'TRACKFILE', this%input_mempath, &
704  found%trackfile)
705  call mem_set_value(trackcsvfile, 'TRACKCSVFILE', this%input_mempath, &
706  found%trackcsvfile)
707  call mem_set_value(this%localz, 'LOCALZ', this%input_mempath, &
708  found%localz)
709  call mem_set_value(this%extend, 'EXTEND', this%input_mempath, &
710  found%extend)
711  call mem_set_value(this%extol, 'EXTOL', this%input_mempath, &
712  found%extol)
713  call mem_set_value(this%rttol, 'RTTOL', this%input_mempath, &
714  found%rttol)
715  call mem_set_value(this%rtfreq, 'RTFREQ', this%input_mempath, &
716  found%rtfreq)
717  call mem_set_value(this%frctrn, 'FRCTRN', this%input_mempath, &
718  found%frctrn)
719  call mem_set_value(this%iexmeth, 'IEXMETH', this%input_mempath, &
720  found%iexmeth)
721  call mem_set_value(this%ichkmeth, 'ICHKMETH', this%input_mempath, &
722  coorcheck_method, found%ichkmeth)
723  call mem_set_value(this%icycwin, 'ICYCWIN', this%input_mempath, found%icycwin)
724 
725  ! update internal state and validate input
726  if (found%idrymeth) then
727  if (this%idrymeth == 0) then
728  write (errmsg, '(a)') 'Unsupported dry tracking method. &
729  &DRY_TRACKING_METHOD must be "DROP", "STOP", or "STAY"'
730  call store_error(errmsg)
731  else
732  ! adjust for method zero indexing
733  this%idrymeth = this%idrymeth - 1
734  end if
735  end if
736 
737  if (found%extol) then
738  if (this%extol <= dzero) &
739  call store_error('EXIT_SOLVE_TOLERANCE MUST BE POSITIVE')
740  end if
741 
742  if (found%rttol) then
743  if (this%rttol <= dzero) &
744  call store_error('RELEASE_TIME_TOLERANCE MUST BE POSITIVE')
745  end if
746 
747  if (found%rtfreq) then
748  if (this%rtfreq <= dzero) &
749  call store_error('RELEASE_TIME_FREQUENCY MUST BE POSITIVE')
750  end if
751 
752  if (found%iexmeth) then
753  if (.not. (this%iexmeth /= 1 .or. this%iexmeth /= 2)) &
754  call store_error('DEV_EXIT_SOLVE_METHOD MUST BE &
755  &1 (BRENT) OR 2 (CHANDRUPATLA)')
756  end if
757 
758  if (found%ichkmeth) then
759  if (this%ichkmeth == 0) then
760  write (errmsg, '(a)') 'Unsupported coordinate check method. &
761  &COORDINATE_CHECK_METHOD must be "NONE" or "EAGER"'
762  call store_error(errmsg)
763  else
764  ! adjust for method zero based indexing
765  this%ichkmeth = this%ichkmeth - 1
766  end if
767  end if
768 
769  if (found%icycwin) then
770  if (this%icycwin < 0) &
771  call store_error('CYCLE_DETECTION_WINDOW MUST BE NON-NEGATIVE')
772  end if
773 
774  ! fileout options
775  if (found%trackfile) then
776  this%itrkout = getunit()
777  call openfile(this%itrkout, this%iout, trackfile, 'DATA(BINARY)', &
778  form, access, filstat_opt='REPLACE', &
779  mode_opt=mnormal)
780  ! open and write ascii header spec file
781  this%itrkhdr = getunit()
782  fname = trim(trackfile)//'.hdr'
783  call openfile(this%itrkhdr, this%iout, fname, 'CSV', &
784  filstat_opt='REPLACE', mode_opt=mnormal)
785  write (this%itrkhdr, '(a,/,a)') trackheader, trackdtypes
786  end if
787 
788  if (found%trackcsvfile) then
789  this%itrkcsv = getunit()
790  call openfile(this%itrkcsv, this%iout, trackcsvfile, 'CSV', &
791  filstat_opt='REPLACE')
792  write (this%itrkcsv, '(a)') trackheader
793  end if
794 
795  ! terminate if any errors were detected
796  if (count_errors() > 0) then
797  call store_error_filename(this%input_fname)
798  end if
799 
800  ! log found options
801  call this%prp_log_options(found, trackfile, trackcsvfile)
802 
803  ! Create release schedule now that we know
804  ! the coincident release time tolerance
805  this%schedule => create_release_schedule(tolerance=this%rttol)
This module contains simulation constants.
Definition: Constants.f90:9
@ mnormal
normal output mode
Definition: Constants.f90:206
integer(i4b), parameter lenvarname
maximum length of a variable name
Definition: Constants.f90:17
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
integer(i4b) function, public getunit()
Get a free unit number.
subroutine, public openfile(iu, iout, fname, ftype, fmtarg_opt, accarg_opt, filstat_opt, mode_opt)
Open a file.
Definition: InputOutput.f90:30
character(len=20) access
Definition: OpenSpec.f90:7
character(len=20) form
Definition: OpenSpec.f90:7
Here is the call graph for this function:

◆ prp_packagedata()

subroutine prtprpmodule::prp_packagedata ( class(prtprptype), intent(inout)  this)

Definition at line 878 of file prt-prp.f90.

880  use geomutilmodule, only: get_node
882  ! dummy
883  class(PrtPrpType), intent(inout) :: this
884  ! local
885  integer(I4B), dimension(:), pointer, contiguous :: irptno
886  integer(I4B), dimension(:, :), pointer, contiguous :: cellids
887  real(DP), dimension(:), pointer, contiguous :: xrpts, yrpts, zrpts
888  type(CharacterStringType), dimension(:), pointer, &
889  contiguous :: boundnames
890  character(len=LENBOUNDNAME) :: bndName, bndNameTemp
891  character(len=9) :: cno
892  integer(I4B), dimension(:), allocatable :: nboundchk
893  integer(I4B), dimension(:), pointer :: cellid
894  integer(I4B) :: n, noder, nodeu, rptno
895 
896  ! set input context pointers
897  call mem_setptr(irptno, 'IRPTNO', this%input_mempath)
898  call mem_setptr(cellids, 'CELLID', this%input_mempath)
899  call mem_setptr(xrpts, 'XRPT', this%input_mempath)
900  call mem_setptr(yrpts, 'YRPT', this%input_mempath)
901  call mem_setptr(zrpts, 'ZRPT', this%input_mempath)
902  call mem_setptr(boundnames, 'BOUNDNAME', this%input_mempath)
903 
904  ! allocate and initialize temporary variables
905  allocate (nboundchk(this%nreleasepoints))
906  do n = 1, this%nreleasepoints
907  nboundchk(n) = 0
908  end do
909 
910  write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%packName)) &
911  //' PACKAGEDATA'
912 
913  do n = 1, size(irptno)
914 
915  rptno = irptno(n)
916 
917  if (rptno < 1 .or. rptno > this%nreleasepoints) then
918  write (errmsg, '(a,i0,a,i0,a)') &
919  'Expected ', this%nreleasepoints, ' release points. &
920  &Points must be numbered from 1 to ', this%nreleasepoints, '.'
921  call store_error(errmsg)
922  cycle
923  end if
924 
925  ! increment nboundchk
926  nboundchk(rptno) = nboundchk(rptno) + 1
927 
928  ! set cellid
929  cellid => cellids(:, n)
930 
931  ! set node user
932  if (this%dis%ndim == 1) then
933  nodeu = cellid(1)
934  elseif (this%dis%ndim == 2) then
935  nodeu = get_node(cellid(1), 1, cellid(2), &
936  this%dis%mshape(1), 1, &
937  this%dis%mshape(2))
938  else
939  nodeu = get_node(cellid(1), cellid(2), cellid(3), &
940  this%dis%mshape(1), &
941  this%dis%mshape(2), &
942  this%dis%mshape(3))
943  end if
944 
945  ! set noder
946  noder = this%dis%get_nodenumber(nodeu, 1)
947  if (noder <= 0) then
948  cycle
949  else
950  this%rptnode(rptno) = noder
951  end if
952 
953  if (this%localz .and. (zrpts(n) < 0 .or. zrpts(n) > 1)) then
954  call store_error('Local z coordinate must fall in the interval [0, 1]')
955  cycle
956  end if
957 
958  ! set coordinates
959  this%rptx(rptno) = xrpts(n)
960  this%rpty(rptno) = yrpts(n)
961  this%rptz(rptno) = zrpts(n)
962 
963  ! set default boundname
964  write (cno, '(i9.9)') rptno
965  bndname = 'PRP'//cno
966 
967  ! read boundnames from file, if provided
968  if (this%inamedbound /= 0) then
969  bndnametemp = boundnames(n)
970  if (bndnametemp /= '') bndname = bndnametemp
971  else
972  bndname = ''
973  end if
974 
975  ! set boundname
976  this%rptname(rptno) = bndname
977  end do
978 
979  write (this%iout, '(1x,a)') &
980  'END OF '//trim(adjustl(this%packName))//' PACKAGEDATA'
981 
982  ! check for duplicate or missing particle release points
983  do n = 1, this%nreleasepoints
984  if (nboundchk(n) == 0) then
985  write (errmsg, '(a,a,1x,i0,a)') 'No data specified for particle ', &
986  'release point', n, '.'
987  call store_error(errmsg)
988  else if (nboundchk(n) > 1) then
989  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
990  'Data for particle release point', n, 'specified', nboundchk(n), &
991  'times.'
992  call store_error(errmsg)
993  end if
994  end do
995 
996  ! terminate if any errors were detected
997  if (count_errors() > 0) then
998  call store_error_filename(this%input_fname)
999  end if
1000 
1001  ! cleanup
1002  deallocate (nboundchk)
integer(i4b) function, public get_node(ilay, irow, icol, nlay, nrow, ncol)
Get node number, given layer, row, and column indices for a structured grid. If any argument is inval...
Definition: GeomUtil.f90:83
This class is used to store a single deferred-length character string. It was designed to work in an ...
Definition: CharString.f90:23
Here is the call graph for this function:

◆ prp_releasetimes()

subroutine prtprpmodule::prp_releasetimes ( class(prtprptype), intent(inout)  this)

Definition at line 1006 of file prt-prp.f90.

1008  ! dummy
1009  class(PrtPrpType), intent(inout) :: this
1010  ! local
1011  real(DP), dimension(:), pointer, contiguous :: time
1012  integer(I4B) :: n, isize
1013  real(DP), allocatable :: times(:)
1014 
1015  if (this%nreleasetimes <= 0) return
1016 
1017  ! allocate times array
1018  allocate (times(this%nreleasetimes))
1019 
1020  ! check if input array was read
1021  call get_isize('TIME', this%input_mempath, isize)
1022 
1023  if (isize <= 0) then
1024  errmsg = "RELEASTIMES block expected when &
1025  &NRELEASETIMES dimension is non-zero."
1026  call store_error(errmsg)
1027  call store_error_filename(this%input_fname)
1028  end if
1029 
1030  ! set input context pointer
1031  call mem_setptr(time, 'TIME', this%input_mempath)
1032 
1033  ! set input data
1034  do n = 1, size(time)
1035  times(n) = time(n)
1036  end do
1037 
1038  ! register times with the release schedule
1039  call this%schedule%time_select%extend(times)
1040 
1041  ! make sure times strictly increase
1042  if (.not. this%schedule%time_select%increasing()) then
1043  errmsg = "RELEASTIMES block entries must strictly increase."
1044  call store_error(errmsg)
1045  call store_error_filename(this%input_fname)
1046  end if
1047 
1048  ! deallocate
1049  deallocate (times)
subroutine, public get_isize(name, mem_path, isize)
@ brief Get the number of elements for this variable
Here is the call graph for this function:

◆ prp_rp()

subroutine prtprpmodule::prp_rp ( class(prtprptype), intent(inout)  this)

Definition at line 555 of file prt-prp.f90.

556  ! modules
557  use tdismodule, only: kper, nper
560  ! dummy variables
561  class(PrtPrpType), intent(inout) :: this
562  ! local variables
563  type(CharacterStringType), dimension(:), contiguous, &
564  pointer :: settings
565  integer(I4B), pointer :: iper, ionper, nlist
566  character(len=LINELENGTH), allocatable :: lines(:)
567  integer(I4B) :: n
568 
569  ! set pointer to last and next period loaded
570  call mem_setptr(iper, 'IPER', this%input_mempath)
571  call mem_setptr(ionper, 'IONPER', this%input_mempath)
572 
573  if (kper == 1 .and. &
574  (iper == 0) .and. &
575  (ionper > nper) .and. &
576  size(this%schedule%time_select%times) == 0) then
577  ! If the user hasn't provided any release settings (neither
578  ! explicit release times, release time frequency, or period
579  ! block release settings), default to a single release at the
580  ! start of the first period's first time step.
581  allocate (lines(1))
582  lines(1) = "FIRST"
583  call this%schedule%advance(lines=lines)
584  deallocate (lines)
585  return
586  else if (iper /= kper) then
587  return
588  end if
589 
590  ! set input context pointers
591  call mem_setptr(nlist, 'NBOUND', this%input_mempath)
592  call mem_setptr(settings, 'SETTING', this%input_mempath)
593 
594  ! allocate and set input
595  allocate (lines(nlist))
596  do n = 1, nlist
597  lines(n) = settings(n)
598  end do
599 
600  ! update schedule
601  if (size(lines) > 0) &
602  call this%schedule%advance(lines=lines)
603 
604  ! cleanup
605  deallocate (lines)
integer(i4b), pointer, public kper
current stress period number
Definition: tdis.f90:23
integer(i4b), pointer, public nper
number of stress period
Definition: tdis.f90:21

◆ prp_set_pointers()

subroutine prtprpmodule::prp_set_pointers ( class(prtprptype this,
integer(i4b), dimension(:), pointer, contiguous  ibound,
integer(i4b), dimension(:), pointer, contiguous  izone 
)
private

Definition at line 196 of file prt-prp.f90.

197  class(PrtPrpType) :: this
198  integer(I4B), dimension(:), pointer, contiguous :: ibound
199  integer(I4B), dimension(:), pointer, contiguous :: izone
200 
201  this%ibound => ibound
202  this%rptzone => izone

◆ release()

subroutine prtprpmodule::release ( class(prtprptype), intent(inout)  this,
integer(i4b), intent(in)  ip,
real(dp), intent(in)  trelease 
)
private

Releasing a particle entails validating the particle's coordinates and settings, transforming its coordinates if needed, initializing the particle's initial tracking time to the given release time, storing the particle in the particle store (from which the PRT model will later retrieve it, apply the tracking method, and check it in again), and accumulating the particle's mass (the total mass released from each release point is calculated for budget reporting).

Parameters
[in,out]thisthis instance
[in]ipparticle index
[in]treleaserelease time

Definition at line 438 of file prt-prp.f90.

439  ! dummy
440  class(PrtPrpType), intent(inout) :: this !< this instance
441  integer(I4B), intent(in) :: ip !< particle index
442  real(DP), intent(in) :: trelease !< release time
443  ! local
444  integer(I4B) :: np
445  type(ParticleType), pointer :: particle
446 
447  call this%initialize_particle(particle, ip, trelease)
448  np = this%nparticles + 1
449  this%nparticles = np
450  call this%particles%put(particle, np)
451  deallocate (particle)
452  this%rptm(ip) = this%rptm(ip) + done ! TODO configurable mass
453 

◆ validate_release_point()

subroutine prtprpmodule::validate_release_point ( class(prtprptype), intent(inout)  this,
integer(i4b), intent(in)  ic,
real(dp), intent(in)  x,
real(dp), intent(in)  y,
real(dp), intent(in)  z 
)
private

Terminate with an error if the release point lies outside the given cell, or if the point is above or below the grid top or bottom, respectively.

Parameters
[in,out]thisthis instance
[in]iccell index
[in]zrelease point

Definition at line 395 of file prt-prp.f90.

396  class(PrtPrpType), intent(inout) :: this !< this instance
397  integer(I4B), intent(in) :: ic !< cell index
398  real(DP), intent(in) :: x, y, z !< release point
399  ! local
400  real(DP), allocatable :: polyverts(:, :)
401 
402  call this%fmi%dis%get_polyverts(ic, polyverts)
403  if (.not. point_in_polygon(x, y, polyverts)) then
404  write (errmsg, '(a,g0,a,g0,a,i0)') &
405  'Error: release point (x=', x, ', y=', y, ') is not in cell ', &
406  this%dis%get_nodeuser(ic)
407  call store_error(errmsg, terminate=.false.)
408  call store_error_filename(this%input_fname)
409  end if
410  if (z > maxval(this%dis%top)) then
411  write (errmsg, '(a,g0,a,g0,a,i0)') &
412  'Error: release point (z=', z, ') is above grid top ', &
413  maxval(this%dis%top)
414  call store_error(errmsg, terminate=.false.)
415  call store_error_filename(this%input_fname)
416  else if (z < minval(this%dis%bot)) then
417  write (errmsg, '(a,g0,a,g0,a,i0)') &
418  'Error: release point (z=', z, ') is below grid bottom ', &
419  minval(this%dis%bot)
420  call store_error(errmsg, terminate=.false.)
421  call store_error_filename(this%input_fname)
422  end if
423  deallocate (polyverts)
Here is the call graph for this function:

Variable Documentation

◆ ftype

character(len=lenftype) prtprpmodule::ftype = 'PRP'
private

Definition at line 34 of file prt-prp.f90.

34  character(len=LENFTYPE) :: ftype = 'PRP'

◆ text

character(len=16) prtprpmodule::text = ' PRP'
private

Definition at line 35 of file prt-prp.f90.

35  character(len=16) :: text = ' PRP'