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, fmi)
	Create a new particle release point package. More...
subroutine	prp_da (this)
	Deallocate memory. More...
subroutine	prp_set_pointers (this, ibound, izone, trackctl)
	@ 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	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, option, found)
	Set options specific to PrtPrpType. More...
subroutine	prp_read_dimensions (this)
	Read package dimensions. More...
subroutine	prp_read_packagedata (this)
	Load package data (release points). More...
subroutine	prp_read_releasetimes (this)
	Load explicitly specified release times. More...
subroutine	prp_load_releasetimefrequency (this)
	Load regularly spaced release times if configured. More...

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

Function/Subroutine Documentation

define_listlabel()

subroutine prtprpmodule::define_listlabel

(

class(prtprptype), intent(inout)

this

)

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

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

log_release()

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

private

Parameters

[in,out]

this

prp

Definition at line 352 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 309 of file prt-prp.f90.

    class(PrtPrpType) :: this
    integer(I4B) :: ip, it
 
    ! 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()
        call this%release(ip, this%schedule%times(it))
      end do
    end do

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 203 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
 
    ! Allocate particle store, starting with the number
    ! of release points (arrays resized if/when needed)
    call allocate_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 235 of file prt-prp.f90.

    class(PrtPrpType) :: this
 
    ! Allocate parent's scalars
    call this%BndType%allocate_scalars()
 
    ! Allocate scalars for this type
    call mem_allocate(this%ilocalz, 'ILOCALZ', this%memoryPath)
    call mem_allocate(this%iextend, 'IEXTEND', 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%idrape, 'IDRAPE', 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%ifrctrn, 'IFRCTRN', this%memoryPath)
    call mem_allocate(this%iexmeth, 'IEXMETH', 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)
    call mem_allocate(this%foundtol, 'FOUNDTOL', this%memoryPath)
 
    ! Set values
    this%ilocalz = 0
    this%iextend = 0
    this%offset = dzero
    this%stoptime = huge(1d0)
    this%stoptraveltime = huge(1d0)
    this%istopweaksink = 0
    this%istopzone = 0
    this%idrape = 0
    this%nreleasepoints = 0
    this%nreleasetimes = 0
    this%nparticles = 0
    this%itrkout = 0
    this%itrkhdr = 0
    this%itrkcsv = 0
    this%irlstls = 0
    this%ifrctrn = 0
    this%iexmeth = 0
    this%extol = dzero
    this%rttol = dsame * dep9
    this%rtfreq = dzero
    this%foundtol = .false.
 

prp_ar()

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

private

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

    ! dummy variables
    class(PrtPrpType), intent(inout) :: this
    ! local variables
    integer(I4B) :: n
 
    call this%obs%obs_ar()
    call this%BndType%allocate_arrays()
    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 597 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,
		type(prtfmitype), pointer	fmi
	)

Definition at line 99 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
    type(PrtFmiType), pointer :: fmi
    ! local
    type(PrtPrpType), pointer :: prpobj
    ! formats
    character(len=*), parameter :: fmtheader = &
      "(1x, /1x, 'PRP PARTICLE RELEASE POINT PACKAGE', &
       &' INPUT READ FROM UNIT ', i0, /)"
 
    ! 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)
    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) inunit

Here is the caller graph for this function:

prp_da()

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

private

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

    class(PrtPrpType) :: this
 
    ! Deallocate parent
    call this%BndType%bnd_da()
 
    ! Deallocate scalars
    call mem_deallocate(this%ilocalz)
    call mem_deallocate(this%iextend)
    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%idrape)
    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%ifrctrn)
    call mem_deallocate(this%iexmeth)
    call mem_deallocate(this%extol)
    call mem_deallocate(this%rttol)
    call mem_deallocate(this%rtfreq)
    call mem_deallocate(this%foundtol)
 
    ! 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 particle store, release time and release step selections
    call this%particles%deallocate(this%memoryPath)
    call this%schedule%deallocate()
    deallocate (this%particles)
    deallocate (this%schedule)

prp_df_obs()

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

private

Definition at line 643 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_load_releasetimefrequency()

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

private

Definition at line 1037 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 this%parser%StoreErrorUnit(terminate=.true.)
    end if
 
    ! deallocate
    deallocate (times)
 

Here is the call graph for this function:

prp_obs_supported()

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

private

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

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

prp_options()

subroutine prtprpmodule::prp_options ( class(prtprptype), intent(inout)  this, character(len=*), intent(inout)  option, logical(lgp), intent(inout)  found  )

private

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

    use openspecmodule, only: access, form
    use constantsmodule, only: maxcharlen, dzero
    use inputoutputmodule, only: urword, getunit, openfile
    use trackfilemodule, only: trackheader, trackdtypes
    ! dummy
    class(PrtPrpType), intent(inout) :: this
    character(len=*), intent(inout) :: option
    logical(LGP), intent(inout) :: found
    ! locals
    character(len=MAXCHARLEN) :: fname
    character(len=MAXCHARLEN) :: keyword
    ! 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)"
 
    select case (option)
    case ('STOPTIME')
      this%stoptime = this%parser%GetDouble()
      found = .true.
    case ('STOPTRAVELTIME')
      this%stoptraveltime = this%parser%GetDouble()
      found = .true.
    case ('STOP_AT_WEAK_SINK')
      this%istopweaksink = 1
      found = .true.
    case ('ISTOPZONE')
      this%istopzone = this%parser%GetInteger()
      found = .true.
    case ('DRAPE')
      this%idrape = 1
      found = .true.
    case ('TRACK')
      call this%parser%GetStringCaps(keyword)
      if (keyword == 'FILEOUT') then
        ! parse filename
        call this%parser%GetString(fname)
        ! open binary output file
        this%itrkout = getunit()
        call openfile(this%itrkout, this%iout, fname, 'DATA(BINARY)', &
                      form, access, filstat_opt='REPLACE', &
                      mode_opt=mnormal)
        write (this%iout, fmttrkbin) trim(adjustl(fname)), this%itrkout
        ! open and write ascii header spec file
        this%itrkhdr = getunit()
        fname = trim(fname)//'.hdr'
        call openfile(this%itrkhdr, this%iout, fname, 'CSV', &
                      filstat_opt='REPLACE', mode_opt=mnormal)
        write (this%itrkhdr, '(a,/,a)') trackheader, trackdtypes
      else
        call store_error('OPTIONAL TRACK KEYWORD MUST BE '// &
                         'FOLLOWED BY FILEOUT')
      end if
      found = .true.
    case ('TRACKCSV')
      call this%parser%GetStringCaps(keyword)
      if (keyword == 'FILEOUT') then
        ! parse filename
        call this%parser%GetString(fname)
        ! open CSV output file and write headers
        this%itrkcsv = getunit()
        call openfile(this%itrkcsv, this%iout, fname, 'CSV', &
                      filstat_opt='REPLACE')
        write (this%iout, fmttrkcsv) trim(adjustl(fname)), this%itrkcsv
        write (this%itrkcsv, '(a)') trackheader
      else
        call store_error('OPTIONAL TRACKCSV KEYWORD MUST BE &
          &FOLLOWED BY FILEOUT')
      end if
      found = .true.
    case ('LOCAL_Z')
      this%ilocalz = 1
      found = .true.
    case ('EXTEND_TRACKING')
      this%iextend = 1
      found = .true.
    case ('EXIT_SOLVE_TOLERANCE')
      this%extol = this%parser%GetDouble()
      if (this%extol <= dzero) &
        call store_error('EXIT_SOLVE_TOLERANCE MUST BE POSITIVE')
      found = .true.
      this%foundtol = .true.
    case ('RELEASE_TIME_TOLERANCE')
      this%rttol = this%parser%GetDouble()
      if (this%rttol <= dzero) &
        call store_error('RELEASE_TIME_TOLERANCE MUST BE POSITIVE')
      found = .true.
    case ('RELEASE_TIME_FREQUENCY')
      this%rtfreq = this%parser%GetDouble()
      if (this%rtfreq <= dzero) &
        call store_error('RELEASE_TIME_FREQUENCY MUST BE POSITIVE')
      found = .true.
    case ('DEV_FORCETERNARY')
      call this%parser%DevOpt()
      this%ifrctrn = 1
      write (this%iout, '(4x,a)') &
        'TRACKING WILL BE DONE USING THE TERNARY METHOD REGARDLESS OF CELL TYPE'
      found = .true.
    case ('DEV_EXIT_SOLVE_METHOD')
      call this%parser%DevOpt()
      this%iexmeth = this%parser%GetInteger()
      if (.not. (this%iexmeth /= 1 .or. this%iexmeth /= 2)) &
        call store_error('DEV_EXIT_SOLVE_METHOD MUST BE &
          &1 (BRENT) OR 2 (CHANDRUPATLA)')
      found = .true.
    case default
      found = .false.
    end select
 
    ! Catch unrecognized options
    if (.not. found) then
      errmsg = "UNKNOWN PRP OPTION '"//trim(keyword)//"'."
      call store_error(errmsg)
      call this%parser%StoreErrorUnit()
    end if
 
    ! Create release schedule now that we know
    ! the coincident release time tolerance
    this%schedule => create_release_schedule(tol=this%rttol)
 

Here is the call graph for this function:

prp_read_dimensions()

subroutine prtprpmodule::prp_read_dimensions

(

class(prtprptype), intent(inout)

this

)

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

    ! dummy
    class(PrtPrpType), intent(inout) :: this
    ! local
    character(len=LINELENGTH) :: errmsg, keyword
    integer(I4B) :: ierr
    logical :: isfound, endOfBlock
 
    if (.not. this%foundtol) &
      call store_error('EXIT_SOLVE_TOLERANCE MISSING, VALUE REQUIRED')
 
    ! get dimension block
    call this%parser%GetBlock('DIMENSIONS', isfound, ierr, &
                              supportopenclose=.true.)
 
    ! parse dimension block if detected
    if (isfound) then
      write (this%iout, '(1x,a)') 'PROCESSING PARTICLE INPUT DIMENSIONS'
      do
        call this%parser%GetNextLine(endofblock)
        if (endofblock) exit
        call this%parser%GetStringCaps(keyword)
        select case (keyword)
        case ('NRELEASEPTS')
          this%nreleasepoints = this%parser%GetInteger()
        case ('NRELEASETIMES')
          this%nreleasetimes = this%parser%GetInteger()
        case default
          write (errmsg, &
                 '(4x,a,a)') '****ERROR. UNKNOWN PARTICLE INPUT DIMENSION: ', &
            trim(keyword)
          call store_error(errmsg)
          call this%parser%StoreErrorUnit()
        end select
      end do
      write (this%iout, '(1x,a)') 'END OF PARTICLE INPUT DIMENSIONS'
    else
      call store_error('ERROR.  REQUIRED DIMENSIONS BLOCK NOT FOUND.')
    end if
 
    ! 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_read_packagedata()
    call this%prp_read_releasetimes()
    call this%prp_load_releasetimefrequency()

Here is the call graph for this function:

prp_read_packagedata()

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

private

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

    ! dummy
    class(PrtPrpType), intent(inout) :: this
    ! local
    character(len=LINELENGTH) :: cellid
    character(len=LENBOUNDNAME) :: bndName, bndNameTemp
    character(len=9) :: cno
    logical :: isfound
    logical :: endOfBlock
    integer(I4B) :: ival
    integer(I4B) :: n
    integer(I4B) :: ierr
    character(len=LENBOUNDNAME), dimension(:), allocatable :: nametxt
    integer(I4B), dimension(:), allocatable :: nboundchk
    integer(I4B), dimension(:), allocatable :: noder
    real(DP), dimension(:), allocatable :: x
    real(DP), dimension(:), allocatable :: y
    real(DP), dimension(:), allocatable :: z
    real(DP), dimension(:), allocatable :: tstop
    ! format
    character(len=*), parameter :: fmttend = &
      "('end time (', G0, ') must be greater than or equal to the              &
     &begin time (', G0, ').')"
 
    ! allocate temporary variables
    allocate (noder(this%nreleasepoints))
    allocate (x(this%nreleasepoints))
    allocate (y(this%nreleasepoints))
    allocate (z(this%nreleasepoints))
    allocate (tstop(this%nreleasepoints))
    allocate (nametxt(this%nreleasepoints))
    allocate (nboundchk(this%nreleasepoints))
 
    ! initialize temporary variables
    do n = 1, this%nreleasepoints
      nboundchk(n) = 0
    end do
 
    ! read particle release point data
    ! get particle release points block
    call this%parser%GetBlock('PACKAGEDATA', isfound, ierr, &
                              supportopenclose=.true.)
 
    ! parse block if detected
    if (isfound) then
      write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%packName)) &
        //' PACKAGEDATA'
      do
        call this%parser%GetNextLine(endofblock)
        if (endofblock) exit
        ival = this%parser%GetInteger()
        n = ival
 
        if (n < 1 .or. n > this%nreleasepoints) then
          write (errmsg, '(a,1x,i0,a)') &
            'Release point number must be greater than 0 and less than ', &
            'or equal to', this%nreleasepoints, '.'
          call store_error(errmsg)
          cycle
        end if
 
        ! increment nboundchk
        nboundchk(n) = nboundchk(n) + 1
 
        ! node number
        call this%parser%GetCellid(this%dis%ndim, cellid)
        noder(n) = this%dis%noder_from_cellid(cellid, this%inunit, this%iout)
 
        ! x, y, z coordinates
        x(n) = this%parser%GetDouble()
        y(n) = this%parser%GetDouble()
        z(n) = this%parser%GetDouble()
 
        if (this%ilocalz > 0 .and. (z(n) < 0 .or. z(n) > 1)) then
          call store_error('Local z coordinate must fall in the interval [0, 1]')
          cycle
        end if
 
        ! set default boundname
        write (cno, '(i9.9)') n
        bndname = 'PRP'//cno
 
        ! read boundnames from file, if provided
        if (this%inamedbound /= 0) then
          call this%parser%GetStringCaps(bndnametemp)
          if (bndnametemp /= '') &
            bndname = bndnametemp
        else
          bndname = ''
        end if
 
        ! store temp boundnames
        nametxt(n) = 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
    else
      call store_error('Required packagedata block not found.')
    end if
 
    ! terminate if any errors were detected
    if (count_errors() > 0) then
      call this%parser%StoreErrorUnit()
    end if
 
    ! fill particle release point data with data stored in temporary local arrays
    do n = 1, this%nreleasepoints
      this%rptnode(n) = noder(n)
      this%rptx(n) = x(n)
      this%rpty(n) = y(n)
      this%rptz(n) = z(n)
      this%rptname(n) = nametxt(n)
    end do
 
    ! deallocate local storage
    deallocate (noder)
    deallocate (x)
    deallocate (y)
    deallocate (z)
    deallocate (tstop)
    deallocate (nametxt)
    deallocate (nboundchk)

Here is the call graph for this function:

prp_read_releasetimes()

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

private

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

    ! dummy
    class(PrtPrpType), intent(inout) :: this
    ! local
    integer(I4B) :: i, ierr
    logical(LGP) :: eob, found, success
    real(DP) :: t
    real(DP), allocatable :: times(:)
 
    ! get releasetimes block
    call this%parser%GetBlock('RELEASETIMES', found, ierr, &
                              supportopenclose=.true., &
                              blockrequired=.false.)
 
    ! raise an error if releasetimes has a dimension
    ! but no block was found, otherwise return early
    if (.not. found) then
      if (this%nreleasetimes <= 0) return
      write (errmsg, '(a, i0)') &
        "Expected RELEASETIMES with length ", this%nreleasetimes
      call store_error(errmsg)
      call this%parser%StoreErrorUnit(terminate=.true.)
    end if
 
    ! allocate times array
    allocate (times(this%nreleasetimes))
 
    ! read times from the block
    write (this%iout, '(/1x,a)') &
      'PROCESSING '//trim(adjustl(this%text))//' RELEASETIMES'
    do i = 1, this%nreleasetimes
      call this%parser%GetNextLine(eob)
      if (eob) exit
      call this%parser%TryGetDouble(t, success)
      if (.not. success) then
        errmsg = "Failed to read double precision value"
        call store_error(errmsg)
        call this%parser%StoreErrorUnit(terminate=.true.)
      end if
      times(i) = t
    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 = "Release times must strictly increase"
      call store_error(errmsg)
      call this%parser%StoreErrorUnit(terminate=.true.)
    end if
 
    ! deallocate
    deallocate (times)
 

Here is the call graph for this function:

prp_rp()

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

private

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

    ! modules
    use tdismodule, only: kper, nper
    use inputoutputmodule, only: urword
    ! dummy variables
    class(PrtPrpType), intent(inout) :: this
    ! local variables
    integer(I4B) :: ierr
    logical(LGP) :: is_found
    logical(LGP) :: end_of_block
    logical(LGP) :: no_blocks
    character(len=LINELENGTH) :: line
    character(len=LINELENGTH), allocatable :: lines(:)
    ! formats
    character(len=*), parameter :: fmtblkerr = &
                      "('Looking for BEGIN PERIOD iper.  &
                      &Found ', a, ' instead.')"
    character(len=*), parameter :: fmt_steps = &
                                   "(6x,'TIME STEP(S) ',50(I0,' '))" ! 50 limit is similar to STEPS in OC
    character(len=*), parameter :: fmt_freq = &
                                   "(6x,'EVERY ',I0,' TIME STEP(S)')"
    character(len=*), parameter :: fmt_fracs = &
                                   "(6x,50(f10.3,' '))"
 
    ! Set ionper to the stress period number for which a new block of data
    ! will be read.
    if (this%inunit == 0) return
 
    ! get stress period data
    no_blocks = .false.
    if (this%ionper < kper) then
      ! get period block
      call this%parser%GetBlock('PERIOD', is_found, ierr, &
                                supportopenclose=.true., &
                                blockrequired=.false.)
      if (is_found) then
        ! read ionper and check for increasing period numbers
        call this%read_check_ionper()
      else
        ! PERIOD block not found
        if (ierr < 0) then
          if (kper == 1) then
            ! End of file found; no period data for the simulation.
            no_blocks = .true.
          else
            ! End of file found; no more period data.
            this%ionper = nper + 1
          end if
        else
          ! Found invalid block
          call this%parser%GetCurrentLine(line)
          write (errmsg, fmtblkerr) adjustl(trim(line))
          call store_error(errmsg, terminate=.true.)
        end if
      end if
    end if
 
    ! 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.
    if (no_blocks .and. &
        kper == 1 .and. &
        size(this%schedule%time_select%times) == 0) then
      allocate (lines(1))
      line = "FIRST"
      lines(1) = line
      call this%schedule%advance(lines=lines)
    else if (this%ionper == kper) then
      ! If the current stress period matches the
      ! block we are reading, parse the setting
      ! and register it with the schedule.
      allocate (lines(0))
      recordloop: do
        call this%parser%GetNextLine(end_of_block)
        if (end_of_block) exit recordloop
        call this%parser%GetCurrentLine(line)
        call expandarray(lines)
        lines(size(lines)) = line
      end do recordloop
      if (size(lines) > 0) &
        call this%schedule%advance(lines=lines)
      deallocate (lines)
    end if
 

Here is the call graph for this function:

prp_set_pointers()

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

private

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

    class(PrtPrpType) :: this
    integer(I4B), dimension(:), pointer, contiguous :: ibound
    integer(I4B), dimension(:), pointer, contiguous :: izone
    type(TrackControlType), pointer :: trackctl
 
    this%ibound => ibound
    this%rptzone => izone
    this%trackctl => trackctl

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 410 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) :: irow, icol, ilay, icpl
    integer(I4B) :: ic, icu
    integer(I4B) :: np
    real(DP) :: x, y, z
    real(DP) :: top, bot, hds
    type(ParticleType), pointer :: particle
 
    ! Increment cumulative particle count
    np = this%nparticles + 1
    this%nparticles = np
 
    ! Get reduced and user node numbers
    ic = this%rptnode(ip)
    icu = this%dis%get_nodeuser(ic)
 
    ! Load coordinates and transform if needed
    x = this%rptx(ip)
    y = this%rpty(ip)
    if (this%ilocalz > 0) then
      top = this%fmi%dis%top(ic)
      bot = this%fmi%dis%bot(ic)
      hds = this%fmi%gwfhead(ic)
      z = bot + this%rptz(ip) * (hds - bot)
    else
      z = this%rptz(ip)
    end if
 
    ! Make sure release point is in the grid/cell
    call this%validate_release_point(ic, x, y, z)
 
    ! Initialize the particle
    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%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
    ! Handle inactive cells
    if (this%ibound(ic) == 0) then
      ! If drape option activated, release in highest active
      ! cell vertically below release point.
      if (this%idrape /= 0) &
        call this%dis%highest_active(ic, this%ibound)
      ! If returned cell is inactive, do not release particle
      if (this%ibound(ic) == 0) &
        particle%istatus = 8 ! permanently unreleased
    end if
    particle%x = x
    particle%y = y
    particle%z = z
    particle%trelease = trelease
    ! Set stopping time to earlier of times specified by 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%idomain(1) = 0
    particle%iboundary(1) = 0
    particle%idomain(2) = ic
    particle%iboundary(2) = 0
    particle%idomain(3) = 0
    particle%iboundary(3) = 0
    particle%ifrctrn = this%ifrctrn
    particle%iexmeth = this%iexmeth
    particle%iextend = this%iextend
    particle%extol = this%extol
 
    ! Store the particle's state and deallocate it
    call this%particles%save_particle(particle, np)
    deallocate (particle)
 
    ! Accumulate particle mass for release point
    this%rptm(ip) = this%rptm(ip) + done
 

Here is the call graph for this function:

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 367 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_unit(this%inunit, terminate=.true.)
    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_unit(this%inunit, terminate=.true.)
    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_unit(this%inunit, terminate=.true.)
    end if
    deallocate (polyverts)

Here is the call graph for this function:

Variable Documentation

ftype

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

private

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

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

text

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

private

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

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