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

Data Types

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

Functions/Subroutines

subroutine, public prp_create (packobj, id, ibcnum, inunit, iout, namemodel, pakname, fmi)
 Create a new particle release point package. More...
 
subroutine prp_da (this)
 Deallocate memory. More...
 
subroutine prp_set_pointers (this, ibound, izone, trackfilectl)
 @ 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.

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

◆ log_release()

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

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

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

◆ prp_ad()

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

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

310  class(PrtPrpType) :: this
311  integer(I4B) :: ip, it
312 
313  ! Notes
314  ! -----
315  ! Each release point can be thought of as
316  ! a gumball machine with infinite supply:
317  ! a point can release an arbitrary number
318  ! of particles, but only one at any time.
319  ! Coincident release times are merged to
320  ! a single time by the release scheduler.
321 
322  ! Reset mass accumulators for this time step.
323  do ip = 1, this%nreleasepoints
324  this%rptm(ip) = dzero
325  end do
326 
327  ! Advance the release schedule and check if
328  ! any releases will be made this time step.
329  call this%schedule%advance()
330  if (.not. this%schedule%any()) return
331 
332  ! Log the schedule to the list file.
333  call this%log_release()
334 
335  ! Expand the particle store. We know from the
336  ! schedule how many particles will be released.
337  call this%particles%resize( &
338  this%particles%num_stored() + &
339  (this%nreleasepoints * this%schedule%count()), &
340  this%memoryPath)
341 
342  ! Release a particle from each point for
343  ! each release time in the current step.
344  do ip = 1, this%nreleasepoints
345  do it = 1, this%schedule%count()
346  call this%release(ip, this%schedule%times(it))
347  end do
348  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.

204  ! dummy
205  class(PrtPrpType) :: this
206  integer(I4B), dimension(:), pointer, contiguous, optional :: nodelist
207  real(DP), dimension(:, :), pointer, contiguous, optional :: auxvar
208  ! local
209  integer(I4B) :: nps
210 
211  ! Allocate particle store, starting with the number
212  ! of release points (arrays resized if/when needed)
213  call allocate_particle_store( &
214  this%particles, &
215  this%nreleasepoints, &
216  this%memoryPath)
217 
218  ! Allocate arrays
219  call mem_allocate(this%rptx, this%nreleasepoints, 'RPTX', this%memoryPath)
220  call mem_allocate(this%rpty, this%nreleasepoints, 'RPTY', this%memoryPath)
221  call mem_allocate(this%rptz, this%nreleasepoints, 'RPTZ', this%memoryPath)
222  call mem_allocate(this%rptm, this%nreleasepoints, 'RPTMASS', this%memoryPath)
223  call mem_allocate(this%rptnode, this%nreleasepoints, 'RPTNODER', &
224  this%memoryPath)
225  call mem_allocate(this%rptname, lenboundname, this%nreleasepoints, &
226  'RPTNAME', this%memoryPath)
227 
228  ! Initialize arrays
229  do nps = 1, this%nreleasepoints
230  this%rptm(nps) = dzero
231  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.

236  class(PrtPrpType) :: this
237 
238  ! Allocate parent's scalars
239  call this%BndType%allocate_scalars()
240 
241  ! Allocate scalars for this type
242  call mem_allocate(this%ilocalz, 'ILOCALZ', this%memoryPath)
243  call mem_allocate(this%iextend, 'IEXTEND', this%memoryPath)
244  call mem_allocate(this%offset, 'OFFSET', this%memoryPath)
245  call mem_allocate(this%stoptime, 'STOPTIME', this%memoryPath)
246  call mem_allocate(this%stoptraveltime, 'STOPTRAVELTIME', this%memoryPath)
247  call mem_allocate(this%istopweaksink, 'ISTOPWEAKSINK', this%memoryPath)
248  call mem_allocate(this%istopzone, 'ISTOPZONE', this%memoryPath)
249  call mem_allocate(this%idrape, 'IDRAPE', this%memoryPath)
250  call mem_allocate(this%nreleasepoints, 'NRELEASEPOINTS', this%memoryPath)
251  call mem_allocate(this%nreleasetimes, 'NRELEASETIMES', this%memoryPath)
252  call mem_allocate(this%nparticles, 'NPARTICLES', this%memoryPath)
253  call mem_allocate(this%itrkout, 'ITRKOUT', this%memoryPath)
254  call mem_allocate(this%itrkhdr, 'ITRKHDR', this%memoryPath)
255  call mem_allocate(this%itrkcsv, 'ITRKCSV', this%memoryPath)
256  call mem_allocate(this%irlstls, 'IRLSTLS', this%memoryPath)
257  call mem_allocate(this%ifrctrn, 'IFRCTRN', this%memoryPath)
258  call mem_allocate(this%iexmeth, 'IEXMETH', this%memoryPath)
259  call mem_allocate(this%extol, 'EXTOL', this%memoryPath)
260  call mem_allocate(this%rttol, 'RTTOL', this%memoryPath)
261  call mem_allocate(this%rtfreq, 'RTFREQ', this%memoryPath)
262  call mem_allocate(this%foundtol, 'FOUNDTOL', this%memoryPath)
263 
264  ! Set values
265  this%ilocalz = 0
266  this%iextend = 0
267  this%offset = dzero
268  this%stoptime = huge(1d0)
269  this%stoptraveltime = huge(1d0)
270  this%istopweaksink = 0
271  this%istopzone = 0
272  this%idrape = 0
273  this%nreleasepoints = 0
274  this%nreleasetimes = 0
275  this%nparticles = 0
276  this%itrkout = 0
277  this%itrkhdr = 0
278  this%itrkcsv = 0
279  this%irlstls = 0
280  this%ifrctrn = 0
281  this%iexmeth = 0
282  this%extol = dzero
283  this%rttol = dsame * dep9
284  this%rtfreq = dzero
285  this%foundtol = .false.
286 

◆ prp_ar()

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

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

291  ! dummy variables
292  class(PrtPrpType), intent(inout) :: this
293  ! local variables
294  integer(I4B) :: n
295 
296  call this%obs%obs_ar()
297  call this%BndType%allocate_arrays()
298  if (this%inamedbound /= 0) then
299  do n = 1, this%nreleasepoints
300  this%boundname(n) = this%rptname(n)
301  end do
302  end if
303  do n = 1, this%nreleasepoints
304  this%nodelist(n) = this%rptnode(n)
305  end do

◆ prp_cq_simrate()

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

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

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

◆ prp_create()

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

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

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

147  class(PrtPrpType) :: this
148 
149  ! Deallocate parent
150  call this%BndType%bnd_da()
151 
152  ! Deallocate scalars
153  call mem_deallocate(this%ilocalz)
154  call mem_deallocate(this%iextend)
155  call mem_deallocate(this%offset)
156  call mem_deallocate(this%stoptime)
157  call mem_deallocate(this%stoptraveltime)
158  call mem_deallocate(this%istopweaksink)
159  call mem_deallocate(this%istopzone)
160  call mem_deallocate(this%idrape)
161  call mem_deallocate(this%nreleasepoints)
162  call mem_deallocate(this%nreleasetimes)
163  call mem_deallocate(this%nparticles)
164  call mem_deallocate(this%itrkout)
165  call mem_deallocate(this%itrkhdr)
166  call mem_deallocate(this%itrkcsv)
167  call mem_deallocate(this%irlstls)
168  call mem_deallocate(this%ifrctrn)
169  call mem_deallocate(this%iexmeth)
170  call mem_deallocate(this%extol)
171  call mem_deallocate(this%rttol)
172  call mem_deallocate(this%rtfreq)
173  call mem_deallocate(this%foundtol)
174 
175  ! Deallocate arrays
176  call mem_deallocate(this%rptx)
177  call mem_deallocate(this%rpty)
178  call mem_deallocate(this%rptz)
179  call mem_deallocate(this%rptnode)
180  call mem_deallocate(this%rptm)
181  call mem_deallocate(this%rptname, 'RPTNAME', this%memoryPath)
182 
183  ! Deallocate particle store, release time and release step selections
184  call this%particles%deallocate(this%memoryPath)
185  call this%schedule%deallocate()
186  deallocate (this%particles)
187  deallocate (this%schedule)

◆ prp_df_obs()

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

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

644  ! dummy
645  class(PrtPrpType) :: this
646  ! local
647  integer(I4B) :: indx
648  call this%obs%StoreObsType('prp', .true., indx)
649  this%obs%obsData(indx)%ProcessIdPtr => defaultobsidprocessor
650 
651  ! Store obs type and assign procedure pointer
652  ! for to-mvr observation type.
653  call this%obs%StoreObsType('to-mvr', .true., indx)
654  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.

1038  ! modules
1039  use tdismodule, only: totalsimtime
1040  ! dummy
1041  class(PrtPrpType), intent(inout) :: this
1042  ! local
1043  real(DP), allocatable :: times(:)
1044 
1045  ! check if a release time frequency is configured
1046  if (this%rtfreq <= dzero) return
1047 
1048  ! create array of regularly-spaced release times
1049  times = arange( &
1050  start=dzero, &
1051  stop=totalsimtime, &
1052  step=this%rtfreq)
1053 
1054  ! register times with release schedule
1055  call this%schedule%time_select%extend(times)
1056 
1057  ! make sure times strictly increase
1058  if (.not. this%schedule%time_select%increasing()) then
1059  errmsg = "Release times must strictly increase"
1060  call store_error(errmsg)
1061  call this%parser%StoreErrorUnit(terminate=.true.)
1062  end if
1063 
1064  ! deallocate
1065  deallocate (times)
1066 
real(dp), pointer, public totalsimtime
time at end of simulation
Definition: tdis.f90:37
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.

638  class(PrtPrpType) :: this
639  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.

659  use openspecmodule, only: access, form
660  use constantsmodule, only: maxcharlen, dzero
663  ! dummy
664  class(PrtPrpType), intent(inout) :: this
665  character(len=*), intent(inout) :: option
666  logical(LGP), intent(inout) :: found
667  ! locals
668  character(len=MAXCHARLEN) :: fname
669  character(len=MAXCHARLEN) :: keyword
670  ! formats
671  character(len=*), parameter :: fmttrkbin = &
672  "(4x, 'PARTICLE TRACKS WILL BE SAVED TO BINARY FILE: ', a, /4x, &
673  &'OPENED ON UNIT: ', I0)"
674  character(len=*), parameter :: fmttrkcsv = &
675  "(4x, 'PARTICLE TRACKS WILL BE SAVED TO CSV FILE: ', a, /4x, &
676  &'OPENED ON UNIT: ', I0)"
677 
678  select case (option)
679  case ('STOPTIME')
680  this%stoptime = this%parser%GetDouble()
681  found = .true.
682  case ('STOPTRAVELTIME')
683  this%stoptraveltime = this%parser%GetDouble()
684  found = .true.
685  case ('STOP_AT_WEAK_SINK')
686  this%istopweaksink = 1
687  found = .true.
688  case ('ISTOPZONE')
689  this%istopzone = this%parser%GetInteger()
690  found = .true.
691  case ('DRAPE')
692  this%idrape = 1
693  found = .true.
694  case ('TRACK')
695  call this%parser%GetStringCaps(keyword)
696  if (keyword == 'FILEOUT') then
697  ! parse filename
698  call this%parser%GetString(fname)
699  ! open binary output file
700  this%itrkout = getunit()
701  call openfile(this%itrkout, this%iout, fname, 'DATA(BINARY)', &
702  form, access, filstat_opt='REPLACE', &
703  mode_opt=mnormal)
704  write (this%iout, fmttrkbin) trim(adjustl(fname)), this%itrkout
705  ! open and write ascii header spec file
706  this%itrkhdr = getunit()
707  fname = trim(fname)//'.hdr'
708  call openfile(this%itrkhdr, this%iout, fname, 'CSV', &
709  filstat_opt='REPLACE', mode_opt=mnormal)
710  write (this%itrkhdr, '(a,/,a)') trackheader, trackdtypes
711  else
712  call store_error('OPTIONAL TRACK KEYWORD MUST BE '// &
713  'FOLLOWED BY FILEOUT')
714  end if
715  found = .true.
716  case ('TRACKCSV')
717  call this%parser%GetStringCaps(keyword)
718  if (keyword == 'FILEOUT') then
719  ! parse filename
720  call this%parser%GetString(fname)
721  ! open CSV output file and write headers
722  this%itrkcsv = getunit()
723  call openfile(this%itrkcsv, this%iout, fname, 'CSV', &
724  filstat_opt='REPLACE')
725  write (this%iout, fmttrkcsv) trim(adjustl(fname)), this%itrkcsv
726  write (this%itrkcsv, '(a)') trackheader
727  else
728  call store_error('OPTIONAL TRACKCSV KEYWORD MUST BE &
729  &FOLLOWED BY FILEOUT')
730  end if
731  found = .true.
732  case ('LOCAL_Z')
733  this%ilocalz = 1
734  found = .true.
735  case ('EXTEND_TRACKING')
736  this%iextend = 1
737  found = .true.
738  case ('EXIT_SOLVE_TOLERANCE')
739  this%extol = this%parser%GetDouble()
740  if (this%extol <= dzero) &
741  call store_error('EXIT_SOLVE_TOLERANCE MUST BE POSITIVE')
742  found = .true.
743  this%foundtol = .true.
744  case ('RELEASE_TIME_TOLERANCE')
745  this%rttol = this%parser%GetDouble()
746  if (this%rttol <= dzero) &
747  call store_error('RELEASE_TIME_TOLERANCE MUST BE POSITIVE')
748  found = .true.
749  case ('RELEASE_TIME_FREQUENCY')
750  this%rtfreq = this%parser%GetDouble()
751  if (this%rtfreq <= dzero) &
752  call store_error('RELEASE_TIME_FREQUENCY MUST BE POSITIVE')
753  found = .true.
754  case ('DEV_FORCETERNARY')
755  call this%parser%DevOpt()
756  this%ifrctrn = 1
757  write (this%iout, '(4x,a)') &
758  'TRACKING WILL BE DONE USING THE TERNARY METHOD REGARDLESS OF CELL TYPE'
759  found = .true.
760  case ('DEV_EXIT_SOLVE_METHOD')
761  call this%parser%DevOpt()
762  this%iexmeth = this%parser%GetInteger()
763  if (.not. (this%iexmeth /= 1 .or. this%iexmeth /= 2)) &
764  call store_error('DEV_EXIT_SOLVE_METHOD MUST BE &
765  &1 (BRENT) OR 2 (CHANDRUPATLA)')
766  found = .true.
767  case default
768  found = .false.
769  end select
770 
771  ! Catch unrecognized options
772  if (.not. found) then
773  errmsg = "UNKNOWN PRP OPTION '"//trim(keyword)//"'."
774  call store_error(errmsg)
775  call this%parser%StoreErrorUnit()
776  end if
777 
778  ! Create release schedule now that we know
779  ! the coincident release time tolerance
780  this%schedule => create_release_schedule(tol=this%rttol)
781 
This module contains simulation constants.
Definition: Constants.f90:9
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
integer(i4b), parameter maxcharlen
maximum length of char string
Definition: Constants.f90:47
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
subroutine, public urword(line, icol, istart, istop, ncode, n, r, iout, in)
Extract a word from a string.
character(len=20) access
Definition: OpenSpec.f90:7
character(len=20) form
Definition: OpenSpec.f90:7
character(len= *), parameter, public trackheader
Definition: TrackData.f90:56
character(len= *), parameter, public trackdtypes
Definition: TrackData.f90:61
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.

786  ! dummy
787  class(PrtPrpType), intent(inout) :: this
788  ! local
789  character(len=LINELENGTH) :: errmsg, keyword
790  integer(I4B) :: ierr
791  logical :: isfound, endOfBlock
792 
793  if (.not. this%foundtol) &
794  call store_error('EXIT_SOLVE_TOLERANCE MISSING, VALUE REQUIRED')
795 
796  ! get dimension block
797  call this%parser%GetBlock('DIMENSIONS', isfound, ierr, &
798  supportopenclose=.true.)
799 
800  ! parse dimension block if detected
801  if (isfound) then
802  write (this%iout, '(1x,a)') 'PROCESSING PARTICLE INPUT DIMENSIONS'
803  do
804  call this%parser%GetNextLine(endofblock)
805  if (endofblock) exit
806  call this%parser%GetStringCaps(keyword)
807  select case (keyword)
808  case ('NRELEASEPTS')
809  this%nreleasepoints = this%parser%GetInteger()
810  case ('NRELEASETIMES')
811  this%nreleasetimes = this%parser%GetInteger()
812  case default
813  write (errmsg, &
814  '(4x,a,a)') '****ERROR. UNKNOWN PARTICLE INPUT DIMENSION: ', &
815  trim(keyword)
816  call store_error(errmsg)
817  call this%parser%StoreErrorUnit()
818  end select
819  end do
820  write (this%iout, '(1x,a)') 'END OF PARTICLE INPUT DIMENSIONS'
821  else
822  call store_error('ERROR. REQUIRED DIMENSIONS BLOCK NOT FOUND.')
823  end if
824 
825  ! set maxbound and nbound to nreleasepts
826  this%maxbound = this%nreleasepoints
827  this%nbound = this%nreleasepoints
828 
829  ! allocate arrays for prp package
830  call this%prp_allocate_arrays()
831 
832  ! read packagedata and releasetimes blocks
833  call this%prp_read_packagedata()
834  call this%prp_read_releasetimes()
835  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.

840  ! dummy
841  class(PrtPrpType), intent(inout) :: this
842  ! local
843  character(len=LINELENGTH) :: cellid
844  character(len=LENBOUNDNAME) :: bndName, bndNameTemp
845  character(len=9) :: cno
846  logical :: isfound
847  logical :: endOfBlock
848  integer(I4B) :: ival
849  integer(I4B) :: n
850  integer(I4B) :: ierr
851  character(len=LENBOUNDNAME), dimension(:), allocatable :: nametxt
852  integer(I4B), dimension(:), allocatable :: nboundchk
853  integer(I4B), dimension(:), allocatable :: noder
854  real(DP), dimension(:), allocatable :: x
855  real(DP), dimension(:), allocatable :: y
856  real(DP), dimension(:), allocatable :: z
857  real(DP), dimension(:), allocatable :: tstop
858  ! format
859  character(len=*), parameter :: fmttend = &
860  "('end time (', G0, ') must be greater than or equal to the &
861  &begin time (', G0, ').')"
862 
863  ! allocate temporary variables
864  allocate (noder(this%nreleasepoints))
865  allocate (x(this%nreleasepoints))
866  allocate (y(this%nreleasepoints))
867  allocate (z(this%nreleasepoints))
868  allocate (tstop(this%nreleasepoints))
869  allocate (nametxt(this%nreleasepoints))
870  allocate (nboundchk(this%nreleasepoints))
871 
872  ! initialize temporary variables
873  do n = 1, this%nreleasepoints
874  nboundchk(n) = 0
875  end do
876 
877  ! read particle release point data
878  ! get particle release points block
879  call this%parser%GetBlock('PACKAGEDATA', isfound, ierr, &
880  supportopenclose=.true.)
881 
882  ! parse block if detected
883  if (isfound) then
884  write (this%iout, '(/1x,a)') 'PROCESSING '//trim(adjustl(this%packName)) &
885  //' PACKAGEDATA'
886  do
887  call this%parser%GetNextLine(endofblock)
888  if (endofblock) exit
889  ival = this%parser%GetInteger()
890  n = ival
891 
892  if (n < 1 .or. n > this%nreleasepoints) then
893  write (errmsg, '(a,1x,i0,a)') &
894  'Release point number must be greater than 0 and less than ', &
895  'or equal to', this%nreleasepoints, '.'
896  call store_error(errmsg)
897  cycle
898  end if
899 
900  ! increment nboundchk
901  nboundchk(n) = nboundchk(n) + 1
902 
903  ! node number
904  call this%parser%GetCellid(this%dis%ndim, cellid)
905  noder(n) = this%dis%noder_from_cellid(cellid, this%inunit, this%iout)
906 
907  ! x, y, z coordinates
908  x(n) = this%parser%GetDouble()
909  y(n) = this%parser%GetDouble()
910  z(n) = this%parser%GetDouble()
911 
912  if (this%ilocalz > 0 .and. (z(n) < 0 .or. z(n) > 1)) then
913  call store_error('Local z coordinate must fall in the interval [0, 1]')
914  cycle
915  end if
916 
917  ! set default boundname
918  write (cno, '(i9.9)') n
919  bndname = 'PRP'//cno
920 
921  ! read boundnames from file, if provided
922  if (this%inamedbound /= 0) then
923  call this%parser%GetStringCaps(bndnametemp)
924  if (bndnametemp /= '') &
925  bndname = bndnametemp
926  else
927  bndname = ''
928  end if
929 
930  ! store temp boundnames
931  nametxt(n) = bndname
932  end do
933 
934  write (this%iout, '(1x,a)') &
935  'END OF '//trim(adjustl(this%packName))//' PACKAGEDATA'
936 
937  ! check for duplicate or missing particle release points
938  do n = 1, this%nreleasepoints
939  if (nboundchk(n) == 0) then
940  write (errmsg, '(a,a,1x,i0,a)') 'No data specified for particle ', &
941  'release point', n, '.'
942  call store_error(errmsg)
943  else if (nboundchk(n) > 1) then
944  write (errmsg, '(a,1x,i0,1x,a,1x,i0,1x,a)') &
945  'Data for particle release point', n, 'specified', nboundchk(n), &
946  'times.'
947  call store_error(errmsg)
948  end if
949  end do
950  else
951  call store_error('Required packagedata block not found.')
952  end if
953 
954  ! terminate if any errors were detected
955  if (count_errors() > 0) then
956  call this%parser%StoreErrorUnit()
957  end if
958 
959  ! fill particle release point data with data stored in temporary local arrays
960  do n = 1, this%nreleasepoints
961  this%rptnode(n) = noder(n)
962  this%rptx(n) = x(n)
963  this%rpty(n) = y(n)
964  this%rptz(n) = z(n)
965  this%rptname(n) = nametxt(n)
966  end do
967 
968  ! deallocate local storage
969  deallocate (noder)
970  deallocate (x)
971  deallocate (y)
972  deallocate (z)
973  deallocate (tstop)
974  deallocate (nametxt)
975  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.

980  ! dummy
981  class(PrtPrpType), intent(inout) :: this
982  ! local
983  integer(I4B) :: i, ierr
984  logical(LGP) :: eob, found, success
985  real(DP) :: t
986  real(DP), allocatable :: times(:)
987 
988  ! get releasetimes block
989  call this%parser%GetBlock('RELEASETIMES', found, ierr, &
990  supportopenclose=.true., &
991  blockrequired=.false.)
992 
993  ! raise an error if releasetimes has a dimension
994  ! but no block was found, otherwise return early
995  if (.not. found) then
996  if (this%nreleasetimes <= 0) return
997  write (errmsg, '(a, i0)') &
998  "Expected RELEASETIMES with length ", this%nreleasetimes
999  call store_error(errmsg)
1000  call this%parser%StoreErrorUnit(terminate=.true.)
1001  end if
1002 
1003  ! allocate times array
1004  allocate (times(this%nreleasetimes))
1005 
1006  ! read times from the block
1007  write (this%iout, '(/1x,a)') &
1008  'PROCESSING '//trim(adjustl(this%text))//' RELEASETIMES'
1009  do i = 1, this%nreleasetimes
1010  call this%parser%GetNextLine(eob)
1011  if (eob) exit
1012  call this%parser%TryGetDouble(t, success)
1013  if (.not. success) then
1014  errmsg = "Failed to read double precision value"
1015  call store_error(errmsg)
1016  call this%parser%StoreErrorUnit(terminate=.true.)
1017  end if
1018  times(i) = t
1019  end do
1020 
1021  ! register times with the release schedule
1022  call this%schedule%time_select%extend(times)
1023 
1024  ! make sure times strictly increase
1025  if (.not. this%schedule%time_select%increasing()) then
1026  errmsg = "Release times must strictly increase"
1027  call store_error(errmsg)
1028  call this%parser%StoreErrorUnit(terminate=.true.)
1029  end if
1030 
1031  ! deallocate
1032  deallocate (times)
1033 
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.

510  ! modules
511  use tdismodule, only: kper, nper
512  use inputoutputmodule, only: urword
513  ! dummy variables
514  class(PrtPrpType), intent(inout) :: this
515  ! local variables
516  integer(I4B) :: ierr
517  logical(LGP) :: is_found
518  logical(LGP) :: end_of_block
519  logical(LGP) :: no_blocks
520  character(len=LINELENGTH) :: line
521  character(len=LINELENGTH), allocatable :: lines(:)
522  ! formats
523  character(len=*), parameter :: fmtblkerr = &
524  "('Looking for BEGIN PERIOD iper. &
525  &Found ', a, ' instead.')"
526  character(len=*), parameter :: fmt_steps = &
527  "(6x,'TIME STEP(S) ',50(I0,' '))" ! 50 limit is similar to STEPS in OC
528  character(len=*), parameter :: fmt_freq = &
529  "(6x,'EVERY ',I0,' TIME STEP(S)')"
530  character(len=*), parameter :: fmt_fracs = &
531  "(6x,50(f10.3,' '))"
532 
533  ! Set ionper to the stress period number for which a new block of data
534  ! will be read.
535  if (this%inunit == 0) return
536 
537  ! get stress period data
538  no_blocks = .false.
539  if (this%ionper < kper) then
540  ! get period block
541  call this%parser%GetBlock('PERIOD', is_found, ierr, &
542  supportopenclose=.true., &
543  blockrequired=.false.)
544  if (is_found) then
545  ! read ionper and check for increasing period numbers
546  call this%read_check_ionper()
547  else
548  ! PERIOD block not found
549  if (ierr < 0) then
550  if (kper == 1) then
551  ! End of file found; no period data for the simulation.
552  no_blocks = .true.
553  else
554  ! End of file found; no more period data.
555  this%ionper = nper + 1
556  end if
557  else
558  ! Found invalid block
559  call this%parser%GetCurrentLine(line)
560  write (errmsg, fmtblkerr) adjustl(trim(line))
561  call store_error(errmsg, terminate=.true.)
562  end if
563  end if
564  end if
565 
566  ! If the user hasn't provided any release settings (neither
567  ! explicit release times, release time frequency, or period
568  ! block release settings), default to a single release at the
569  ! start of the first period's first time step.
570  if (no_blocks .and. &
571  kper == 1 .and. &
572  size(this%schedule%time_select%times) == 0) then
573  allocate (lines(1))
574  line = "FIRST"
575  lines(1) = line
576  call this%schedule%advance(lines=lines)
577  else if (this%ionper == kper) then
578  ! If the current stress period matches the
579  ! block we are reading, parse the setting
580  ! and register it with the schedule.
581  allocate (lines(0))
582  recordloop: do
583  call this%parser%GetNextLine(end_of_block)
584  if (end_of_block) exit recordloop
585  call this%parser%GetCurrentLine(line)
586  call expandarray(lines)
587  lines(size(lines)) = line
588  end do recordloop
589  if (size(lines) > 0) &
590  call this%schedule%advance(lines=lines)
591  deallocate (lines)
592  end if
593 
integer(i4b), pointer, public kper
current stress period number
Definition: tdis.f90:23
integer(i4b), pointer, public nper
number of stress period
Definition: tdis.f90:21
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(trackfilecontroltype), pointer  trackfilectl 
)
private

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

192  class(PrtPrpType) :: this
193  integer(I4B), dimension(:), pointer, contiguous :: ibound
194  integer(I4B), dimension(:), pointer, contiguous :: izone
195  type(TrackFileControlType), pointer :: trackfilectl
196 
197  this%ibound => ibound
198  this%rptzone => izone
199  this%tracks => trackfilectl

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

411  ! dummy
412  class(PrtPrpType), intent(inout) :: this !< this instance
413  integer(I4B), intent(in) :: ip !< particle index
414  real(DP), intent(in) :: trelease !< release time
415  ! local
416  integer(I4B) :: irow, icol, ilay, icpl
417  integer(I4B) :: ic, icu
418  integer(I4B) :: np
419  real(DP) :: x, y, z
420  real(DP) :: top, bot, hds
421  type(ParticleType), pointer :: particle
422 
423  ! Increment cumulative particle count
424  np = this%nparticles + 1
425  this%nparticles = np
426 
427  ! Get reduced and user node numbers
428  ic = this%rptnode(ip)
429  icu = this%dis%get_nodeuser(ic)
430 
431  ! Load coordinates and transform if needed
432  x = this%rptx(ip)
433  y = this%rpty(ip)
434  if (this%ilocalz > 0) then
435  top = this%fmi%dis%top(ic)
436  bot = this%fmi%dis%bot(ic)
437  hds = this%fmi%gwfhead(ic)
438  z = bot + this%rptz(ip) * (hds - bot)
439  else
440  z = this%rptz(ip)
441  end if
442 
443  ! Make sure release point is in the grid/cell
444  call this%validate_release_point(ic, x, y, z)
445 
446  ! Initialize the particle
447  call create_particle(particle)
448  if (size(this%boundname) /= 0) then
449  particle%name = this%boundname(ip)
450  else
451  particle%name = ''
452  end if
453  particle%irpt = ip
454  particle%istopweaksink = this%istopweaksink
455  particle%istopzone = this%istopzone
456  particle%icu = icu
457  select type (dis => this%dis)
458  type is (distype)
459  call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, irow, icol, ilay)
460  type is (disvtype)
461  call get_jk(icu, dis%ncpl, dis%nlay, icpl, ilay)
462  end select
463  particle%ilay = ilay
464  particle%izone = this%rptzone(ic)
465  particle%istatus = 0
466  ! Handle inactive cells
467  if (this%ibound(ic) == 0) then
468  ! If drape option activated, release in highest active
469  ! cell vertically below release point.
470  if (this%idrape /= 0) &
471  call this%dis%highest_active(ic, this%ibound)
472  ! If returned cell is inactive, do not release particle
473  if (this%ibound(ic) == 0) &
474  particle%istatus = 8 ! permanently unreleased
475  end if
476  particle%x = x
477  particle%y = y
478  particle%z = z
479  particle%trelease = trelease
480  ! Set stopping time to earlier of times specified by STOPTIME and STOPTRAVELTIME
481  if (this%stoptraveltime == huge(1d0)) then
482  particle%tstop = this%stoptime
483  else
484  particle%tstop = particle%trelease + this%stoptraveltime
485  if (this%stoptime < particle%tstop) particle%tstop = this%stoptime
486  end if
487  particle%ttrack = particle%trelease
488  particle%idomain(1) = 0
489  particle%iboundary(1) = 0
490  particle%idomain(2) = ic
491  particle%iboundary(2) = 0
492  particle%idomain(3) = 0
493  particle%iboundary(3) = 0
494  particle%ifrctrn = this%ifrctrn
495  particle%iexmeth = this%iexmeth
496  particle%iextend = this%iextend
497  particle%extol = this%extol
498 
499  ! Store the particle's state and deallocate it
500  call this%particles%save_particle(particle, np)
501  deallocate (particle)
502 
503  ! Accumulate particle mass for release point
504  this%rptm(ip) = this%rptm(ip) + done
505 
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]thisthis instance
[in]iccell index
[in]zrelease point

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

368  class(PrtPrpType), intent(inout) :: this !< this instance
369  integer(I4B), intent(in) :: ic !< cell index
370  real(DP), intent(in) :: x, y, z !< release point
371  ! local
372  real(DP), allocatable :: polyverts(:, :)
373 
374  call this%fmi%dis%get_polyverts(ic, polyverts)
375  if (.not. point_in_polygon(x, y, polyverts)) then
376  write (errmsg, '(a,g0,a,g0,a,i0)') &
377  'Error: release point (x=', x, ', y=', y, ') is not in cell ', &
378  this%dis%get_nodeuser(ic)
379  call store_error(errmsg, terminate=.false.)
380  call store_error_unit(this%inunit, terminate=.true.)
381  end if
382  if (z > maxval(this%dis%top)) then
383  write (errmsg, '(a,g0,a,g0,a,i0)') &
384  'Error: release point (z=', z, ') is above grid top ', &
385  maxval(this%dis%top)
386  call store_error(errmsg, terminate=.false.)
387  call store_error_unit(this%inunit, terminate=.true.)
388  else if (z < minval(this%dis%bot)) then
389  write (errmsg, '(a,g0,a,g0,a,i0)') &
390  'Error: release point (z=', z, ') is below grid bottom ', &
391  minval(this%dis%bot)
392  call store_error(errmsg, terminate=.false.)
393  call store_error_unit(this%inunit, terminate=.true.)
394  end if
395  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.

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

◆ text

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

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

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