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 prp_commit (this)
 Commit staged particle state to the "final" store. 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 824 of file prt-prp.f90.

825  class(PrtPrpType), intent(inout) :: this
826  ! not implemented, not used

◆ exg_prp_ad()

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

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

515  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 396 of file prt-prp.f90.

399  class(ExgPrtPrpType) :: this
400  integer(I4B), dimension(:), pointer, contiguous, optional :: nodelist
401  real(DP), dimension(:, :), pointer, contiguous, optional :: auxvar
402  ! local
403  integer(I4B), dimension(:, :), pointer, contiguous :: cellid
404  integer(I4B), dimension(:), pointer, contiguous :: nodeulist
405  type(CharacterStringType), dimension(:), pointer, contiguous :: boundname
406  real(DP), dimension(:, :), pointer, contiguous :: auxvar_input
407 
408  call mem_allocate(cellid, this%dis%ndim, 0, 'CELLID', this%input_mempath)
409  call mem_allocate(nodeulist, 0, 'NODEULIST', this%input_mempath)
410  call mem_allocate(boundname, lenboundname, 0, 'BOUNDNAME', &
411  this%input_mempath)
412  call mem_allocate(auxvar_input, 0, 0, 'AUXVAR', this%input_mempath)
413 
414  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 373 of file prt-prp.f90.

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

◆ exg_prp_ar()

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

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

419  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 531 of file prt-prp.f90.

532  use budgetmodule, only: budgettype
533  class(ExgPrtPrpType) :: this
534  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 522 of file prt-prp.f90.

523  class(ExgPrtPrpType) :: this
524  real(DP), dimension(:), intent(in) :: hnew
525  real(DP), dimension(:), intent(inout) :: flowja
526  integer(I4B), intent(in) :: imover

◆ exg_prp_dimensions()

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

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

1077  class(ExgPrtPrpType), intent(inout) :: this
1078 
1079  this%nreleasepoints = 0
1080  this%nreleasetimes = 0
1081  this%maxbound = 0
1082  this%nbound = 0
1083 
1084  call this%prp_allocate_arrays()

◆ exg_prp_options()

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

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

1004  class(ExgPrtPrpType), intent(inout) :: this
1005  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 539 of file prt-prp.f90.

540  class(ExgPrtPrpType) :: this
541  integer(I4B), intent(in) :: icbcfl
542  integer(I4B), intent(in) :: ibudfl
543  integer(I4B), intent(in) :: icbcun
544  integer(I4B), dimension(:), optional, intent(in) :: imap

◆ exg_prp_rp()

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

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

786  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 624 of file prt-prp.f90.

626  class(PrtPrpType), intent(inout) :: this !< this instance
627  type(ParticleType), pointer, intent(inout) :: particle !< the particle
628  integer(I4B), intent(in) :: ip !< particle index
629  real(DP), intent(in) :: trelease !< release time
630  ! local
631  logical(LGP) :: draped
632  integer(I4B) :: irow, icol, ilay, icpl
633  integer(I4B) :: ic, icu, ic_old
634  real(DP) :: x, y, z
635  ! formats
636  character(len=*), parameter :: fmticterr = &
637  "('Error in ',a,': Flow model interface does not contain ICELLTYPE. &
638  &ICELLTYPE is required for PRT to distinguish convertible cells &
639  &from confined cells if LOCAL_Z release coordinates are provided. &
640  &Make sure a GWFGRID entry is configured in the PRT FMI package.')"
641 
642  ic = this%rptnode(ip)
643 
644  call create_particle(particle)
645 
646  if (size(this%boundname) /= 0) then
647  particle%name = this%boundname(ip)
648  else
649  particle%name = ''
650  end if
651 
652  particle%irpt = ip
653  particle%istopweaksink = this%istopweaksink
654  particle%istopzone = this%istopzone
655  particle%idrymeth = this%idrymeth
656  particle%istatus = 0 ! status 0 until tracking starts
657 
658  ! If the cell is inactive, either drape the particle
659  ! to the top-most active cell beneath it if drape is
660  ! enabled, or else terminate permanently unreleased.
661  draped = .false.
662  if (this%ibound(ic) == 0) then
663  ic_old = ic
664  if (this%drape) then
665  call this%dis%highest_active(ic, this%ibound)
666  draped = ic /= ic_old
667  if (.not. draped .and. this%ibound(ic) == 0) then
668  ! negative unreleased status signals to the
669  ! tracking method that we haven't yet saved
670  ! a termination record, it needs to do so.
671  particle%istatus = -1 * term_unreleased
672  end if
673  else
674  particle%istatus = -1 * term_unreleased
675  end if
676  end if
677 
678  icu = this%dis%get_nodeuser(ic)
679  particle%icu = icu
680  select type (dis => this%dis)
681  type is (distype)
682  call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, irow, icol, ilay)
683  type is (disvtype)
684  call get_jk(icu, dis%ncpl, dis%nlay, icpl, ilay)
685  end select
686  particle%ilay = ilay
687  particle%izone = this%rptzone(ic)
688 
689  ! if the particle was draped to this cell, set the z coord to
690  ! the effective top of the cell. if it was not draped, and is
691  ! a local z coord, calculate the corresponding model z coord.
692  if (draped) then
693  z = this%fmi%dis%bot(ic) + &
694  this%fmi%gwfsat(ic) * &
695  (this%fmi%dis%top(ic) - this%fmi%dis%bot(ic))
696  else if (this%localz) then
697  z = this%fmi%dis%bot(ic) + &
698  this%rptz(ip) * &
699  this%fmi%gwfsat(ic) * &
700  (this%fmi%dis%top(ic) - this%fmi%dis%bot(ic))
701  else
702  z = this%rptz(ip)
703  end if
704 
705  x = this%rptx(ip)
706  y = this%rpty(ip)
707 
708  if (this%ichkmeth > 0) &
709  call this%validate_release_point(ic, x, y, z)
710 
711  particle%x = x
712  particle%y = y
713  particle%z = z
714  particle%trelease = trelease
715 
716  ! Set stop time to earlier of STOPTIME and STOPTRAVELTIME
717  if (this%stoptraveltime == huge(1d0)) then
718  particle%tstop = this%stoptime
719  else
720  particle%tstop = particle%trelease + this%stoptraveltime
721  if (this%stoptime < particle%tstop) particle%tstop = this%stoptime
722  end if
723 
724  particle%ttrack = particle%trelease
725  particle%itrdomain(level_model) = 0
726  particle%iboundary(level_model) = 0
727  particle%itrdomain(level_feature) = ic
728  particle%iboundary(level_feature) = 0
729  particle%itrdomain(level_subfeature) = 0
730  particle%iboundary(level_subfeature) = 0
731  particle%frctrn = this%frctrn
732  particle%iexmeth = this%iexmeth
733  particle%extend = this%extend
734  particle%icycwin = this%icycwin
735  particle%extol = this%extol
@ term_unreleased
terminated permanently unreleased
Definition: Particle.f90:38
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 548 of file prt-prp.f90.

549  class(PrtPrpType), intent(inout) :: this !< prp
550  if (this%iprpak > 0) then
551  write (this%iout, "(1x,/1x,a,1x,i0)") &
552  'PARTICLE RELEASE FOR PRP', this%ibcnum
553  call this%schedule%log(this%iout)
554  end if

◆ prp_ad()

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

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

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

257  ! dummy
258  class(PrtPrpType) :: this
259  integer(I4B), dimension(:), pointer, contiguous, optional :: nodelist
260  real(DP), dimension(:, :), pointer, contiguous, optional :: auxvar
261  ! local
262  integer(I4B) :: nps
263 
264  call this%BndExtType%allocate_arrays()
265 
266  ! Allocate particle stores
267  call create_particle_store( &
268  this%particles, 0, &
269  this%memoryPath)
270  call create_particle_store( &
271  this%particles_staging, 0, &
272  trim(this%memoryPath)//'-STAGING')
273 
274  ! Allocate arrays
275  call mem_allocate(this%rptx, this%nreleasepoints, 'RPTX', this%memoryPath)
276  call mem_allocate(this%rpty, this%nreleasepoints, 'RPTY', this%memoryPath)
277  call mem_allocate(this%rptz, this%nreleasepoints, 'RPTZ', this%memoryPath)
278  call mem_allocate(this%rptm, this%nreleasepoints, 'RPTMASS', &
279  this%memoryPath)
280  call mem_allocate(this%rptnode, this%nreleasepoints, 'RPTNODER', &
281  this%memoryPath)
282  call mem_allocate(this%rptname, lenboundname, this%nreleasepoints, &
283  'RPTNAME', this%memoryPath)
284 
285  ! Initialize arrays
286  do nps = 1, this%nreleasepoints
287  this%rptm(nps) = dzero
288  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 292 of file prt-prp.f90.

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

◆ prp_ar()

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

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

353  class(PrtPrpType), intent(inout) :: this
354  integer(I4B) :: n
355 
356  call this%obs%obs_ar()
357 
358  if (this%inamedbound /= 0) then
359  do n = 1, this%nreleasepoints
360  this%boundname(n) = this%rptname(n)
361  end do
362  end if
363  do n = 1, this%nreleasepoints
364  this%nodelist(n) = this%rptnode(n)
365  end do

◆ prp_commit()

subroutine prtprpmodule::prp_commit ( class(prtprptype this)

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

506  class(PrtPrpType) :: this
507  call this%particles%resize( &
508  this%particles_staging%num_stored(), &
509  this%memoryPath)
510  call this%particles%copy_from(this%particles_staging)

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

791  ! modules
792  use tdismodule, only: delt
793  ! dummy variables
794  class(PrtPrpType) :: this
795  real(DP), dimension(:), intent(in) :: hnew
796  real(DP), dimension(:), intent(inout) :: flowja !< flow between package and model
797  integer(I4B), intent(in) :: imover !< flag indicating if the mover package is active
798  ! local variables
799  integer(I4B) :: i
800  integer(I4B) :: node
801  integer(I4B) :: idiag
802  real(DP) :: rrate
803 
804  ! If no boundaries, skip flow calculations.
805  if (this%nbound <= 0) return
806 
807  ! Loop through each boundary calculating flow.
808  do i = 1, this%nbound
809  node = this%nodelist(i)
810  rrate = dzero
811  ! If cell is no-flow or constant-head, then ignore it.
812  if (node > 0) then
813  ! Calculate the flow rate into the cell.
814  idiag = this%dis%con%ia(node)
815  rrate = this%rptm(i) * (done / delt) ! reciprocal of tstp length
816  flowja(idiag) = flowja(idiag) + rrate
817  end if
818 
819  ! Save simulated value to simvals array.
820  this%simvals(i) = rrate
821  end do
real(dp), pointer, public delt
length of the current time step
Definition: tdis.f90:32

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

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

◆ prp_da()

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

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

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

◆ prp_df_obs()

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

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

837  ! dummy
838  class(PrtPrpType) :: this
839  ! local
840  integer(I4B) :: indx
841  call this%obs%StoreObsType('prp', .true., indx)
842  this%obs%obsData(indx)%ProcessIdPtr => defaultobsidprocessor
843 
844  ! Store obs type and assign procedure pointer
845  ! for to-mvr observation type.
846  call this%obs%StoreObsType('to-mvr', .true., indx)
847  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 1045 of file prt-prp.f90.

1046  ! modules
1049  ! dummy variables
1050  class(PrtPrpType), intent(inout) :: this
1051  ! local variables
1052  type(PrtPrpParamFoundType) :: found
1053 
1054  call mem_set_value(this%nreleasepoints, 'NRELEASEPTS', this%input_mempath, &
1055  found%nreleasepts)
1056  call mem_set_value(this%nreleasetimes, 'NRELEASETIMES', this%input_mempath, &
1057  found%nreleasetimes)
1058 
1059  write (this%iout, '(1x,a)') 'PROCESSING PARTICLE INPUT DIMENSIONS'
1060  write (this%iout, '(4x,a,i0)') 'NRELEASEPTS = ', this%nreleasepoints
1061  write (this%iout, '(4x,a,i0)') 'NRELEASETIMES = ', this%nreleasetimes
1062  write (this%iout, '(1x,a)') 'END OF PARTICLE INPUT DIMENSIONS'
1063 
1064  ! set maxbound and nbound to nreleasepts
1065  this%maxbound = this%nreleasepoints
1066  this%nbound = this%nreleasepoints
1067 
1068  call this%prp_allocate_arrays()
1069  call this%prp_packagedata()
1070  call this%prp_releasetimes()
1071  call this%prp_load_releasetimefrequency()

◆ prp_load_releasetimefrequency()

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

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

1279  ! modules
1280  use tdismodule, only: totalsimtime
1281  ! dummy
1282  class(PrtPrpType), intent(inout) :: this
1283  ! local
1284  real(DP), allocatable :: times(:)
1285 
1286  ! check if a release time frequency is configured
1287  if (this%rtfreq <= dzero) return
1288 
1289  ! create array of regularly-spaced release times
1290  times = arange( &
1291  start=dzero, &
1292  stop=totalsimtime, &
1293  step=this%rtfreq)
1294 
1295  ! register times with release schedule
1296  call this%schedule%time_select%extend(times)
1297 
1298  ! make sure times strictly increase
1299  if (.not. this%schedule%time_select%increasing()) then
1300  errmsg = "Release times must strictly increase"
1301  call store_error(errmsg)
1302  call store_error_filename(this%input_fname)
1303  end if
1304 
1305  ! deallocate
1306  deallocate (times)
1307 
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 1009 of file prt-prp.f90.

1010  ! -- modules
1012  ! -- dummy variables
1013  class(PrtPrpType), intent(inout) :: this
1014  type(PrtPrpParamFoundType), intent(in) :: found
1015  character(len=*), intent(in) :: trackfile
1016  character(len=*), intent(in) :: trackcsvfile
1017  ! -- local variables
1018  ! formats
1019  character(len=*), parameter :: fmttrkbin = &
1020  "(4x, 'PARTICLE TRACKS WILL BE SAVED TO BINARY FILE: ', a, /4x, &
1021  &'OPENED ON UNIT: ', I0)"
1022  character(len=*), parameter :: fmttrkcsv = &
1023  "(4x, 'PARTICLE TRACKS WILL BE SAVED TO CSV FILE: ', a, /4x, &
1024  &'OPENED ON UNIT: ', I0)"
1025 
1026  write (this%iout, '(1x,a)') 'PROCESSING PARTICLE INPUT DIMENSIONS'
1027 
1028  if (found%frctrn) then
1029  write (this%iout, '(4x,a)') &
1030  'IF DISV, TRACKING WILL USE THE TERNARY METHOD REGARDLESS OF CELL TYPE'
1031  end if
1032 
1033  if (found%trackfile) then
1034  write (this%iout, fmttrkbin) trim(adjustl(trackfile)), this%itrkout
1035  end if
1036 
1037  if (found%trackcsvfile) then
1038  write (this%iout, fmttrkcsv) trim(adjustl(trackcsvfile)), this%itrkcsv
1039  end if
1040 
1041  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 830 of file prt-prp.f90.

831  class(PrtPrpType) :: this
832  prp_obs_supported = .true.

◆ prp_options()

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

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

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

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

◆ prp_releasetimes()

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

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

1233  ! dummy
1234  class(PrtPrpType), intent(inout) :: this
1235  ! local
1236  real(DP), dimension(:), pointer, contiguous :: time
1237  integer(I4B) :: n, isize
1238  real(DP), allocatable :: times(:)
1239 
1240  if (this%nreleasetimes <= 0) return
1241 
1242  ! allocate times array
1243  allocate (times(this%nreleasetimes))
1244 
1245  ! check if input array was read
1246  call get_isize('TIME', this%input_mempath, isize)
1247 
1248  if (isize <= 0) then
1249  errmsg = "RELEASTIMES block expected when &
1250  &NRELEASETIMES dimension is non-zero."
1251  call store_error(errmsg)
1252  call store_error_filename(this%input_fname)
1253  end if
1254 
1255  ! set input context pointer
1256  call mem_setptr(time, 'TIME', this%input_mempath)
1257 
1258  ! set input data
1259  do n = 1, size(time)
1260  times(n) = time(n)
1261  end do
1262 
1263  ! register times with the release schedule
1264  call this%schedule%time_select%extend(times)
1265 
1266  ! make sure times strictly increase
1267  if (.not. this%schedule%time_select%increasing()) then
1268  errmsg = "RELEASTIMES block entries must strictly increase."
1269  call store_error(errmsg)
1270  call store_error_filename(this%input_fname)
1271  end if
1272 
1273  ! deallocate
1274  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 739 of file prt-prp.f90.

740  ! modules
741  use tdismodule, only: kper, nper
744  ! dummy variables
745  class(PrtPrpType), intent(inout) :: this
746  ! local variables
747  type(CharacterStringType), dimension(:), contiguous, &
748  pointer :: settings
749  integer(I4B), pointer :: iper, ionper, nlist
750  integer(I4B) :: n
751 
752  ! set pointer to last and next period loaded
753  call mem_setptr(iper, 'IPER', this%input_mempath)
754  call mem_setptr(ionper, 'IONPER', this%input_mempath)
755 
756  if (kper == 1 .and. &
757  (iper == 0) .and. &
758  (ionper > nper) .and. &
759  size(this%schedule%time_select%times) == 0) then
760  ! If the user hasn't provided any release settings (neither
761  ! explicit release times, release time frequency, nor period
762  ! block release settings), default to a single release at the
763  ! start of the simulation (t=0). Add t=0 directly to the time
764  ! selection rather than time step selection because the latter
765  ! would fill forward, releasing at the start of every period.
766  call this%schedule%time_select%extend([dzero])
767  return
768  else if (iper /= kper) then
769  return
770  end if
771 
772  ! set input context pointers
773  call mem_setptr(nlist, 'NBOUND', this%input_mempath)
774  call mem_setptr(settings, 'SETTING', this%input_mempath)
775 
776  ! Store period block configuration for fill-forward.
777  if (allocated(this%period_block_lines)) deallocate (this%period_block_lines)
778  allocate (this%period_block_lines(nlist))
779  do n = 1, nlist
780  this%period_block_lines(n) = settings(n)
781  end do
integer(i4b), pointer, public nper
number of stress period
Definition: tdis.f90:24

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

247  class(PrtPrpType) :: this
248  integer(I4B), dimension(:), pointer, contiguous :: ibound
249  integer(I4B), dimension(:), pointer, contiguous :: izone
250 
251  this%ibound => ibound
252  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 606 of file prt-prp.f90.

607  ! dummy
608  class(PrtPrpType), intent(inout) :: this !< this instance
609  integer(I4B), intent(in) :: ip !< particle index
610  real(DP), intent(in) :: trelease !< release time
611  ! local
612  integer(I4B) :: np
613  type(ParticleType), pointer :: particle
614 
615  call this%initialize_particle(particle, ip, trelease)
616  np = this%nparticles + 1
617  this%nparticles = np
618  call this%particles_staging%put(particle, np)
619  deallocate (particle)
620  this%rptm(ip) = this%rptm(ip) + done ! TODO configurable mass
621 

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

564  class(PrtPrpType), intent(inout) :: this !< this instance
565  integer(I4B), intent(in) :: ic !< cell index
566  real(DP), intent(in) :: x, y, z !< release point
567  ! local
568  real(DP), allocatable :: polyverts(:, :)
569 
570  call this%fmi%dis%get_polyverts(ic, polyverts)
571  if (.not. point_in_polygon(x, y, polyverts)) then
572  write (errmsg, '(a,g0,a,g0,a,i0)') &
573  'Error: release point (x=', x, ', y=', y, ') is not in cell ', &
574  this%dis%get_nodeuser(ic)
575  call store_error(errmsg, terminate=.false.)
576  call store_error_filename(this%input_fname)
577  end if
578  if (z > maxval(this%dis%top)) then
579  write (errmsg, '(a,g0,a,g0,a,i0)') &
580  'Error: release point (z=', z, ') is above grid top ', &
581  maxval(this%dis%top)
582  call store_error(errmsg, terminate=.false.)
583  call store_error_filename(this%input_fname)
584  else if (z < minval(this%dis%bot)) then
585  write (errmsg, '(a,g0,a,g0,a,i0)') &
586  'Error: release point (z=', z, ') is below grid bottom ', &
587  minval(this%dis%bot)
588  call store_error(errmsg, terminate=.false.)
589  call store_error_filename(this%input_fname)
590  end if
591  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'