MODFLOW 6  version 6.8.0.dev0
USGS Modular Hydrologic Model
Data Types | Enumerations | Functions/Subroutines | Variables
particlemodule Module Reference

Data Types

type  particletype
 Particle tracked by the PRT model. More...
 
type  particlestoretype
 Structure of arrays to store particles. More...
 

Enumerations

enum  {
  active = 1 , term_boundary = 2 , term_weaksink = 3 , term_no_exits = 5 ,
  term_stopzone = 6 , term_inactive = 7 , term_unreleased = 8 , term_no_exits_sub = 9 ,
  term_timeout = 10
}
 Particle status enumeration. More...
 

Functions/Subroutines

subroutine create_particle (particle)
 Create a new particle. More...
 
subroutine create_particle_store (store, np, mempath)
 Allocate particle store. More...
 
subroutine destroy (this, mempath)
 Destroy particle store after use. More...
 
subroutine destroy_particle (particle)
 Destroy a particle after use. More...
 
subroutine resize (this, np, mempath)
 Reallocate particle storage to the given size. More...
 
subroutine get (this, particle, imdl, iprp, ip)
 Load a particle from the particle store. More...
 
subroutine put (this, particle, ip)
 Save a particle's state to the particle store. More...
 
subroutine copy_from (this, source)
 Copy all particle state from source to this store. The two particle stores must already be of the same size. More...
 
subroutine transform_coords (this, xorigin, yorigin, zorigin, sinrot, cosrot, invert)
 Transform particle coordinates. More...
 
subroutine reset_transform (this)
 Reset particle coordinate transformation properties. More...
 
subroutine get_model_coords (this, x, y, z)
 Return the particle's model coordinates, inverting any applied transformation if needed. The particle's state is not altered. More...
 
integer function num_stored (this)
 Return the number of particles. More...
 
character(len=:) function, allocatable get_id (this)
 Get a string identifier for the particle. More...
 

Variables

integer, parameter max_level = 4
 Tracking "levels" defined in method modules. Currently only 3 used. More...
 

Enumeration Type Documentation

◆ anonymous enum

anonymous enum

Particles begin in status 1 (active) at release time. Status may only increase over time. Status values greater than one imply termination. A particle may terminate for several reasons, all mutually exclusive. A particle's final tracking status will always be greater than one.

Status codes 0-3 and 5-8 correspond directly to MODPATH 7 status codes. Code 4 does not apply to PRT because PRT does not distinguish forwards from backwards tracking. Status code 9 provides more specific, subcell- level information about a particle which terminates due to no outflow. Code 10 distinguishes particles which have "timed out" upon reaching a user-specified stop time or the end of the simulation.

Enumerator
active 
term_boundary 

terminated at a boundary face

term_weaksink 

terminated in a weak sink cell

term_no_exits 

terminated in a cell with no exit face

term_stopzone 

terminated in a cell with a stop zone number

term_inactive 

terminated in an inactive cell

term_unreleased 

terminated permanently unreleased

term_no_exits_sub 

terminated in a subcell with no exit face

term_timeout 

terminated at stop time or end of simulation

Definition at line 31 of file Particle.f90.

Function/Subroutine Documentation

◆ copy_from()

subroutine particlemodule::copy_from ( class(particlestoretype), intent(inout)  this,
class(particlestoretype), intent(in)  source 
)

Definition at line 346 of file Particle.f90.

347  ! dummy
348  class(ParticleStoreType), intent(inout) :: this
349  class(ParticleStoreType), intent(in) :: source
350  ! local
351  integer(I4B) :: np
352 
353  np = source%num_stored()
354 
355  if (this%num_stored() /= np) then
356  write (errmsg, '(a,i0,a,i0,a)') &
357  'Cannot copy particle store, size mismatch (source=', np, &
358  ', dest=', this%num_stored(), ')'
359  call pstop(1, errmsg)
360  end if
361 
362  this%name(:) = source%name(:)
363  this%imdl(:) = source%imdl(:)
364  this%iprp(:) = source%iprp(:)
365  this%irpt(:) = source%irpt(:)
366  this%extend(:) = source%extend(:)
367  this%frctrn(:) = source%frctrn(:)
368  this%istopweaksink(:) = source%istopweaksink(:)
369  this%istopzone(:) = source%istopzone(:)
370  this%idrymeth(:) = source%idrymeth(:)
371  this%iexmeth(:) = source%iexmeth(:)
372  this%icycwin(:) = source%icycwin(:)
373  this%extol(:) = source%extol(:)
374  this%itrdomain(:, :) = source%itrdomain(:, :)
375  this%iboundary(:, :) = source%iboundary(:, :)
376  this%icu(:) = source%icu(:)
377  this%ilay(:) = source%ilay(:)
378  this%izone(:) = source%izone(:)
379  this%istatus(:) = source%istatus(:)
380  this%x(:) = source%x(:)
381  this%y(:) = source%y(:)
382  this%z(:) = source%z(:)
383  this%trelease(:) = source%trelease(:)
384  this%tstop(:) = source%tstop(:)
385  this%ttrack(:) = source%ttrack(:)
Here is the call graph for this function:

◆ create_particle()

subroutine particlemodule::create_particle ( type(particletype), pointer  particle)

Definition at line 152 of file Particle.f90.

153  type(ParticleType), pointer :: particle !< particle
154  allocate (particle)
155  allocate (particle%history)
Here is the caller graph for this function:

◆ create_particle_store()

subroutine particlemodule::create_particle_store ( type(particlestoretype), pointer  store,
integer(i4b), intent(in)  np,
character(*), intent(in)  mempath 
)
Parameters
[in]npnumber of particles
[in]mempathpath to memory

Definition at line 159 of file Particle.f90.

160  type(ParticleStoreType), pointer :: store !< store
161  integer(I4B), intent(in) :: np !< number of particles
162  character(*), intent(in) :: mempath !< path to memory
163 
164  allocate (store)
165  call mem_allocate(store%imdl, np, 'PLIMDL', mempath)
166  call mem_allocate(store%irpt, np, 'PLIRPT', mempath)
167  call mem_allocate(store%iprp, np, 'PLIPRP', mempath)
168  call mem_allocate(store%name, lenboundname, np, 'PLNAME', mempath)
169  call mem_allocate(store%icu, np, 'PLICU', mempath)
170  call mem_allocate(store%ilay, np, 'PLILAY', mempath)
171  call mem_allocate(store%izone, np, 'PLIZONE', mempath)
172  call mem_allocate(store%istatus, np, 'PLISTATUS', mempath)
173  call mem_allocate(store%x, np, 'PLX', mempath)
174  call mem_allocate(store%y, np, 'PLY', mempath)
175  call mem_allocate(store%z, np, 'PLZ', mempath)
176  call mem_allocate(store%trelease, np, 'PLTRELEASE', mempath)
177  call mem_allocate(store%tstop, np, 'PLTSTOP', mempath)
178  call mem_allocate(store%ttrack, np, 'PLTTRACK', mempath)
179  call mem_allocate(store%istopweaksink, np, 'PLISTOPWEAKSINK', mempath)
180  call mem_allocate(store%istopzone, np, 'PLISTOPZONE', mempath)
181  call mem_allocate(store%idrymeth, np, 'PLIDRYMETH', mempath)
182  call mem_allocate(store%frctrn, np, 'PLFRCTRN', mempath)
183  call mem_allocate(store%iexmeth, np, 'PLIEXMETH', mempath)
184  call mem_allocate(store%extol, np, 'PLEXTOL', mempath)
185  call mem_allocate(store%extend, np, 'PLEXTEND', mempath)
186  call mem_allocate(store%icycwin, np, 'PLICYCWIN', mempath)
187  call mem_allocate(store%itrdomain, np, max_level, 'PLIDOMAIN', mempath)
188  call mem_allocate(store%iboundary, np, max_level, 'PLIBOUNDARY', mempath)
Here is the caller graph for this function:

◆ destroy()

subroutine particlemodule::destroy ( class(particlestoretype), intent(inout)  this,
character(*), intent(in)  mempath 
)
Parameters
[in,out]thisstore
[in]mempathpath to memory

Definition at line 192 of file Particle.f90.

193  class(ParticleStoreType), intent(inout) :: this !< store
194  character(*), intent(in) :: mempath !< path to memory
195 
196  call mem_deallocate(this%imdl, 'PLIMDL', mempath)
197  call mem_deallocate(this%iprp, 'PLIPRP', mempath)
198  call mem_deallocate(this%irpt, 'PLIRPT', mempath)
199  call mem_deallocate(this%name, 'PLNAME', mempath)
200  call mem_deallocate(this%icu, 'PLICU', mempath)
201  call mem_deallocate(this%ilay, 'PLILAY', mempath)
202  call mem_deallocate(this%izone, 'PLIZONE', mempath)
203  call mem_deallocate(this%istatus, 'PLISTATUS', mempath)
204  call mem_deallocate(this%x, 'PLX', mempath)
205  call mem_deallocate(this%y, 'PLY', mempath)
206  call mem_deallocate(this%z, 'PLZ', mempath)
207  call mem_deallocate(this%trelease, 'PLTRELEASE', mempath)
208  call mem_deallocate(this%tstop, 'PLTSTOP', mempath)
209  call mem_deallocate(this%ttrack, 'PLTTRACK', mempath)
210  call mem_deallocate(this%istopweaksink, 'PLISTOPWEAKSINK', mempath)
211  call mem_deallocate(this%istopzone, 'PLISTOPZONE', mempath)
212  call mem_deallocate(this%idrymeth, 'PLIDRYMETH', mempath)
213  call mem_deallocate(this%frctrn, 'PLFRCTRN', mempath)
214  call mem_deallocate(this%iexmeth, 'PLIEXMETH', mempath)
215  call mem_deallocate(this%extol, 'PLEXTOL', mempath)
216  call mem_deallocate(this%extend, 'PLEXTEND', mempath)
217  call mem_deallocate(this%icycwin, 'PLICYCWIN', mempath)
218  call mem_deallocate(this%itrdomain, 'PLIDOMAIN', mempath)
219  call mem_deallocate(this%iboundary, 'PLIBOUNDARY', mempath)

◆ destroy_particle()

subroutine particlemodule::destroy_particle ( class(particletype), intent(inout)  particle)

Definition at line 223 of file Particle.f90.

224  class(ParticleType), intent(inout) :: particle !< particle
225  deallocate (particle%history)

◆ get()

subroutine particlemodule::get ( class(particlestoretype), intent(inout)  this,
class(particletype), intent(inout)  particle,
integer(i4b), intent(in)  imdl,
integer(i4b), intent(in)  iprp,
integer(i4b), intent(in)  ip 
)

This routine is used to initialize a particle for tracking. The advancing flag and coordinate transformation are reset.

Parameters
[in,out]thisparticle store
[in]imdlindex of model particle originated in
[in]iprpindex of particle release package particle originated in
[in]ipindex into the particle list

Definition at line 266 of file Particle.f90.

267  class(ParticleStoreType), intent(inout) :: this !< particle store
268  class(ParticleType), intent(inout) :: particle !< particle
269  integer(I4B), intent(in) :: imdl !< index of model particle originated in
270  integer(I4B), intent(in) :: iprp !< index of particle release package particle originated in
271  integer(I4B), intent(in) :: ip !< index into the particle list
272 
273  call particle%reset_transform()
274  call particle%history%Clear()
275  particle%imdl = imdl
276  particle%iprp = iprp
277  particle%irpt = this%irpt(ip)
278  particle%ip = ip
279  particle%name = this%name(ip)
280  particle%istopweaksink = this%istopweaksink(ip)
281  particle%istopzone = this%istopzone(ip)
282  particle%idrymeth = this%idrymeth(ip)
283  particle%icu = this%icu(ip)
284  particle%ilay = this%ilay(ip)
285  particle%izone = this%izone(ip)
286  particle%istatus = this%istatus(ip)
287  particle%x = this%x(ip)
288  particle%y = this%y(ip)
289  particle%z = this%z(ip)
290  particle%trelease = this%trelease(ip)
291  particle%tstop = this%tstop(ip)
292  particle%ttrack = this%ttrack(ip)
293  particle%advancing = .true.
294  particle%itrdomain(1:max_level) = &
295  this%itrdomain(ip, 1:max_level)
296  particle%itrdomain(1) = imdl
297  particle%iboundary(1:max_level) = &
298  this%iboundary(ip, 1:max_level)
299  particle%frctrn = this%frctrn(ip)
300  particle%iexmeth = this%iexmeth(ip)
301  particle%extol = this%extol(ip)
302  particle%extend = this%extend(ip)
303  particle%icycwin = this%icycwin(ip)

◆ get_id()

character(len=:) function, allocatable particlemodule::get_id ( class(particletype), intent(in)  this)

Definition at line 464 of file Particle.f90.

465  class(ParticleType), intent(in) :: this
466  character(len=:), allocatable :: str
467  ! local
468  character(len=LINELENGTH) :: temp
469 
470  write (temp, '(I0,1a,I0,1a,I0,1a,G0)') &
471  this%imdl, this%iprp, this%irpt, this%trelease
472  str = trim(adjustl(temp))

◆ get_model_coords()

subroutine particlemodule::get_model_coords ( class(particletype), intent(in)  this,
real(dp), intent(out)  x,
real(dp), intent(out)  y,
real(dp), intent(out)  z 
)
Parameters
[in]thisparticle
[out]xx coordinate
[out]yy coordinate
[out]zz coordinate

Definition at line 439 of file Particle.f90.

440  use geomutilmodule, only: transform, compose
441  class(ParticleType), intent(in) :: this !< particle
442  real(DP), intent(out) :: x !< x coordinate
443  real(DP), intent(out) :: y !< y coordinate
444  real(DP), intent(out) :: z !< z coordinate
445 
446  if (this%transformed) then
447  call transform(this%x, this%y, this%z, x, y, z, &
448  this%xorigin, this%yorigin, this%zorigin, &
449  this%sinrot, this%cosrot, invert=.true.)
450  else
451  x = this%x
452  y = this%y
453  z = this%z
454  end if
geomutilmodule::transform
subroutine, public transform(xin, yin, zin, xout, yout, zout, xorigin, yorigin, zorigin, sinrot, cosrot, invert)
Apply a 3D translation and optional 2D rotation to coordinates.
Definition: GeomUtil.f90:183
geomutilmodule::compose
subroutine, public compose(xorigin, yorigin, zorigin, sinrot, cosrot, xorigin_new, yorigin_new, zorigin_new, sinrot_new, cosrot_new, invert)
Apply a 3D translation and 2D rotation to an existing transformation.
Definition: GeomUtil.f90:243
Here is the call graph for this function:

◆ num_stored()

integer function particlemodule::num_stored ( class(particlestoretype)  this)

Definition at line 458 of file Particle.f90.

459  class(ParticleStoreType) :: this
460  n = size(this%imdl)

◆ put()

subroutine particlemodule::put ( class(particlestoretype), intent(inout)  this,
class(particletype), intent(in)  particle,
integer(i4b), intent(in)  ip 
)
Parameters
[in,out]thisparticle storage
[in]ipparticle index

Definition at line 307 of file Particle.f90.

308  class(ParticleStoreType), intent(inout) :: this !< particle storage
309  class(ParticleType), intent(in) :: particle !< particle
310  integer(I4B), intent(in) :: ip !< particle index
311 
312  this%imdl(ip) = particle%imdl
313  this%iprp(ip) = particle%iprp
314  this%irpt(ip) = particle%irpt
315  this%name(ip) = particle%name
316  this%istopweaksink(ip) = particle%istopweaksink
317  this%istopzone(ip) = particle%istopzone
318  this%idrymeth(ip) = particle%idrymeth
319  this%icu(ip) = particle%icu
320  this%ilay(ip) = particle%ilay
321  this%izone(ip) = particle%izone
322  this%istatus(ip) = particle%istatus
323  this%x(ip) = particle%x
324  this%y(ip) = particle%y
325  this%z(ip) = particle%z
326  this%trelease(ip) = particle%trelease
327  this%tstop(ip) = particle%tstop
328  this%ttrack(ip) = particle%ttrack
329  this%itrdomain( &
330  ip, &
331  1:max_level) = &
332  particle%itrdomain(1:max_level)
333  this%iboundary( &
334  ip, &
335  1:max_level) = &
336  particle%iboundary(1:max_level)
337  this%frctrn(ip) = particle%frctrn
338  this%iexmeth(ip) = particle%iexmeth
339  this%extol(ip) = particle%extol
340  this%extend(ip) = particle%extend
341  this%icycwin(ip) = particle%icycwin

◆ reset_transform()

subroutine particlemodule::reset_transform ( class(particletype), intent(inout)  this)
Parameters
[in,out]thisparticle

Definition at line 424 of file Particle.f90.

425  class(ParticleType), intent(inout) :: this !< particle
426 
427  this%xorigin = dzero
428  this%yorigin = dzero
429  this%zorigin = dzero
430  this%sinrot = dzero
431  this%cosrot = done
432  this%cosrot = done
433  this%transformed = .false.

◆ resize()

subroutine particlemodule::resize ( class(particlestoretype), intent(inout)  this,
integer(i4b), intent(in)  np,
character(*), intent(in)  mempath 
)
Parameters
[in,out]thisparticle store
[in]npnumber of particles
[in]mempathpath to memory

Definition at line 229 of file Particle.f90.

230  ! dummy
231  class(ParticleStoreType), intent(inout) :: this !< particle store
232  integer(I4B), intent(in) :: np !< number of particles
233  character(*), intent(in) :: mempath !< path to memory
234 
235  call mem_reallocate(this%imdl, np, 'PLIMDL', mempath)
236  call mem_reallocate(this%iprp, np, 'PLIPRP', mempath)
237  call mem_reallocate(this%irpt, np, 'PLIRPT', mempath)
238  call mem_reallocate(this%name, lenboundname, np, 'PLNAME', mempath)
239  call mem_reallocate(this%icu, np, 'PLICU', mempath)
240  call mem_reallocate(this%ilay, np, 'PLILAY', mempath)
241  call mem_reallocate(this%izone, np, 'PLIZONE', mempath)
242  call mem_reallocate(this%istatus, np, 'PLISTATUS', mempath)
243  call mem_reallocate(this%x, np, 'PLX', mempath)
244  call mem_reallocate(this%y, np, 'PLY', mempath)
245  call mem_reallocate(this%z, np, 'PLZ', mempath)
246  call mem_reallocate(this%trelease, np, 'PLTRELEASE', mempath)
247  call mem_reallocate(this%tstop, np, 'PLTSTOP', mempath)
248  call mem_reallocate(this%ttrack, np, 'PLTTRACK', mempath)
249  call mem_reallocate(this%istopweaksink, np, 'PLISTOPWEAKSINK', mempath)
250  call mem_reallocate(this%istopzone, np, 'PLISTOPZONE', mempath)
251  call mem_reallocate(this%idrymeth, np, 'PLIDRYMETH', mempath)
252  call mem_reallocate(this%frctrn, np, 'PLFRCTRN', mempath)
253  call mem_reallocate(this%iexmeth, np, 'PLIEXMETH', mempath)
254  call mem_reallocate(this%extol, np, 'PLEXTOL', mempath)
255  call mem_reallocate(this%extend, np, 'PLEXTEND', mempath)
256  call mem_reallocate(this%icycwin, np, 'PLICYCWIN', mempath)
257  call mem_reallocate(this%itrdomain, np, max_level, 'PLIDOMAIN', mempath)
258  call mem_reallocate(this%iboundary, np, max_level, 'PLIBOUNDARY', mempath)

◆ transform_coords()

subroutine particlemodule::transform_coords ( class(particletype), intent(inout)  this,
real(dp), intent(in), optional  xorigin,
real(dp), intent(in), optional  yorigin,
real(dp), intent(in), optional  zorigin,
real(dp), intent(in), optional  sinrot,
real(dp), intent(in), optional  cosrot,
logical(lgp), intent(in), optional  invert 
)

Apply a translation and/or rotation to particle coordinates. No rescaling. It's also possible to invert a transformation. Be sure to reset the transformation after using it.

Parameters
[in,out]thisparticle
[in]xoriginx coordinate of origin
[in]yoriginy coordinate of origin
[in]zoriginz coordinate of origin
[in]sinrotsine of rotation angle
[in]cosrotcosine of rotation angle
[in]invertwhether to invert

Definition at line 394 of file Particle.f90.

396  use geomutilmodule, only: transform, compose
397  class(ParticleType), intent(inout) :: this !< particle
398  real(DP), intent(in), optional :: xorigin !< x coordinate of origin
399  real(DP), intent(in), optional :: yorigin !< y coordinate of origin
400  real(DP), intent(in), optional :: zorigin !< z coordinate of origin
401  real(DP), intent(in), optional :: sinrot !< sine of rotation angle
402  real(DP), intent(in), optional :: cosrot !< cosine of rotation angle
403  logical(LGP), intent(in), optional :: invert !< whether to invert
404 
405  call transform(this%x, this%y, this%z, &
406  this%x, this%y, this%z, &
407  xorigin, yorigin, zorigin, &
408  sinrot, cosrot, invert)
409 
410  ! Compose is needed only because we have to untransform
411  ! coordinates: we may need to report a particle event
412  ! at any point in the tracking method hierarchy, so we
413  ! need to know how to map the particle position back to
414  ! model coords from coords local to the current domain.
415  call compose(this%xorigin, this%yorigin, this%zorigin, &
416  this%sinrot, this%cosrot, &
417  xorigin, yorigin, zorigin, &
418  sinrot, cosrot, invert)
419 
420  this%transformed = .true.
Here is the call graph for this function:

Variable Documentation

◆ max_level

integer, parameter particlemodule::max_level = 4

Definition at line 15 of file Particle.f90.

15  integer, parameter :: MAX_LEVEL = 4