MODFLOW 6  version 6.8.0.dev0
USGS Modular Hydrologic Model
prtprpmodule Module Reference

Data Types

type  prtprptype
 Particle release point (PRP) package. More...
 
type  exgprtprptype
 Exchange PRP package. A variant of the normal PRP package that doesn't read from input files but instead receives particle transfers from coupled models while preserving the pattern where PRP packages own particles. Call it "Particle Registry Package"? More...
 

Functions/Subroutines

subroutine, public prp_create (packobj, id, ibcnum, inunit, iout, namemodel, pakname, fmi, input_mempath)
 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 exg_prp_allocate_scalars (this)
 Allocate scalars for exchange PRP. More...
 
subroutine exg_prp_allocate_arrays (this, nodelist, auxvar)
 Allocate arrays for exchange PRP. More...
 
subroutine exg_prp_ar (this)
 @ brief No-op AR method override for exchange PRP. More...
 
subroutine prp_ad (this)
 Advance a time step and release particles if scheduled. More...
 
subroutine exg_prp_ad (this)
 No-op AD method override for exchange PRP. More...
 
subroutine exg_prp_cq_simrate (this, hnew, flowja, imover)
 No-op flow calculation for exchange PRP. More...
 
subroutine exg_prp_bd (this, model_budget)
 No-op budget method for exchange PRP. Likewise about the STORAGE term accounting. More...
 
subroutine exg_prp_ot_model_flows (this, icbcfl, ibudfl, icbcun, imap)
 No-op flow output method for exchange PRP. No contribution to budget, no need to write output. 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 exg_prp_rp (this)
 @ brief No-op RP method override for exchange PRP. 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 exg_prp_options (this)
 @ brief No-op options method override for exchange PRP. Just creates an empty release schedule. 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 exg_prp_dimensions (this)
 @ brief Dimensions method override for exchange PRP. Just set all dimensions to zero and allocate arrays. 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'
 
real(dp), parameter default_exit_solve_tolerance = DEM5
 

Function/Subroutine Documentation

◆ define_listlabel()

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

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

811  class(PrtPrpType), intent(inout) :: this
812  ! not implemented, not used

◆ exg_prp_ad()

subroutine prtprpmodule::exg_prp_ad ( class(exgprtprptype this)

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

492  class(ExgPrtPrpType) :: this

◆ exg_prp_allocate_arrays()

subroutine prtprpmodule::exg_prp_allocate_arrays ( class(exgprtprptype this,
integer(i4b), dimension(:), optional, pointer, contiguous  nodelist,
real(dp), dimension(:, :), optional, pointer, contiguous  auxvar 
)

BndExtType expects certain array variables to exist in the input context (CELLID, NODEULIST, BOUNDNAME, AUXVAR). This method manually creates zero-sized arrays before calling the parent's allocate_arrays.

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

394  class(ExgPrtPrpType) :: this
395  integer(I4B), dimension(:), pointer, contiguous, optional :: nodelist
396  real(DP), dimension(:, :), pointer, contiguous, optional :: auxvar
397  ! local
398  integer(I4B), dimension(:, :), pointer, contiguous :: cellid
399  integer(I4B), dimension(:), pointer, contiguous :: nodeulist
400  type(CharacterStringType), dimension(:), pointer, contiguous :: boundname
401  real(DP), dimension(:, :), pointer, contiguous :: auxvar_input
402 
403  call mem_allocate(cellid, this%dis%ndim, 0, 'CELLID', this%input_mempath)
404  call mem_allocate(nodeulist, 0, 'NODEULIST', this%input_mempath)
405  call mem_allocate(boundname, lenboundname, 0, 'BOUNDNAME', &
406  this%input_mempath)
407  call mem_allocate(auxvar_input, 0, 0, 'AUXVAR', this%input_mempath)
408 
409  call this%PrtPrpType%prp_allocate_arrays(nodelist, auxvar)
This class is used to store a single deferred-length character string. It was designed to work in an ...
Definition: CharString.f90:23

◆ exg_prp_allocate_scalars()

subroutine prtprpmodule::exg_prp_allocate_scalars ( class(exgprtprptype this)
private

The exchange PRP is a headless package (no input file) but BndExtType expects certain variables to exist in the input context (IPER, IONPER) so we need to manually create them before calling the parent procedure.

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

370  class(ExgPrtPrpType) :: this
371  integer(I4B), pointer :: iper, ionper
372 
373  this%input_mempath = trim(this%memoryPath)//'-INPUT'
374 
375  call mem_allocate(iper, 'IPER', this%input_mempath)
376  call mem_allocate(ionper, 'IONPER', this%input_mempath)
377 
378  ! set iper = 0. this forces BndExtType%bnd_rp to
379  ! return early, since iper will never match kper.
380  iper = 0
381  ionper = 0
382 
383  call this%PrtPrpType%prp_allocate_scalars()

◆ exg_prp_ar()

subroutine prtprpmodule::exg_prp_ar ( class(exgprtprptype), intent(inout)  this)

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

414  class(ExgPrtPrpType), intent(inout) :: this

◆ exg_prp_bd()

subroutine prtprpmodule::exg_prp_bd ( class(exgprtprptype this,
type(budgettype), intent(inout)  model_budget 
)
private

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

509  use budgetmodule, only: budgettype
510  class(ExgPrtPrpType) :: this
511  type(BudgetType), intent(inout) :: model_budget
This module contains the BudgetModule.
Definition: Budget.f90:20
Derived type for the Budget object.
Definition: Budget.f90:39

◆ exg_prp_cq_simrate()

subroutine prtprpmodule::exg_prp_cq_simrate ( class(exgprtprptype this,
real(dp), dimension(:), intent(in)  hnew,
real(dp), dimension(:), intent(inout)  flowja,
integer(i4b), intent(in)  imover 
)
private

Exchange PRP particles are transferred from other models and already active, so their mass is already accounted for in the STORAGE term.

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

500  class(ExgPrtPrpType) :: this
501  real(DP), dimension(:), intent(in) :: hnew
502  real(DP), dimension(:), intent(inout) :: flowja
503  integer(I4B), intent(in) :: imover

◆ exg_prp_dimensions()

subroutine prtprpmodule::exg_prp_dimensions ( class(exgprtprptype), intent(inout)  this)

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

1063  class(ExgPrtPrpType), intent(inout) :: this
1064 
1065  this%nreleasepoints = 0
1066  this%nreleasetimes = 0
1067  this%maxbound = 0
1068  this%nbound = 0
1069 
1070  call this%prp_allocate_arrays()

◆ exg_prp_options()

subroutine prtprpmodule::exg_prp_options ( class(exgprtprptype), intent(inout)  this)

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

990  class(ExgPrtPrpType), intent(inout) :: this
991  this%schedule => create_release_schedule(tolerance=this%rttol)
Here is the call graph for this function:

◆ exg_prp_ot_model_flows()

subroutine prtprpmodule::exg_prp_ot_model_flows ( class(exgprtprptype this,
integer(i4b), intent(in)  icbcfl,
integer(i4b), intent(in)  ibudfl,
integer(i4b), intent(in)  icbcun,
integer(i4b), dimension(:), intent(in), optional  imap 
)

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

517  class(ExgPrtPrpType) :: this
518  integer(I4B), intent(in) :: icbcfl
519  integer(I4B), intent(in) :: ibudfl
520  integer(I4B), intent(in) :: icbcun
521  integer(I4B), dimension(:), optional, intent(in) :: imap

◆ exg_prp_rp()

subroutine prtprpmodule::exg_prp_rp ( class(exgprtprptype), intent(inout)  this)

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

772  class(ExgPrtPrpType), intent(inout) :: this

◆ 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 601 of file prt-prp.f90.

603  class(PrtPrpType), intent(inout) :: this !< this instance
604  type(ParticleType), pointer, intent(inout) :: particle !< the particle
605  integer(I4B), intent(in) :: ip !< particle index
606  real(DP), intent(in) :: trelease !< release time
607  ! local
608  integer(I4B) :: irow, icol, ilay, icpl
609  integer(I4B) :: ic, icu, ic_old
610  real(DP) :: x, y, z
611  real(DP) :: top, bot, hds
612  ! formats
613  character(len=*), parameter :: fmticterr = &
614  "('Error in ',a,': Flow model interface does not contain ICELLTYPE. &
615  &ICELLTYPE is required for PRT to distinguish convertible cells &
616  &from confined cells if LOCAL_Z release coordinates are provided. &
617  &Make sure a GWFGRID entry is configured in the PRT FMI package.')"
618 
619  ic = this%rptnode(ip)
620  icu = this%dis%get_nodeuser(ic)
621 
622  call create_particle(particle)
623 
624  if (size(this%boundname) /= 0) then
625  particle%name = this%boundname(ip)
626  else
627  particle%name = ''
628  end if
629 
630  particle%irpt = ip
631  particle%istopweaksink = this%istopweaksink
632  particle%istopzone = this%istopzone
633  particle%idrymeth = this%idrymeth
634  particle%icu = icu
635 
636  select type (dis => this%dis)
637  type is (distype)
638  call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, irow, icol, ilay)
639  type is (disvtype)
640  call get_jk(icu, dis%ncpl, dis%nlay, icpl, ilay)
641  end select
642  particle%ilay = ilay
643  particle%izone = this%rptzone(ic)
644  particle%istatus = 0 ! status 0 until tracking starts
645  ! If the cell is inactive, either drape the particle
646  ! to the top-most active cell beneath it if drape is
647  ! enabled, or else terminate permanently unreleased.
648  if (this%ibound(ic) == 0) then
649  ic_old = ic
650  if (this%drape) then
651  call this%dis%highest_active(ic, this%ibound)
652  if (ic == ic_old .or. this%ibound(ic) == 0) then
653  ! negative unreleased status signals to the
654  ! tracking method that we haven't yet saved
655  ! a termination record, it needs to do so.
656  particle%istatus = -1 * term_unreleased
657  end if
658  else
659  particle%istatus = -1 * term_unreleased
660  end if
661  end if
662 
663  ! load coordinates
664  x = this%rptx(ip)
665  y = this%rpty(ip)
666  if (this%localz) then
667  ! make sure FMI has cell type array. we need
668  ! it to distinguish convertible and confined
669  ! cells if release z coordinates are local
670  if (this%fmi%igwfceltyp /= 1) then
671  write (errmsg, fmticterr) trim(this%text)
672  call store_error(errmsg, terminate=.true.)
673  end if
674 
675  ! calculate model z coord from local z coord.
676  ! if cell is confined (icelltype == 0) use the
677  ! actual cell height (geometric top - bottom).
678  ! otherwise use head as cell top, clamping to
679  ! the cell bottom if head is below the bottom
680  ! and to geometric cell top if head is above top
681  top = this%fmi%dis%top(ic)
682  bot = this%fmi%dis%bot(ic)
683  if (this%fmi%gwfceltyp(icu) /= 0) then
684  hds = this%fmi%gwfhead(ic)
685  top = min(top, hds)
686  top = max(top, bot)
687  end if
688  z = bot + this%rptz(ip) * (top - bot)
689  else
690  z = this%rptz(ip)
691  end if
692 
693  if (this%ichkmeth > 0) &
694  call this%validate_release_point(ic, x, y, z)
695 
696  particle%x = x
697  particle%y = y
698  particle%z = z
699  particle%trelease = trelease
700 
701  ! Set stop time to earlier of STOPTIME and STOPTRAVELTIME
702  if (this%stoptraveltime == huge(1d0)) then
703  particle%tstop = this%stoptime
704  else
705  particle%tstop = particle%trelease + this%stoptraveltime
706  if (this%stoptime < particle%tstop) particle%tstop = this%stoptime
707  end if
708 
709  particle%ttrack = particle%trelease
710  particle%itrdomain(level_model) = 0
711  particle%iboundary(level_model) = 0
712  particle%itrdomain(level_feature) = ic
713  particle%iboundary(level_feature) = 0
714  particle%itrdomain(level_subfeature) = 0
715  particle%iboundary(level_subfeature) = 0
716  particle%frctrn = this%frctrn
717  particle%iexmeth = this%iexmeth
718  particle%extend = this%extend
719  particle%icycwin = this%icycwin
720  particle%extol = this%extol
@ term_unreleased
terminated permanently unreleased
Definition: Particle.f90:36
Here is the call graph for this function:

◆ log_release()

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

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

526  class(PrtPrpType), intent(inout) :: this !< prp
527  if (this%iprpak > 0) then
528  write (this%iout, "(1x,/1x,a,1x,i0)") &
529  'PARTICLE RELEASE FOR PRP', this%ibcnum
530  call this%schedule%log(this%iout)
531  end if

◆ prp_ad()

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

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

419  use tdismodule, only: totalsimtime, kstp, kper
420  class(PrtPrpType) :: this
421  integer(I4B) :: ip, it
422  real(DP) :: t
423 
424  ! Notes
425  ! -----
426  ! Each release point can be thought of as
427  ! a gumball machine with infinite supply:
428  ! a point can release an arbitrary number
429  ! of particles, but only one at any time.
430  ! Coincident release times are merged to
431  ! a single time by the release scheduler.
432 
433  ! Reset mass accumulators for this time step.
434  do ip = 1, this%nreleasepoints
435  this%rptm(ip) = dzero
436  end do
437 
438  ! Advance the release schedule. At the start of each period (kstp==1),
439  ! apply period block configuration if available and not yet applied.
440  ! This handles both new configuration and fill-forward periods.
441  ! For subsequent time steps, just advance without arguments to
442  ! advance the time selection object to the current time step.
443  if (kstp == 1 .and. &
444  kper /= this%applied_kper .and. &
445  allocated(this%period_block_lines)) then
446  call this%schedule%advance(lines=this%period_block_lines)
447  this%applied_kper = kper
448  else
449  call this%schedule%advance()
450  end if
451 
452  ! Check if any releases will be made this time step.
453  if (.not. this%schedule%any()) return
454 
455  ! Log the schedule to the list file.
456  call this%log_release()
457 
458  ! Expand the particle store. We know from the
459  ! schedule how many particles will be released.
460  call this%particles%resize( &
461  this%particles%num_stored() + &
462  (this%nreleasepoints * this%schedule%count()), &
463  this%memoryPath)
464 
465  ! Release a particle from each point for
466  ! each release time in the current step.
467  do ip = 1, this%nreleasepoints
468  do it = 1, this%schedule%count()
469  t = this%schedule%times(it)
470  ! Skip the release time if it's before the simulation
471  ! starts, or if no `extend_tracking`, after it ends.
472  if (t < dzero) then
473  write (warnmsg, '(a,g0,a)') &
474  'Skipping negative release time (t=', t, ').'
475  call store_warning(warnmsg)
476  cycle
477  else if (t > totalsimtime .and. .not. this%extend) then
478  write (warnmsg, '(a,g0,a)') &
479  'Skipping release time falling after the end of the &
480  &simulation (t=', t, '). Enable EXTEND_TRACKING to &
481  &release particles after the simulation end time.'
482  call store_warning(warnmsg)
483  cycle
484  end if
485  call this%release(ip, t)
486  end do
487  end do
real(dp), pointer, public totalsimtime
time at end of simulation
Definition: tdis.f90:37
integer(i4b), pointer, public kstp
current time step number
Definition: tdis.f90:24
integer(i4b), pointer, public kper
current stress period number
Definition: tdis.f90:23
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 252 of file prt-prp.f90.

253  ! dummy
254  class(PrtPrpType) :: this
255  integer(I4B), dimension(:), pointer, contiguous, optional :: nodelist
256  real(DP), dimension(:, :), pointer, contiguous, optional :: auxvar
257  ! local
258  integer(I4B) :: nps
259 
260  call this%BndExtType%allocate_arrays()
261 
262  ! Allocate particle store, starting with the number
263  ! of release points (arrays resized if/when needed)
264  call create_particle_store( &
265  this%particles, &
266  this%nreleasepoints, &
267  this%memoryPath)
268 
269  ! Allocate arrays
270  call mem_allocate(this%rptx, this%nreleasepoints, 'RPTX', this%memoryPath)
271  call mem_allocate(this%rpty, this%nreleasepoints, 'RPTY', this%memoryPath)
272  call mem_allocate(this%rptz, this%nreleasepoints, 'RPTZ', this%memoryPath)
273  call mem_allocate(this%rptm, this%nreleasepoints, 'RPTMASS', &
274  this%memoryPath)
275  call mem_allocate(this%rptnode, this%nreleasepoints, 'RPTNODER', &
276  this%memoryPath)
277  call mem_allocate(this%rptname, lenboundname, this%nreleasepoints, &
278  'RPTNAME', this%memoryPath)
279 
280  ! Initialize arrays
281  do nps = 1, this%nreleasepoints
282  this%rptm(nps) = dzero
283  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 287 of file prt-prp.f90.

288  class(PrtPrpType) :: this
289 
290  ! Allocate parent's scalars
291  call this%BndExtType%allocate_scalars()
292 
293  ! Allocate scalars for this type
294  call mem_allocate(this%localz, 'LOCALZ', this%memoryPath)
295  call mem_allocate(this%extend, 'EXTEND', this%memoryPath)
296  call mem_allocate(this%offset, 'OFFSET', this%memoryPath)
297  call mem_allocate(this%stoptime, 'STOPTIME', this%memoryPath)
298  call mem_allocate(this%stoptraveltime, 'STOPTRAVELTIME', this%memoryPath)
299  call mem_allocate(this%istopweaksink, 'ISTOPWEAKSINK', this%memoryPath)
300  call mem_allocate(this%istopzone, 'ISTOPZONE', this%memoryPath)
301  call mem_allocate(this%drape, 'DRAPE', this%memoryPath)
302  call mem_allocate(this%idrymeth, 'IDRYMETH', this%memoryPath)
303  call mem_allocate(this%nreleasepoints, 'NRELEASEPOINTS', this%memoryPath)
304  call mem_allocate(this%nreleasetimes, 'NRELEASETIMES', this%memoryPath)
305  call mem_allocate(this%nparticles, 'NPARTICLES', this%memoryPath)
306  call mem_allocate(this%itrkout, 'ITRKOUT', this%memoryPath)
307  call mem_allocate(this%itrkhdr, 'ITRKHDR', this%memoryPath)
308  call mem_allocate(this%itrkcsv, 'ITRKCSV', this%memoryPath)
309  call mem_allocate(this%irlstls, 'IRLSTLS', this%memoryPath)
310  call mem_allocate(this%frctrn, 'FRCTRN', this%memoryPath)
311  call mem_allocate(this%iexmeth, 'IEXMETH', this%memoryPath)
312  call mem_allocate(this%ichkmeth, 'ICHKMETH', this%memoryPath)
313  call mem_allocate(this%icycwin, 'ICYCWIN', this%memoryPath)
314  call mem_allocate(this%extol, 'EXTOL', this%memoryPath)
315  call mem_allocate(this%rttol, 'RTTOL', this%memoryPath)
316  call mem_allocate(this%rtfreq, 'RTFREQ', this%memoryPath)
317 
318  ! Set values
319  this%localz = .false.
320  this%extend = .false.
321  this%offset = dzero
322  this%stoptime = huge(1d0)
323  this%stoptraveltime = huge(1d0)
324  this%istopweaksink = 0
325  this%istopzone = 0
326  this%drape = .false.
327  this%idrymeth = 0
328  this%nreleasepoints = 0
329  this%nreleasetimes = 0
330  this%nparticles = 0
331  this%itrkout = 0
332  this%itrkhdr = 0
333  this%itrkcsv = 0
334  this%irlstls = 0
335  this%frctrn = .false.
336  this%iexmeth = 0
337  this%ichkmeth = 1
338  this%icycwin = 0
339  this%extol = default_exit_solve_tolerance
340  this%rttol = dsame * dep9
341  this%rtfreq = dzero
342  this%applied_kper = 0
343 

◆ prp_ar()

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

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

348  class(PrtPrpType), intent(inout) :: this
349  integer(I4B) :: n
350 
351  call this%obs%obs_ar()
352 
353  if (this%inamedbound /= 0) then
354  do n = 1, this%nreleasepoints
355  this%boundname(n) = this%rptname(n)
356  end do
357  end if
358  do n = 1, this%nreleasepoints
359  this%nodelist(n) = this%rptnode(n)
360  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 
)
private
Parameters
[in,out]flowjaflow between package and model
[in]imoverflag indicating if the mover package is active

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

777  ! modules
778  use tdismodule, only: delt
779  ! dummy variables
780  class(PrtPrpType) :: this
781  real(DP), dimension(:), intent(in) :: hnew
782  real(DP), dimension(:), intent(inout) :: flowja !< flow between package and model
783  integer(I4B), intent(in) :: imover !< flag indicating if the mover package is active
784  ! local variables
785  integer(I4B) :: i
786  integer(I4B) :: node
787  integer(I4B) :: idiag
788  real(DP) :: rrate
789 
790  ! If no boundaries, skip flow calculations.
791  if (this%nbound <= 0) return
792 
793  ! Loop through each boundary calculating flow.
794  do i = 1, this%nbound
795  node = this%nodelist(i)
796  rrate = dzero
797  ! If cell is no-flow or constant-head, then ignore it.
798  if (node > 0) then
799  ! Calculate the flow rate into the cell.
800  idiag = this%dis%con%ia(node)
801  rrate = this%rptm(i) * (done / delt) ! reciprocal of tstp length
802  flowja(idiag) = flowja(idiag) + rrate
803  end if
804 
805  ! Save simulated value to simvals array.
806  this%simvals(i) = rrate
807  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,
type(prtfmitype), pointer  fmi,
character(len=*), intent(in), optional  input_mempath 
)

Creates either a standard PRP (reads from input file) or an exchange PRP (programmatically populated). The type is determined by whether input_mempath is provided: if present, standard; if absent, exchange.

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

127  ! dummy
128  class(BndType), pointer :: packobj
129  integer(I4B), intent(in) :: id
130  integer(I4B), intent(in) :: ibcnum
131  integer(I4B), intent(in) :: inunit
132  integer(I4B), intent(in) :: iout
133  character(len=*), intent(in) :: namemodel
134  character(len=*), intent(in) :: pakname
135  character(len=*), intent(in), optional :: input_mempath
136  type(PrtFmiType), pointer :: fmi
137  ! local
138  type(PrtPrpType), pointer :: prpobj
139  type(ExgPrtPrpType), pointer :: exgprpobj
140  ! formats
141  character(len=*), parameter :: fmtheader = &
142  "(1x, /1x, 'PRP PARTICLE RELEASE POINT PACKAGE', &
143  &' INPUT READ FROM MEMPATH: ', a, /)"
144  character(len=*), parameter :: fmtexgheader = &
145  "(1x, /1x, 'PRP-EXG EXCHANGE PARTICLE RELEASE POINT PACKAGE', &
146  &' (PROGRAMMATIC INPUT)', /)"
147 
148  if (present(input_mempath)) then
149  ! standard PRP
150  allocate (prpobj)
151  packobj => prpobj
152 
153  call packobj%set_names(ibcnum, namemodel, pakname, ftype, input_mempath)
154  prpobj%text = text
155 
156  call prpobj%prp_allocate_scalars()
157  call packobj%pack_initialize()
158 
159  packobj%inunit = inunit
160  packobj%iout = iout
161  packobj%id = id
162  packobj%ibcnum = ibcnum
163  packobj%ncolbnd = 4
164  packobj%iscloc = 1
165  prpobj%fmi => fmi
166 
167  if (inunit > 0) write (iout, fmtheader) input_mempath
168  else
169  ! exchange PRP
170  allocate (exgprpobj)
171  packobj => exgprpobj
172 
173  call packobj%set_names(ibcnum, namemodel, pakname, ftype)
174  exgprpobj%text = text
175 
176  call exgprpobj%prp_allocate_scalars()
177  call packobj%pack_initialize()
178 
179  packobj%inunit = inunit
180  packobj%iout = iout
181  packobj%id = id
182  packobj%ibcnum = ibcnum
183  packobj%ncolbnd = 4
184  packobj%iscloc = 1
185  exgprpobj%fmi => fmi
186 
187  if (iout > 0) write (iout, fmtexgheader)
188  end if
Here is the caller graph for this function:

◆ prp_da()

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

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

193  class(PrtPrpType) :: this
194 
195  ! Deallocate parent
196  call this%BndExtType%bnd_da()
197 
198  ! Deallocate scalars
199  call mem_deallocate(this%localz)
200  call mem_deallocate(this%extend)
201  call mem_deallocate(this%offset)
202  call mem_deallocate(this%stoptime)
203  call mem_deallocate(this%stoptraveltime)
204  call mem_deallocate(this%istopweaksink)
205  call mem_deallocate(this%istopzone)
206  call mem_deallocate(this%drape)
207  call mem_deallocate(this%idrymeth)
208  call mem_deallocate(this%nreleasepoints)
209  call mem_deallocate(this%nreleasetimes)
210  call mem_deallocate(this%nparticles)
211  call mem_deallocate(this%itrkout)
212  call mem_deallocate(this%itrkhdr)
213  call mem_deallocate(this%itrkcsv)
214  call mem_deallocate(this%irlstls)
215  call mem_deallocate(this%frctrn)
216  call mem_deallocate(this%iexmeth)
217  call mem_deallocate(this%ichkmeth)
218  call mem_deallocate(this%icycwin)
219  call mem_deallocate(this%extol)
220  call mem_deallocate(this%rttol)
221  call mem_deallocate(this%rtfreq)
222 
223  ! Deallocate arrays
224  call mem_deallocate(this%rptx)
225  call mem_deallocate(this%rpty)
226  call mem_deallocate(this%rptz)
227  call mem_deallocate(this%rptnode)
228  call mem_deallocate(this%rptm)
229  call mem_deallocate(this%rptname, 'RPTNAME', this%memoryPath)
230 
231  ! Deallocate period block storage
232  if (allocated(this%period_block_lines)) deallocate (this%period_block_lines)
233 
234  ! Deallocate objects
235  call this%particles%destroy(this%memoryPath)
236  call this%schedule%destroy()
237  deallocate (this%particles)
238  deallocate (this%schedule)

◆ prp_df_obs()

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

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

823  ! dummy
824  class(PrtPrpType) :: this
825  ! local
826  integer(I4B) :: indx
827  call this%obs%StoreObsType('prp', .true., indx)
828  this%obs%obsData(indx)%ProcessIdPtr => defaultobsidprocessor
829 
830  ! Store obs type and assign procedure pointer
831  ! for to-mvr observation type.
832  call this%obs%StoreObsType('to-mvr', .true., indx)
833  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 1031 of file prt-prp.f90.

1032  ! modules
1035  ! dummy variables
1036  class(PrtPrpType), intent(inout) :: this
1037  ! local variables
1038  type(PrtPrpParamFoundType) :: found
1039 
1040  call mem_set_value(this%nreleasepoints, 'NRELEASEPTS', this%input_mempath, &
1041  found%nreleasepts)
1042  call mem_set_value(this%nreleasetimes, 'NRELEASETIMES', this%input_mempath, &
1043  found%nreleasetimes)
1044 
1045  write (this%iout, '(1x,a)') 'PROCESSING PARTICLE INPUT DIMENSIONS'
1046  write (this%iout, '(4x,a,i0)') 'NRELEASEPTS = ', this%nreleasepoints
1047  write (this%iout, '(4x,a,i0)') 'NRELEASETIMES = ', this%nreleasetimes
1048  write (this%iout, '(1x,a)') 'END OF PARTICLE INPUT DIMENSIONS'
1049 
1050  ! set maxbound and nbound to nreleasepts
1051  this%maxbound = this%nreleasepoints
1052  this%nbound = this%nreleasepoints
1053 
1054  call this%prp_allocate_arrays()
1055  call this%prp_packagedata()
1056  call this%prp_releasetimes()
1057  call this%prp_load_releasetimefrequency()

◆ prp_load_releasetimefrequency()

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

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

1257  ! modules
1258  use tdismodule, only: totalsimtime
1259  ! dummy
1260  class(PrtPrpType), intent(inout) :: this
1261  ! local
1262  real(DP), allocatable :: times(:)
1263 
1264  ! check if a release time frequency is configured
1265  if (this%rtfreq <= dzero) return
1266 
1267  ! create array of regularly-spaced release times
1268  times = arange( &
1269  start=dzero, &
1270  stop=totalsimtime, &
1271  step=this%rtfreq)
1272 
1273  ! register times with release schedule
1274  call this%schedule%time_select%extend(times)
1275 
1276  ! make sure times strictly increase
1277  if (.not. this%schedule%time_select%increasing()) then
1278  errmsg = "Release times must strictly increase"
1279  call store_error(errmsg)
1280  call store_error_filename(this%input_fname)
1281  end if
1282 
1283  ! deallocate
1284  deallocate (times)
1285 
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 
)
private

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

996  ! -- modules
998  ! -- dummy variables
999  class(PrtPrpType), intent(inout) :: this
1000  type(PrtPrpParamFoundType), intent(in) :: found
1001  character(len=*), intent(in) :: trackfile
1002  character(len=*), intent(in) :: trackcsvfile
1003  ! -- local variables
1004  ! formats
1005  character(len=*), parameter :: fmttrkbin = &
1006  "(4x, 'PARTICLE TRACKS WILL BE SAVED TO BINARY FILE: ', a, /4x, &
1007  &'OPENED ON UNIT: ', I0)"
1008  character(len=*), parameter :: fmttrkcsv = &
1009  "(4x, 'PARTICLE TRACKS WILL BE SAVED TO CSV FILE: ', a, /4x, &
1010  &'OPENED ON UNIT: ', I0)"
1011 
1012  write (this%iout, '(1x,a)') 'PROCESSING PARTICLE INPUT DIMENSIONS'
1013 
1014  if (found%frctrn) then
1015  write (this%iout, '(4x,a)') &
1016  'IF DISV, TRACKING WILL USE THE TERNARY METHOD REGARDLESS OF CELL TYPE'
1017  end if
1018 
1019  if (found%trackfile) then
1020  write (this%iout, fmttrkbin) trim(adjustl(trackfile)), this%itrkout
1021  end if
1022 
1023  if (found%trackcsvfile) then
1024  write (this%iout, fmttrkcsv) trim(adjustl(trackcsvfile)), this%itrkcsv
1025  end if
1026 
1027  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 816 of file prt-prp.f90.

817  class(PrtPrpType) :: this
818  prp_obs_supported = .true.

◆ prp_options()

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

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

838  ! -- modules
841  use openspecmodule, only: access, form
844  ! -- dummy variables
845  class(PrtPrpType), intent(inout) :: this
846  ! -- local variables
847  character(len=LENVARNAME), dimension(3) :: drytrack_method = &
848  &[character(len=LENVARNAME) :: 'DROP', 'STOP', 'STAY']
849  character(len=lenvarname), dimension(2) :: coorcheck_method = &
850  &[character(len=LENVARNAME) :: 'NONE', 'EAGER']
851  character(len=LINELENGTH) :: trackfile, trackcsvfile, fname
852  type(PrtPrpParamFoundType) :: found
853  character(len=*), parameter :: fmtextolwrn = &
854  "('WARNING: EXIT_SOLVE_TOLERANCE is set to ',g10.3,' &
855  &which is much greater than the default value of ',g10.3,'. &
856  &The tolerance that strikes the best balance between accuracy &
857  &and runtime is problem-dependent. Since the variable being &
858  &solved varies from 0 to 1, tolerance values much less than 1 &
859  &typically give the best results.')"
860 
861  ! source base class options
862  call this%BndExtType%source_options()
863 
864  ! update defaults from input context
865  call mem_set_value(this%stoptime, 'STOPTIME', this%input_mempath, &
866  found%stoptime)
867  call mem_set_value(this%stoptraveltime, 'STOPTRAVELTIME', &
868  this%input_mempath, found%stoptraveltime)
869  call mem_set_value(this%istopweaksink, 'ISTOPWEAKSINK', this%input_mempath, &
870  found%istopweaksink)
871  call mem_set_value(this%istopzone, 'ISTOPZONE', this%input_mempath, &
872  found%istopzone)
873  call mem_set_value(this%drape, 'DRAPE', this%input_mempath, &
874  found%drape)
875  call mem_set_value(this%idrymeth, 'IDRYMETH', this%input_mempath, &
876  drytrack_method, found%idrymeth)
877  call mem_set_value(trackfile, 'TRACKFILE', this%input_mempath, &
878  found%trackfile)
879  call mem_set_value(trackcsvfile, 'TRACKCSVFILE', this%input_mempath, &
880  found%trackcsvfile)
881  call mem_set_value(this%localz, 'LOCALZ', this%input_mempath, &
882  found%localz)
883  call mem_set_value(this%extend, 'EXTEND', this%input_mempath, &
884  found%extend)
885  call mem_set_value(this%extol, 'EXTOL', this%input_mempath, &
886  found%extol)
887  call mem_set_value(this%rttol, 'RTTOL', this%input_mempath, &
888  found%rttol)
889  call mem_set_value(this%rtfreq, 'RTFREQ', this%input_mempath, &
890  found%rtfreq)
891  call mem_set_value(this%frctrn, 'FRCTRN', this%input_mempath, &
892  found%frctrn)
893  call mem_set_value(this%iexmeth, 'IEXMETH', this%input_mempath, &
894  found%iexmeth)
895  call mem_set_value(this%ichkmeth, 'ICHKMETH', this%input_mempath, &
896  coorcheck_method, found%ichkmeth)
897  call mem_set_value(this%icycwin, 'ICYCWIN', this%input_mempath, found%icycwin)
898 
899  ! update internal state and validate input
900  if (found%idrymeth) then
901  if (this%idrymeth == 0) then
902  write (errmsg, '(a)') 'Unsupported dry tracking method. &
903  &DRY_TRACKING_METHOD must be "DROP", "STOP", or "STAY"'
904  call store_error(errmsg)
905  else
906  ! adjust for method zero indexing
907  this%idrymeth = this%idrymeth - 1
908  end if
909  end if
910 
911  if (found%extol) then
912  if (this%extol <= dzero) &
913  call store_error('EXIT_SOLVE_TOLERANCE MUST BE POSITIVE')
914  if (this%extol > dem2) then
915  write (warnmsg, fmt=fmtextolwrn) &
916  this%extol, default_exit_solve_tolerance
917  call store_warning(warnmsg)
918  end if
919  end if
920 
921  if (found%rttol) then
922  if (this%rttol <= dzero) &
923  call store_error('RELEASE_TIME_TOLERANCE MUST BE POSITIVE')
924  end if
925 
926  if (found%rtfreq) then
927  if (this%rtfreq <= dzero) &
928  call store_error('RELEASE_TIME_FREQUENCY MUST BE POSITIVE')
929  end if
930 
931  if (found%iexmeth) then
932  if (.not. (this%iexmeth /= 1 .or. this%iexmeth /= 2)) &
933  call store_error('DEV_EXIT_SOLVE_METHOD MUST BE &
934  &1 (BRENT) OR 2 (CHANDRUPATLA)')
935  end if
936 
937  if (found%ichkmeth) then
938  if (this%ichkmeth == 0) then
939  write (errmsg, '(a)') 'Unsupported coordinate check method. &
940  &COORDINATE_CHECK_METHOD must be "NONE" or "EAGER"'
941  call store_error(errmsg)
942  else
943  ! adjust for method zero based indexing
944  this%ichkmeth = this%ichkmeth - 1
945  end if
946  end if
947 
948  if (found%icycwin) then
949  if (this%icycwin < 0) &
950  call store_error('CYCLE_DETECTION_WINDOW MUST BE NON-NEGATIVE')
951  end if
952 
953  ! fileout options
954  if (found%trackfile) then
955  this%itrkout = getunit()
956  call openfile(this%itrkout, this%iout, trackfile, 'DATA(BINARY)', &
957  form, access, filstat_opt='REPLACE', &
958  mode_opt=mnormal)
959  ! open and write ascii header spec file
960  this%itrkhdr = getunit()
961  fname = trim(trackfile)//'.hdr'
962  call openfile(this%itrkhdr, this%iout, fname, 'CSV', &
963  filstat_opt='REPLACE', mode_opt=mnormal)
964  write (this%itrkhdr, '(a,/,a)') trackheader, trackdtypes
965  end if
966 
967  if (found%trackcsvfile) then
968  this%itrkcsv = getunit()
969  call openfile(this%itrkcsv, this%iout, trackcsvfile, 'CSV', &
970  filstat_opt='REPLACE')
971  write (this%itrkcsv, '(a)') trackheader
972  end if
973 
974  ! terminate if any errors were detected
975  if (count_errors() > 0) then
976  call store_error_filename(this%input_fname)
977  end if
978 
979  ! log found options
980  call this%prp_log_options(found, trackfile, trackcsvfile)
981 
982  ! Create release schedule now that we know
983  ! the coincident release time tolerance
984  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)
private

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

1075  use memorymanagermodule, only: mem_setptr
1076  use geomutilmodule, only: get_node
1078  ! dummy
1079  class(PrtPrpType), intent(inout) :: this
1080  ! local
1081  integer(I4B), dimension(:), pointer, contiguous :: irptno
1082  integer(I4B), dimension(:, :), pointer, contiguous :: cellids
1083  real(DP), dimension(:), pointer, contiguous :: xrpts, yrpts, zrpts
1084  type(CharacterStringType), dimension(:), pointer, &
1085  contiguous :: boundnames
1086  character(len=LENBOUNDNAME) :: bndName, bndNameTemp
1087  character(len=9) :: cno
1088  character(len=20) :: cellidstr
1089  integer(I4B), dimension(:), allocatable :: nboundchk
1090  integer(I4B), dimension(:), pointer :: cellid
1091  integer(I4B) :: n, noder, nodeu, rptno
1092 
1093  ! set input context pointers
1094  call mem_setptr(irptno, 'IRPTNO', this%input_mempath)
1095  call mem_setptr(cellids, 'CELLID', this%input_mempath)
1096  call mem_setptr(xrpts, 'XRPT', this%input_mempath)
1097  call mem_setptr(yrpts, 'YRPT', this%input_mempath)
1098  call mem_setptr(zrpts, 'ZRPT', this%input_mempath)
1099  call mem_setptr(boundnames, 'BOUNDNAME', this%input_mempath)
1100 
1101  ! allocate and initialize temporary variables
1102  allocate (nboundchk(this%nreleasepoints))
1103  do n = 1, this%nreleasepoints
1104  nboundchk(n) = 0
1105  end do
1106 
1107  write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%packName)) &
1108  //' PACKAGEDATA'
1109 
1110  do n = 1, size(irptno)
1111 
1112  rptno = irptno(n)
1113 
1114  if (rptno < 1 .or. rptno > this%nreleasepoints) then
1115  write (errmsg, '(a,i0,a,i0,a)') &
1116  'Expected ', this%nreleasepoints, ' release points. &
1117  &Points must be numbered from 1 to ', this%nreleasepoints, '.'
1118  call store_error(errmsg)
1119  cycle
1120  end if
1121 
1122  ! increment nboundchk
1123  nboundchk(rptno) = nboundchk(rptno) + 1
1124 
1125  ! set cellid
1126  cellid => cellids(:, n)
1127 
1128  ! set node user
1129  if (this%dis%ndim == 1) then
1130  nodeu = cellid(1)
1131  elseif (this%dis%ndim == 2) then
1132  nodeu = get_node(cellid(1), 1, cellid(2), &
1133  this%dis%mshape(1), 1, &
1134  this%dis%mshape(2))
1135  else
1136  nodeu = get_node(cellid(1), cellid(2), cellid(3), &
1137  this%dis%mshape(1), &
1138  this%dis%mshape(2), &
1139  this%dis%mshape(3))
1140  end if
1141 
1142  ! set noder
1143  noder = this%dis%get_nodenumber(nodeu, 1)
1144  if (noder <= 0) then
1145  call this%dis%nodeu_to_string(nodeu, cellidstr)
1146  write (errmsg, '(a)') &
1147  'Particle release point configured for nonexistent cell: '// &
1148  trim(adjustl(cellidstr))//'. This cell has IDOMAIN <= 0 and '&
1149  &'therefore does not exist in the model grid.'
1150  call store_error(errmsg)
1151  cycle
1152  else
1153  this%rptnode(rptno) = noder
1154  end if
1155 
1156  if (this%localz .and. (zrpts(n) < 0 .or. zrpts(n) > 1)) then
1157  call store_error('Local z coordinate must fall in the interval [0, 1]')
1158  cycle
1159  end if
1160 
1161  ! set coordinates
1162  this%rptx(rptno) = xrpts(n)
1163  this%rpty(rptno) = yrpts(n)
1164  this%rptz(rptno) = zrpts(n)
1165 
1166  ! set default boundname
1167  write (cno, '(i9.9)') rptno
1168  bndname = 'PRP'//cno
1169 
1170  ! read boundnames from file, if provided
1171  if (this%inamedbound /= 0) then
1172  bndnametemp = boundnames(n)
1173  if (bndnametemp /= '') bndname = bndnametemp
1174  else
1175  bndname = ''
1176  end if
1177 
1178  ! set boundname
1179  this%rptname(rptno) = bndname
1180  end do
1181 
1182  write (this%iout, '(1x,a)') &
1183  'END OF '//trim(adjustl(this%packName))//' PACKAGEDATA'
1184 
1185  ! check for duplicate or missing particle release points
1186  do n = 1, this%nreleasepoints
1187  if (nboundchk(n) == 0) then
1188  write (errmsg, '(a,a,1x,i0,a)') 'No data specified for particle ', &
1189  'release point', n, '.'
1190  call store_error(errmsg)
1191  else if (nboundchk(n) > 1) then
1192  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
1193  'Data for particle release point', n, 'specified', nboundchk(n), &
1194  'times.'
1195  call store_error(errmsg)
1196  end if
1197  end do
1198 
1199  ! terminate if any errors were detected
1200  if (count_errors() > 0) then
1201  call store_error_filename(this%input_fname)
1202  end if
1203 
1204  ! cleanup
1205  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
Here is the call graph for this function:

◆ prp_releasetimes()

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

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

1211  ! dummy
1212  class(PrtPrpType), intent(inout) :: this
1213  ! local
1214  real(DP), dimension(:), pointer, contiguous :: time
1215  integer(I4B) :: n, isize
1216  real(DP), allocatable :: times(:)
1217 
1218  if (this%nreleasetimes <= 0) return
1219 
1220  ! allocate times array
1221  allocate (times(this%nreleasetimes))
1222 
1223  ! check if input array was read
1224  call get_isize('TIME', this%input_mempath, isize)
1225 
1226  if (isize <= 0) then
1227  errmsg = "RELEASTIMES block expected when &
1228  &NRELEASETIMES dimension is non-zero."
1229  call store_error(errmsg)
1230  call store_error_filename(this%input_fname)
1231  end if
1232 
1233  ! set input context pointer
1234  call mem_setptr(time, 'TIME', this%input_mempath)
1235 
1236  ! set input data
1237  do n = 1, size(time)
1238  times(n) = time(n)
1239  end do
1240 
1241  ! register times with the release schedule
1242  call this%schedule%time_select%extend(times)
1243 
1244  ! make sure times strictly increase
1245  if (.not. this%schedule%time_select%increasing()) then
1246  errmsg = "RELEASTIMES block entries must strictly increase."
1247  call store_error(errmsg)
1248  call store_error_filename(this%input_fname)
1249  end if
1250 
1251  ! deallocate
1252  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 724 of file prt-prp.f90.

725  ! modules
726  use tdismodule, only: kper, nper
729  ! dummy variables
730  class(PrtPrpType), intent(inout) :: this
731  ! local variables
732  type(CharacterStringType), dimension(:), contiguous, &
733  pointer :: settings
734  integer(I4B), pointer :: iper, ionper, nlist
735  integer(I4B) :: n
736 
737  ! set pointer to last and next period loaded
738  call mem_setptr(iper, 'IPER', this%input_mempath)
739  call mem_setptr(ionper, 'IONPER', this%input_mempath)
740 
741  if (kper == 1 .and. &
742  (iper == 0) .and. &
743  (ionper > nper) .and. &
744  size(this%schedule%time_select%times) == 0) then
745  ! If the user hasn't provided any release settings (neither
746  ! explicit release times, release time frequency, or period
747  ! block release settings), default to a single release at the
748  ! start of the first period's first time step.
749  ! Store default configuration; advance() will be called in prp_ad().
750  if (allocated(this%period_block_lines)) deallocate (this%period_block_lines)
751  allocate (this%period_block_lines(1))
752  this%period_block_lines(1) = "FIRST"
753  return
754  else if (iper /= kper) then
755  return
756  end if
757 
758  ! set input context pointers
759  call mem_setptr(nlist, 'NBOUND', this%input_mempath)
760  call mem_setptr(settings, 'SETTING', this%input_mempath)
761 
762  ! Store period block configuration for fill-forward.
763  if (allocated(this%period_block_lines)) deallocate (this%period_block_lines)
764  allocate (this%period_block_lines(nlist))
765  do n = 1, nlist
766  this%period_block_lines(n) = settings(n)
767  end do
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 242 of file prt-prp.f90.

243  class(PrtPrpType) :: this
244  integer(I4B), dimension(:), pointer, contiguous :: ibound
245  integer(I4B), dimension(:), pointer, contiguous :: izone
246 
247  this%ibound => ibound
248  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 583 of file prt-prp.f90.

584  ! dummy
585  class(PrtPrpType), intent(inout) :: this !< this instance
586  integer(I4B), intent(in) :: ip !< particle index
587  real(DP), intent(in) :: trelease !< release time
588  ! local
589  integer(I4B) :: np
590  type(ParticleType), pointer :: particle
591 
592  call this%initialize_particle(particle, ip, trelease)
593  np = this%nparticles + 1
594  this%nparticles = np
595  call this%particles%put(particle, np)
596  deallocate (particle)
597  this%rptm(ip) = this%rptm(ip) + done ! TODO configurable mass
598 

◆ 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 540 of file prt-prp.f90.

541  class(PrtPrpType), intent(inout) :: this !< this instance
542  integer(I4B), intent(in) :: ic !< cell index
543  real(DP), intent(in) :: x, y, z !< release point
544  ! local
545  real(DP), allocatable :: polyverts(:, :)
546 
547  call this%fmi%dis%get_polyverts(ic, polyverts)
548  if (.not. point_in_polygon(x, y, polyverts)) then
549  write (errmsg, '(a,g0,a,g0,a,i0)') &
550  'Error: release point (x=', x, ', y=', y, ') is not in cell ', &
551  this%dis%get_nodeuser(ic)
552  call store_error(errmsg, terminate=.false.)
553  call store_error_filename(this%input_fname)
554  end if
555  if (z > maxval(this%dis%top)) then
556  write (errmsg, '(a,g0,a,g0,a,i0)') &
557  'Error: release point (z=', z, ') is above grid top ', &
558  maxval(this%dis%top)
559  call store_error(errmsg, terminate=.false.)
560  call store_error_filename(this%input_fname)
561  else if (z < minval(this%dis%bot)) then
562  write (errmsg, '(a,g0,a,g0,a,i0)') &
563  'Error: release point (z=', z, ') is below grid bottom ', &
564  minval(this%dis%bot)
565  call store_error(errmsg, terminate=.false.)
566  call store_error_filename(this%input_fname)
567  end if
568  deallocate (polyverts)
Here is the call graph for this function:

Variable Documentation

◆ default_exit_solve_tolerance

real(dp), parameter prtprpmodule::default_exit_solve_tolerance = DEM5
private

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

36  real(DP), parameter :: DEFAULT_EXIT_SOLVE_TOLERANCE = dem5

◆ 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'