prtprpmodule Module Reference

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

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

Variables
character(len=lenftype)	ftype = 'PRP'
character(len=16)	text = ' PRP'
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 665 of file prt-prp.f90.

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

initialize_particle()

subroutine prtprpmodule::initialize_particle ( class(prtprptype), intent(inout)  this, type(particletype), intent(inout), pointer  particle, integer(i4b), intent(in)  ip, real(dp), intent(in)  trelease  )

private

Parameters

[in,out]	this	this instance
[in,out]	particle	the particle
[in]	ip	particle index
[in]	trelease	release time

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

    use particlemodule, only: term_unreleased
    class(PrtPrpType), intent(inout) :: this !< this instance
    type(ParticleType), pointer, intent(inout) :: particle !< the particle
    integer(I4B), intent(in) :: ip !< particle index
    real(DP), intent(in) :: trelease !< release time
    ! local
    integer(I4B) :: irow, icol, ilay, icpl
    integer(I4B) :: ic, icu, ic_old
    real(DP) :: x, y, z
    real(DP) :: top, bot, hds
    ! formats
    character(len=*), parameter :: fmticterr = &
      "('Error in ',a,': Flow model interface does not contain ICELLTYPE. &
      &ICELLTYPE is required for PRT to distinguish convertible cells &
      &from confined cells if LOCAL_Z release coordinates are provided. &
      &Make sure a GWFGRID entry is configured in the PRT FMI package.')"
 
    ic = this%rptnode(ip)
    icu = this%dis%get_nodeuser(ic)
 
    call create_particle(particle)
 
    if (size(this%boundname) /= 0) then
      particle%name = this%boundname(ip)
    else
      particle%name = ''
    end if
 
    particle%irpt = ip
    particle%istopweaksink = this%istopweaksink
    particle%istopzone = this%istopzone
    particle%idrymeth = this%idrymeth
    particle%icu = icu
 
    select type (dis => this%dis)
    type is (distype)
      call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, irow, icol, ilay)
    type is (disvtype)
      call get_jk(icu, dis%ncpl, dis%nlay, icpl, ilay)
    end select
    particle%ilay = ilay
    particle%izone = this%rptzone(ic)
    particle%istatus = 0 ! status 0 until tracking starts
    ! If the cell is inactive, either drape the particle
    ! to the top-most active cell beneath it if drape is
    ! enabled, or else terminate permanently unreleased.
    if (this%ibound(ic) == 0) then
      ic_old = ic
      if (this%drape) then
        call this%dis%highest_active(ic, this%ibound)
        if (ic == ic_old .or. this%ibound(ic) == 0) then
          ! negative unreleased status signals to the
          ! tracking method that we haven't yet saved
          ! a termination record, it needs to do so.
          particle%istatus = -1 * term_unreleased
        end if
      else
        particle%istatus = -1 * term_unreleased
      end if
    end if
 
    ! load coordinates
    x = this%rptx(ip)
    y = this%rpty(ip)
    if (this%localz) then
      ! make sure FMI has cell type array. we need
      ! it to distinguish convertible and confined
      ! cells if release z coordinates are local
      if (this%fmi%igwfceltyp /= 1) then
        write (errmsg, fmticterr) trim(this%text)
        call store_error(errmsg, terminate=.true.)
      end if
 
      ! calculate model z coord from local z coord.
      ! if cell is confined (icelltype == 0) use the
      ! actual cell height (geometric top - bottom).
      ! otherwise use head as cell top, clamping to
      ! the cell bottom if head is below the bottom
      top = this%fmi%dis%top(ic)
      bot = this%fmi%dis%bot(ic)
      hds = this%fmi%gwfhead(ic)
      if (this%fmi%gwfceltyp(icu) /= 0) top = hds
      if (top < bot) top = bot
      z = bot + this%rptz(ip) * (top - bot)
    else
      z = this%rptz(ip)
    end if
 
    if (this%ichkmeth > 0) &
      call this%validate_release_point(ic, x, y, z)
 
    particle%x = x
    particle%y = y
    particle%z = z
    particle%trelease = trelease
 
    ! Set stop time to earlier of STOPTIME and STOPTRAVELTIME
    if (this%stoptraveltime == huge(1d0)) then
      particle%tstop = this%stoptime
    else
      particle%tstop = particle%trelease + this%stoptraveltime
      if (this%stoptime < particle%tstop) particle%tstop = this%stoptime
    end if
 
    particle%ttrack = particle%trelease
    particle%itrdomain(level_model) = 0
    particle%iboundary(level_model) = 0
    particle%itrdomain(level_feature) = ic
    particle%iboundary(level_feature) = 0
    particle%itrdomain(level_subfeature) = 0
    particle%iboundary(level_subfeature) = 0
    particle%frctrn = this%frctrn
    particle%iexmeth = this%iexmeth
    particle%extend = this%extend
    particle%icycwin = this%icycwin
    particle%extol = this%extol

Here is the call graph for this function:

log_release()

subroutine prtprpmodule::log_release

(

class(prtprptype), intent(inout)

this

)

Parameters

[in,out]

this

prp

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

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

prp_ad()

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

private

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

    use tdismodule, only: totalsimtime
    class(PrtPrpType) :: this
    integer(I4B) :: ip, it
    real(DP) :: t
 
    ! Notes
    ! -----
    ! Each release point can be thought of as
    ! a gumball machine with infinite supply:
    ! a point can release an arbitrary number
    ! of particles, but only one at any time.
    ! Coincident release times are merged to
    ! a single time by the release scheduler.
 
    ! Reset mass accumulators for this time step.
    do ip = 1, this%nreleasepoints
      this%rptm(ip) = dzero
    end do
 
    ! Advance the release schedule and check if
    ! any releases will be made this time step.
    call this%schedule%advance()
    if (.not. this%schedule%any()) return
 
    ! Log the schedule to the list file.
    call this%log_release()
 
    ! Expand the particle store. We know from the
    ! schedule how many particles will be released.
    call this%particles%resize( &
      this%particles%num_stored() + &
      (this%nreleasepoints * this%schedule%count()), &
      this%memoryPath)
 
    ! Release a particle from each point for
    ! each release time in the current step.
    do ip = 1, this%nreleasepoints
      do it = 1, this%schedule%count()
        t = this%schedule%times(it)
        ! Skip the release time if it's before the simulation
        ! starts, or if no `extend_tracking`, after it ends.
        if (t < dzero) then
          write (warnmsg, '(a,g0,a)') &
            'Skipping negative release time (t=', t, ').'
          call store_warning(warnmsg)
          cycle
        else if (t > totalsimtime .and. .not. this%extend) then
          write (warnmsg, '(a,g0,a)') &
            'Skipping release time falling after the end of the &
            &simulation (t=', t, '). Enable EXTEND_TRACKING to &
            &release particles after the simulation end time.'
          call store_warning(warnmsg)
          cycle
        end if
        call this%release(ip, t)
      end do
    end do

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

    ! dummy
    class(PrtPrpType) :: this
    integer(I4B), dimension(:), pointer, contiguous, optional :: nodelist
    real(DP), dimension(:, :), pointer, contiguous, optional :: auxvar
    ! local
    integer(I4B) :: nps
 
    call this%BndExtType%allocate_arrays()
 
    ! Allocate particle store, starting with the number
    ! of release points (arrays resized if/when needed)
    call create_particle_store( &
      this%particles, &
      this%nreleasepoints, &
      this%memoryPath)
 
    ! Allocate arrays
    call mem_allocate(this%rptx, this%nreleasepoints, 'RPTX', this%memoryPath)
    call mem_allocate(this%rpty, this%nreleasepoints, 'RPTY', this%memoryPath)
    call mem_allocate(this%rptz, this%nreleasepoints, 'RPTZ', this%memoryPath)
    call mem_allocate(this%rptm, this%nreleasepoints, 'RPTMASS', &
                      this%memoryPath)
    call mem_allocate(this%rptnode, this%nreleasepoints, 'RPTNODER', &
                      this%memoryPath)
    call mem_allocate(this%rptname, lenboundname, this%nreleasepoints, &
                      'RPTNAME', this%memoryPath)
 
    ! Initialize arrays
    do nps = 1, this%nreleasepoints
      this%rptm(nps) = dzero
    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 242 of file prt-prp.f90.

    class(PrtPrpType) :: this
 
    ! Allocate parent's scalars
    call this%BndExtType%allocate_scalars()
 
    ! Allocate scalars for this type
    call mem_allocate(this%localz, 'LOCALZ', this%memoryPath)
    call mem_allocate(this%extend, 'EXTEND', this%memoryPath)
    call mem_allocate(this%offset, 'OFFSET', this%memoryPath)
    call mem_allocate(this%stoptime, 'STOPTIME', this%memoryPath)
    call mem_allocate(this%stoptraveltime, 'STOPTRAVELTIME', this%memoryPath)
    call mem_allocate(this%istopweaksink, 'ISTOPWEAKSINK', this%memoryPath)
    call mem_allocate(this%istopzone, 'ISTOPZONE', this%memoryPath)
    call mem_allocate(this%drape, 'DRAPE', this%memoryPath)
    call mem_allocate(this%idrymeth, 'IDRYMETH', this%memoryPath)
    call mem_allocate(this%nreleasepoints, 'NRELEASEPOINTS', this%memoryPath)
    call mem_allocate(this%nreleasetimes, 'NRELEASETIMES', this%memoryPath)
    call mem_allocate(this%nparticles, 'NPARTICLES', this%memoryPath)
    call mem_allocate(this%itrkout, 'ITRKOUT', this%memoryPath)
    call mem_allocate(this%itrkhdr, 'ITRKHDR', this%memoryPath)
    call mem_allocate(this%itrkcsv, 'ITRKCSV', this%memoryPath)
    call mem_allocate(this%irlstls, 'IRLSTLS', this%memoryPath)
    call mem_allocate(this%frctrn, 'FRCTRN', this%memoryPath)
    call mem_allocate(this%iexmeth, 'IEXMETH', this%memoryPath)
    call mem_allocate(this%ichkmeth, 'ICHKMETH', this%memoryPath)
    call mem_allocate(this%icycwin, 'ICYCWIN', this%memoryPath)
    call mem_allocate(this%extol, 'EXTOL', this%memoryPath)
    call mem_allocate(this%rttol, 'RTTOL', this%memoryPath)
    call mem_allocate(this%rtfreq, 'RTFREQ', this%memoryPath)
 
    ! Set values
    this%localz = .false.
    this%extend = .false.
    this%offset = dzero
    this%stoptime = huge(1d0)
    this%stoptraveltime = huge(1d0)
    this%istopweaksink = 0
    this%istopzone = 0
    this%drape = .false.
    this%idrymeth = 0
    this%nreleasepoints = 0
    this%nreleasetimes = 0
    this%nparticles = 0
    this%itrkout = 0
    this%itrkhdr = 0
    this%itrkcsv = 0
    this%irlstls = 0
    this%frctrn = .false.
    this%iexmeth = 0
    this%ichkmeth = 1
    this%icycwin = 0
    this%extol = default_exit_solve_tolerance
    this%rttol = dsame * dep9
    this%rtfreq = dzero
 

prp_ar()

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

private

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

    ! dummy variables
    class(PrtPrpType), intent(inout) :: this
    ! local variables
    integer(I4B) :: n
 
    call this%obs%obs_ar()
 
    if (this%inamedbound /= 0) then
      do n = 1, this%nreleasepoints
        this%boundname(n) = this%rptname(n)
      end do
    end if
    do n = 1, this%nreleasepoints
      this%nodelist(n) = this%rptnode(n)
    end do

prp_cq_simrate()

subroutine prtprpmodule::prp_cq_simrate	(	class(prtprptype)	this,
		real(dp), dimension(:), intent(in)	hnew,
		real(dp), dimension(:), intent(inout)	flowja,
		integer(i4b), intent(in)	imover
	)

Parameters

[in,out]	flowja	flow between package and model
[in]	imover	flag indicating if the mover package is active

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

    ! modules
    use tdismodule, only: delt
    ! dummy variables
    class(PrtPrpType) :: this
    real(DP), dimension(:), intent(in) :: hnew
    real(DP), dimension(:), intent(inout) :: flowja !< flow between package and model
    integer(I4B), intent(in) :: imover !< flag indicating if the mover package is active
    ! local variables
    integer(I4B) :: i
    integer(I4B) :: node
    integer(I4B) :: idiag
    real(DP) :: rrate
 
    ! If no boundaries, skip flow calculations.
    if (this%nbound <= 0) return
 
    ! Loop through each boundary calculating flow.
    do i = 1, this%nbound
      node = this%nodelist(i)
      rrate = dzero
      ! If cell is no-flow or constant-head, then ignore it.
      if (node > 0) then
        ! Calculate the flow rate into the cell.
        idiag = this%dis%con%ia(node)
        rrate = this%rptm(i) * (done / delt) ! reciprocal of tstp length
        flowja(idiag) = flowja(idiag) + rrate
      end if
 
      ! Save simulated value to simvals array.
      this%simvals(i) = rrate
    end do

prp_create()

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

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

    ! dummy
    class(BndType), pointer :: packobj
    integer(I4B), intent(in) :: id
    integer(I4B), intent(in) :: ibcnum
    integer(I4B), intent(in) :: inunit
    integer(I4B), intent(in) :: iout
    character(len=*), intent(in) :: namemodel
    character(len=*), intent(in) :: pakname
    character(len=*), intent(in) :: input_mempath
    type(PrtFmiType), pointer :: fmi
    ! local
    type(PrtPrpType), pointer :: prpobj
    ! formats
    character(len=*), parameter :: fmtheader = &
      "(1x, /1x, 'PRP PARTICLE RELEASE POINT PACKAGE', &
       &' INPUT READ FROM MEMPATH: ', a, /)"
 
    ! allocate the object and assign values to object variables
    allocate (prpobj)
    packobj => prpobj
 
    ! create name and memory path
    call packobj%set_names(ibcnum, namemodel, pakname, ftype, input_mempath)
    prpobj%text = text
 
    ! allocate scalars
    call prpobj%prp_allocate_scalars()
 
    ! initialize package
    call packobj%pack_initialize()
 
    packobj%inunit = inunit
    packobj%iout = iout
    packobj%id = id
    packobj%ibcnum = ibcnum
    packobj%ncolbnd = 4
    packobj%iscloc = 1
 
    ! store pointer to flow model interface
    prpobj%fmi => fmi
 
    ! if prp is enabled, print a message identifying it
    if (inunit > 0) write (iout, fmtheader) input_mempath

Here is the caller graph for this function:

prp_da()

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

private

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

    class(PrtPrpType) :: this
 
    ! Deallocate parent
    call this%BndExtType%bnd_da()
 
    ! Deallocate scalars
    call mem_deallocate(this%localz)
    call mem_deallocate(this%extend)
    call mem_deallocate(this%offset)
    call mem_deallocate(this%stoptime)
    call mem_deallocate(this%stoptraveltime)
    call mem_deallocate(this%istopweaksink)
    call mem_deallocate(this%istopzone)
    call mem_deallocate(this%drape)
    call mem_deallocate(this%idrymeth)
    call mem_deallocate(this%nreleasepoints)
    call mem_deallocate(this%nreleasetimes)
    call mem_deallocate(this%nparticles)
    call mem_deallocate(this%itrkout)
    call mem_deallocate(this%itrkhdr)
    call mem_deallocate(this%itrkcsv)
    call mem_deallocate(this%irlstls)
    call mem_deallocate(this%frctrn)
    call mem_deallocate(this%iexmeth)
    call mem_deallocate(this%ichkmeth)
    call mem_deallocate(this%icycwin)
    call mem_deallocate(this%extol)
    call mem_deallocate(this%rttol)
    call mem_deallocate(this%rtfreq)
 
    ! Deallocate arrays
    call mem_deallocate(this%rptx)
    call mem_deallocate(this%rpty)
    call mem_deallocate(this%rptz)
    call mem_deallocate(this%rptnode)
    call mem_deallocate(this%rptm)
    call mem_deallocate(this%rptname, 'RPTNAME', this%memoryPath)
 
    ! Deallocate objects
    call this%particles%destroy(this%memoryPath)
    call this%schedule%destroy()
    deallocate (this%particles)
    deallocate (this%schedule)

prp_df_obs()

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

private

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

    ! dummy
    class(PrtPrpType) :: this
    ! local
    integer(I4B) :: indx
    call this%obs%StoreObsType('prp', .true., indx)
    this%obs%obsData(indx)%ProcessIdPtr => defaultobsidprocessor
 
    ! Store obs type and assign procedure pointer
    ! for to-mvr observation type.
    call this%obs%StoreObsType('to-mvr', .true., indx)
    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 879 of file prt-prp.f90.

    ! -- modules
    use memorymanagerextmodule, only: mem_set_value
    use prtprpinputmodule, only: prtprpparamfoundtype
    ! -- dummy variables
    class(PrtPrpType), intent(inout) :: this
    ! -- local variables
    type(PrtPrpParamFoundType) :: found
 
    call mem_set_value(this%nreleasepoints, 'NRELEASEPTS', this%input_mempath, &
                       found%nreleasepts)
    call mem_set_value(this%nreleasetimes, 'NRELEASETIMES', this%input_mempath, &
                       found%nreleasetimes)
 
    write (this%iout, '(1x,a)') 'PROCESSING PARTICLE INPUT DIMENSIONS'
    write (this%iout, '(4x,a,i0)') 'NRELEASEPTS = ', this%nreleasepoints
    write (this%iout, '(4x,a,i0)') 'NRELEASETIMES = ', this%nreleasetimes
    write (this%iout, '(1x,a)') 'END OF PARTICLE INPUT DIMENSIONS'
 
    ! set maxbound and nbound to nreleasepts
    this%maxbound = this%nreleasepoints
    this%nbound = this%nreleasepoints
 
    ! allocate arrays for prp package
    call this%prp_allocate_arrays()
 
    ! read packagedata and releasetimes blocks
    call this%prp_packagedata()
    call this%prp_releasetimes()
    call this%prp_load_releasetimefrequency()

prp_load_releasetimefrequency()

subroutine prtprpmodule::prp_load_releasetimefrequency

(

class(prtprptype), intent(inout)

this

)

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

    ! modules
    use tdismodule, only: totalsimtime
    ! dummy
    class(PrtPrpType), intent(inout) :: this
    ! local
    real(DP), allocatable :: times(:)
 
    ! check if a release time frequency is configured
    if (this%rtfreq <= dzero) return
 
    ! create array of regularly-spaced release times
    times = arange( &
            start=dzero, &
            stop=totalsimtime, &
            step=this%rtfreq)
 
    ! register times with release schedule
    call this%schedule%time_select%extend(times)
 
    ! make sure times strictly increase
    if (.not. this%schedule%time_select%increasing()) then
      errmsg = "Release times must strictly increase"
      call store_error(errmsg)
      call store_error_filename(this%input_fname)
    end if
 
    ! deallocate
    deallocate (times)
 

Here is the call graph for this function:

prp_log_options()

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

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

    ! -- modules
    use prtprpinputmodule, only: prtprpparamfoundtype
    ! -- dummy variables
    class(PrtPrpType), intent(inout) :: this
    type(PrtPrpParamFoundType), intent(in) :: found
    character(len=*), intent(in) :: trackfile
    character(len=*), intent(in) :: trackcsvfile
    ! -- local variables
    ! formats
    character(len=*), parameter :: fmttrkbin = &
      "(4x, 'PARTICLE TRACKS WILL BE SAVED TO BINARY FILE: ', a, /4x, &
    &'OPENED ON UNIT: ', I0)"
    character(len=*), parameter :: fmttrkcsv = &
      "(4x, 'PARTICLE TRACKS WILL BE SAVED TO CSV FILE: ', a, /4x, &
    &'OPENED ON UNIT: ', I0)"
 
    write (this%iout, '(1x,a)') 'PROCESSING PARTICLE INPUT DIMENSIONS'
 
    if (found%frctrn) then
      write (this%iout, '(4x,a)') &
        'IF DISV, TRACKING WILL USE THE TERNARY METHOD REGARDLESS OF CELL TYPE'
    end if
 
    if (found%trackfile) then
      write (this%iout, fmttrkbin) trim(adjustl(trackfile)), this%itrkout
    end if
 
    if (found%trackcsvfile) then
      write (this%iout, fmttrkcsv) trim(adjustl(trackcsvfile)), this%itrkcsv
    end if
 
    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 671 of file prt-prp.f90.

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

prp_options()

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

private

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

    ! -- modules
    use constantsmodule, only: lenvarname, dzero, mnormal
    use memorymanagerextmodule, only: mem_set_value
    use openspecmodule, only: access, form
    use inputoutputmodule, only: getunit, openfile
    use prtprpinputmodule, only: prtprpparamfoundtype
    ! -- dummy variables
    class(PrtPrpType), intent(inout) :: this
    ! -- local variables
    character(len=LENVARNAME), dimension(3) :: drytrack_method = &
      &[character(len=LENVARNAME) :: 'DROP', 'STOP', 'STAY']
    character(len=lenvarname), dimension(2) :: coorcheck_method = &
      &[character(len=LENVARNAME) :: 'NONE', 'EAGER']
    character(len=LINELENGTH) :: trackfile, trackcsvfile, fname
    type(PrtPrpParamFoundType) :: found
    character(len=*), parameter :: fmtextolwrn = &
      "('WARNING: EXIT_SOLVE_TOLERANCE is set to ',g10.3,' &
      &which is much greater than the default value of ',g10.3,'. &
      &The tolerance that strikes the best balance between accuracy &
      &and runtime is problem-dependent. Since the variable being &
      &solved varies from 0 to 1, tolerance values much less than 1 &
      &typically give the best results.')"
 
    ! -- source base class options
    call this%BndExtType%source_options()
 
    ! -- update defaults from input context
    call mem_set_value(this%stoptime, 'STOPTIME', this%input_mempath, &
                       found%stoptime)
    call mem_set_value(this%stoptraveltime, 'STOPTRAVELTIME', &
                       this%input_mempath, found%stoptraveltime)
    call mem_set_value(this%istopweaksink, 'ISTOPWEAKSINK', this%input_mempath, &
                       found%istopweaksink)
    call mem_set_value(this%istopzone, 'ISTOPZONE', this%input_mempath, &
                       found%istopzone)
    call mem_set_value(this%drape, 'DRAPE', this%input_mempath, &
                       found%drape)
    call mem_set_value(this%idrymeth, 'IDRYMETH', this%input_mempath, &
                       drytrack_method, found%idrymeth)
    call mem_set_value(trackfile, 'TRACKFILE', this%input_mempath, &
                       found%trackfile)
    call mem_set_value(trackcsvfile, 'TRACKCSVFILE', this%input_mempath, &
                       found%trackcsvfile)
    call mem_set_value(this%localz, 'LOCALZ', this%input_mempath, &
                       found%localz)
    call mem_set_value(this%extend, 'EXTEND', this%input_mempath, &
                       found%extend)
    call mem_set_value(this%extol, 'EXTOL', this%input_mempath, &
                       found%extol)
    call mem_set_value(this%rttol, 'RTTOL', this%input_mempath, &
                       found%rttol)
    call mem_set_value(this%rtfreq, 'RTFREQ', this%input_mempath, &
                       found%rtfreq)
    call mem_set_value(this%frctrn, 'FRCTRN', this%input_mempath, &
                       found%frctrn)
    call mem_set_value(this%iexmeth, 'IEXMETH', this%input_mempath, &
                       found%iexmeth)
    call mem_set_value(this%ichkmeth, 'ICHKMETH', this%input_mempath, &
                       coorcheck_method, found%ichkmeth)
    call mem_set_value(this%icycwin, 'ICYCWIN', this%input_mempath, found%icycwin)
 
    ! update internal state and validate input
    if (found%idrymeth) then
      if (this%idrymeth == 0) then
        write (errmsg, '(a)') 'Unsupported dry tracking method. &
          &DRY_TRACKING_METHOD must be "DROP", "STOP", or "STAY"'
        call store_error(errmsg)
      else
        ! adjust for method zero indexing
        this%idrymeth = this%idrymeth - 1
      end if
    end if
 
    if (found%extol) then
      if (this%extol <= dzero) &
        call store_error('EXIT_SOLVE_TOLERANCE MUST BE POSITIVE')
      if (this%extol > dem2) then
        write (warnmsg, fmt=fmtextolwrn) &
          this%extol, default_exit_solve_tolerance
        call store_warning(warnmsg)
      end if
    end if
 
    if (found%rttol) then
      if (this%rttol <= dzero) &
        call store_error('RELEASE_TIME_TOLERANCE MUST BE POSITIVE')
    end if
 
    if (found%rtfreq) then
      if (this%rtfreq <= dzero) &
        call store_error('RELEASE_TIME_FREQUENCY MUST BE POSITIVE')
    end if
 
    if (found%iexmeth) then
      if (.not. (this%iexmeth /= 1 .or. this%iexmeth /= 2)) &
        call store_error('DEV_EXIT_SOLVE_METHOD MUST BE &
          &1 (BRENT) OR 2 (CHANDRUPATLA)')
    end if
 
    if (found%ichkmeth) then
      if (this%ichkmeth == 0) then
        write (errmsg, '(a)') 'Unsupported coordinate check method. &
          &COORDINATE_CHECK_METHOD must be "NONE" or "EAGER"'
        call store_error(errmsg)
      else
        ! adjust for method zero based indexing
        this%ichkmeth = this%ichkmeth - 1
      end if
    end if
 
    if (found%icycwin) then
      if (this%icycwin < 0) &
        call store_error('CYCLE_DETECTION_WINDOW MUST BE NON-NEGATIVE')
    end if
 
    ! fileout options
    if (found%trackfile) then
      this%itrkout = getunit()
      call openfile(this%itrkout, this%iout, trackfile, 'DATA(BINARY)', &
                    form, access, filstat_opt='REPLACE', &
                    mode_opt=mnormal)
      ! open and write ascii header spec file
      this%itrkhdr = getunit()
      fname = trim(trackfile)//'.hdr'
      call openfile(this%itrkhdr, this%iout, fname, 'CSV', &
                    filstat_opt='REPLACE', mode_opt=mnormal)
      write (this%itrkhdr, '(a,/,a)') trackheader, trackdtypes
    end if
 
    if (found%trackcsvfile) then
      this%itrkcsv = getunit()
      call openfile(this%itrkcsv, this%iout, trackcsvfile, 'CSV', &
                    filstat_opt='REPLACE')
      write (this%itrkcsv, '(a)') trackheader
    end if
 
    ! terminate if any errors were detected
    if (count_errors() > 0) then
      call store_error_filename(this%input_fname)
    end if
 
    ! log found options
    call this%prp_log_options(found, trackfile, trackcsvfile)
 
    ! Create release schedule now that we know
    ! the coincident release time tolerance
    this%schedule => create_release_schedule(tolerance=this%rttol)

Here is the call graph for this function:

prp_packagedata()

subroutine prtprpmodule::prp_packagedata

(

class(prtprptype), intent(inout)

this

)

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

    use memorymanagermodule, only: mem_setptr
    use geomutilmodule, only: get_node
    use characterstringmodule, only: characterstringtype
    ! dummy
    class(PrtPrpType), intent(inout) :: this
    ! local
    integer(I4B), dimension(:), pointer, contiguous :: irptno
    integer(I4B), dimension(:, :), pointer, contiguous :: cellids
    real(DP), dimension(:), pointer, contiguous :: xrpts, yrpts, zrpts
    type(CharacterStringType), dimension(:), pointer, &
      contiguous :: boundnames
    character(len=LENBOUNDNAME) :: bndName, bndNameTemp
    character(len=9) :: cno
    character(len=20) :: cellidstr
    integer(I4B), dimension(:), allocatable :: nboundchk
    integer(I4B), dimension(:), pointer :: cellid
    integer(I4B) :: n, noder, nodeu, rptno
 
    ! set input context pointers
    call mem_setptr(irptno, 'IRPTNO', this%input_mempath)
    call mem_setptr(cellids, 'CELLID', this%input_mempath)
    call mem_setptr(xrpts, 'XRPT', this%input_mempath)
    call mem_setptr(yrpts, 'YRPT', this%input_mempath)
    call mem_setptr(zrpts, 'ZRPT', this%input_mempath)
    call mem_setptr(boundnames, 'BOUNDNAME', this%input_mempath)
 
    ! allocate and initialize temporary variables
    allocate (nboundchk(this%nreleasepoints))
    do n = 1, this%nreleasepoints
      nboundchk(n) = 0
    end do
 
    write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%packName)) &
      //' PACKAGEDATA'
 
    do n = 1, size(irptno)
 
      rptno = irptno(n)
 
      if (rptno < 1 .or. rptno > this%nreleasepoints) then
        write (errmsg, '(a,i0,a,i0,a)') &
          'Expected ', this%nreleasepoints, ' release points. &
          &Points must be numbered from 1 to ', this%nreleasepoints, '.'
        call store_error(errmsg)
        cycle
      end if
 
      ! increment nboundchk
      nboundchk(rptno) = nboundchk(rptno) + 1
 
      ! set cellid
      cellid => cellids(:, n)
 
      ! set node user
      if (this%dis%ndim == 1) then
        nodeu = cellid(1)
      elseif (this%dis%ndim == 2) then
        nodeu = get_node(cellid(1), 1, cellid(2), &
                         this%dis%mshape(1), 1, &
                         this%dis%mshape(2))
      else
        nodeu = get_node(cellid(1), cellid(2), cellid(3), &
                         this%dis%mshape(1), &
                         this%dis%mshape(2), &
                         this%dis%mshape(3))
      end if
 
      ! set noder
      noder = this%dis%get_nodenumber(nodeu, 1)
      if (noder <= 0) then
        call this%dis%nodeu_to_string(nodeu, cellidstr)
        write (errmsg, '(a)') &
          'Particle release point configured for nonexistent cell: '// &
          trim(adjustl(cellidstr))//'. This cell has IDOMAIN <= 0 and '&
          &'therefore does not exist in the model grid.'
        call store_error(errmsg)
        cycle
      else
        this%rptnode(rptno) = noder
      end if
 
      if (this%localz .and. (zrpts(n) < 0 .or. zrpts(n) > 1)) then
        call store_error('Local z coordinate must fall in the interval [0, 1]')
        cycle
      end if
 
      ! set coordinates
      this%rptx(rptno) = xrpts(n)
      this%rpty(rptno) = yrpts(n)
      this%rptz(rptno) = zrpts(n)
 
      ! set default boundname
      write (cno, '(i9.9)') rptno
      bndname = 'PRP'//cno
 
      ! read boundnames from file, if provided
      if (this%inamedbound /= 0) then
        bndnametemp = boundnames(n)
        if (bndnametemp /= '') bndname = bndnametemp
      else
        bndname = ''
      end if
 
      ! set boundname
      this%rptname(rptno) = bndname
    end do
 
    write (this%iout, '(1x,a)') &
      'END OF '//trim(adjustl(this%packName))//' PACKAGEDATA'
 
    ! check for duplicate or missing particle release points
    do n = 1, this%nreleasepoints
      if (nboundchk(n) == 0) then
        write (errmsg, '(a,a,1x,i0,a)') 'No data specified for particle ', &
          'release point', n, '.'
        call store_error(errmsg)
      else if (nboundchk(n) > 1) then
        write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
          'Data for particle release point', n, 'specified', nboundchk(n), &
          'times.'
        call store_error(errmsg)
      end if
    end do
 
    ! terminate if any errors were detected
    if (count_errors() > 0) then
      call store_error_filename(this%input_fname)
    end if
 
    ! cleanup
    deallocate (nboundchk)

Here is the call graph for this function:

prp_releasetimes()

subroutine prtprpmodule::prp_releasetimes

(

class(prtprptype), intent(inout)

this

)

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

    use memorymanagermodule, only: mem_setptr, get_isize
    ! dummy
    class(PrtPrpType), intent(inout) :: this
    ! local
    real(DP), dimension(:), pointer, contiguous :: time
    integer(I4B) :: n, isize
    real(DP), allocatable :: times(:)
 
    if (this%nreleasetimes <= 0) return
 
    ! allocate times array
    allocate (times(this%nreleasetimes))
 
    ! check if input array was read
    call get_isize('TIME', this%input_mempath, isize)
 
    if (isize <= 0) then
      errmsg = "RELEASTIMES block expected when &
        &NRELEASETIMES dimension is non-zero."
      call store_error(errmsg)
      call store_error_filename(this%input_fname)
    end if
 
    ! set input context pointer
    call mem_setptr(time, 'TIME', this%input_mempath)
 
    ! set input data
    do n = 1, size(time)
      times(n) = time(n)
    end do
 
    ! register times with the release schedule
    call this%schedule%time_select%extend(times)
 
    ! make sure times strictly increase
    if (.not. this%schedule%time_select%increasing()) then
      errmsg = "RELEASTIMES block entries must strictly increase."
      call store_error(errmsg)
      call store_error_filename(this%input_fname)
    end if
 
    ! deallocate
    deallocate (times)

Here is the call graph for this function:

prp_rp()

subroutine prtprpmodule::prp_rp

(

class(prtprptype), intent(inout)

this

)

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

    ! modules
    use tdismodule, only: kper, nper
    use memorymanagermodule, only: mem_setptr
    use characterstringmodule, only: characterstringtype
    ! dummy variables
    class(PrtPrpType), intent(inout) :: this
    ! local variables
    type(CharacterStringType), dimension(:), contiguous, &
      pointer :: settings
    integer(I4B), pointer :: iper, ionper, nlist
    character(len=LINELENGTH), allocatable :: lines(:)
    integer(I4B) :: n
 
    ! set pointer to last and next period loaded
    call mem_setptr(iper, 'IPER', this%input_mempath)
    call mem_setptr(ionper, 'IONPER', this%input_mempath)
 
    if (kper == 1 .and. &
        (iper == 0) .and. &
        (ionper > nper) .and. &
        size(this%schedule%time_select%times) == 0) then
      ! If the user hasn't provided any release settings (neither
      ! explicit release times, release time frequency, or period
      ! block release settings), default to a single release at the
      ! start of the first period's first time step.
      allocate (lines(1))
      lines(1) = "FIRST"
      call this%schedule%advance(lines=lines)
      deallocate (lines)
      return
    else if (iper /= kper) then
      return
    end if
 
    ! set input context pointers
    call mem_setptr(nlist, 'NBOUND', this%input_mempath)
    call mem_setptr(settings, 'SETTING', this%input_mempath)
 
    ! allocate and set input
    allocate (lines(nlist))
    do n = 1, nlist
      lines(n) = settings(n)
    end do
 
    ! update schedule
    if (size(lines) > 0) &
      call this%schedule%advance(lines=lines)
 
    ! cleanup
    deallocate (lines)

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

    class(PrtPrpType) :: this
    integer(I4B), dimension(:), pointer, contiguous :: ibound
    integer(I4B), dimension(:), pointer, contiguous :: izone
 
    this%ibound => ibound
    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]	this	this instance
[in]	ip	particle index
[in]	trelease	release time

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

    ! dummy
    class(PrtPrpType), intent(inout) :: this !< this instance
    integer(I4B), intent(in) :: ip !< particle index
    real(DP), intent(in) :: trelease !< release time
    ! local
    integer(I4B) :: np
    type(ParticleType), pointer :: particle
 
    call this%initialize_particle(particle, ip, trelease)
    np = this%nparticles + 1
    this%nparticles = np
    call this%particles%put(particle, np)
    deallocate (particle)
    this%rptm(ip) = this%rptm(ip) + done ! TODO configurable mass
 

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]	this	this instance
[in]	ic	cell index
[in]	z	release point

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

    class(PrtPrpType), intent(inout) :: this !< this instance
    integer(I4B), intent(in) :: ic !< cell index
    real(DP), intent(in) :: x, y, z !< release point
    ! local
    real(DP), allocatable :: polyverts(:, :)
 
    call this%fmi%dis%get_polyverts(ic, polyverts)
    if (.not. point_in_polygon(x, y, polyverts)) then
      write (errmsg, '(a,g0,a,g0,a,i0)') &
        'Error: release point (x=', x, ', y=', y, ') is not in cell ', &
        this%dis%get_nodeuser(ic)
      call store_error(errmsg, terminate=.false.)
      call store_error_filename(this%input_fname)
    end if
    if (z > maxval(this%dis%top)) then
      write (errmsg, '(a,g0,a,g0,a,i0)') &
        'Error: release point (z=', z, ') is above grid top ', &
        maxval(this%dis%top)
      call store_error(errmsg, terminate=.false.)
      call store_error_filename(this%input_fname)
    else if (z < minval(this%dis%bot)) then
      write (errmsg, '(a,g0,a,g0,a,i0)') &
        'Error: release point (z=', z, ') is below grid bottom ', &
        minval(this%dis%bot)
      call store_error(errmsg, terminate=.false.)
      call store_error_filename(this%input_fname)
    end if
    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.

  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.

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

text

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

private

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

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