MODFLOW 6  version 6.7.0.dev3
USGS Modular Hydrologic Model
Particle.f90
Go to the documentation of this file.
2 
3  use kindmodule, only: dp, i4b, lgp
4  use listmodule, only: listtype
8  implicit none
9  public
10 
11  !> Tracking "levels" defined in method modules. Currently only 3 used.
12  integer, parameter :: max_level = 4
13 
14  !> @brief Particle status enumeration.
15  !!
16  !! Particles begin in status 1 (active) at release time. Status may only
17  !! increase over time. Status values greater than one imply termination.
18  !! A particle may terminate for several reasons, all mutually exclusive.
19  !! A particle's final tracking status will always be greater than one.
20  !!
21  !! Status codes 0-3 and 5-8 correspond directly to MODPATH 7 status codes.
22  !! Code 4 does not apply to PRT because PRT does not distinguish forwards
23  !! from backwards tracking. Status code 9 provides more specific, subcell-
24  !! level information about a particle which terminates due to no outflow.
25  !! Code 10 distinguishes particles which have "timed out" upon reaching a
26  !! user-specified stop time or the end of the simulation.
27  !<
28  enum, bind(C)
29  enumerator :: active = 1
30  enumerator :: term_boundary = 2 !< terminated at a boundary face
31  enumerator :: term_weaksink = 3 !< terminated in a weak sink cell
32  enumerator :: term_no_exits = 5 !< terminated in a cell with no exit face
33  enumerator :: term_stopzone = 6 !< terminated in a cell with a stop zone number
34  enumerator :: term_inactive = 7 !< terminated in an inactive cell
35  enumerator :: term_unreleased = 8 !< terminated permanently unreleased
36  enumerator :: term_no_exits_sub = 9 !< terminated in a subcell with no exit face
37  enumerator :: term_timeout = 10 !< terminated at stop time or end of simulation
38  end enum
39 
40  !> @brief Particle tracked by the PRT model.
41  !!
42  !! Record-type to conveniently shuffle a particle's
43  !! state to/from storage before/after its trajectory
44  !! is solved for each time step.
45  !!
46  !! Particle coordinates may be local to the cell or
47  !! global/model. Routines are provided to convert a
48  !! particle's global coordinates to/from cell-local
49  !! coordinates for tracking through cell subdomains.
50  !!
51  !! Particles are identified by composite key, i.e.,
52  !! combinations of properties imdl, iprp, irpt, and
53  !! trelease. An optional label may be provided, but
54  !! need not be unique
55  !<
57  private
58  ! identity
59  character(len=LENBOUNDNAME), public :: name = '' !< optional particle name
60  integer(I4B), public :: imdl !< index of model the particle originated in
61  integer(I4B), public :: iprp !< index of release package the particle is from
62  integer(I4B), public :: irpt !< index of release point the particle is from
63  integer(I4B), public :: ip !< index of particle in the particle list
64  ! options
65  logical(LGP), public :: extend !< whether to extend tracking beyond the end of the simulation
66  logical(LGP), public :: frctrn !< whether to force solving the particle with the ternary method
67  integer(I4B), public :: istopweaksink !< weak sink option (0: do not stop, 1: stop)
68  integer(I4B), public :: istopzone !< stop zone number
69  integer(I4B), public :: idrymeth !< dry tracking method
70  integer(I4B), public :: iexmeth !< method for iterative solution of particle exit location and time in generalized Pollock's method
71  integer(I4B), public :: icycwin !< cycle detection window size
72  real(dp), public :: extol !< tolerance for iterative solution of particle exit location and time in generalized Pollock's method
73  ! state
74  integer(I4B), public :: itrdomain(max_level) !< tracking domain indices
75  integer(I4B), public :: iboundary(max_level) !< tracking domain boundary indices
76  integer(I4B), public :: icu !< user cell number
77  integer(I4B), public :: ilay !< grid layer
78  integer(I4B), public :: izone !< current zone number
79  integer(I4B), public :: istatus !< tracking status
80  real(dp), public :: x !< x coordinate
81  real(dp), public :: y !< y coordinate
82  real(dp), public :: z !< z coordinate
83  real(dp), public :: trelease !< release time
84  real(dp), public :: tstop !< stop time
85  real(dp), public :: ttrack !< time tracked so far
86  real(dp), public :: xorigin !< x origin for coordinate transformation from model to local
87  real(dp), public :: yorigin !< y origin for coordinate transformation from model to local
88  real(dp), public :: zorigin !< z origin for coordinate transformation from model to local
89  real(dp), public :: sinrot !< sine of rotation angle for coordinate transformation from model to local
90  real(dp), public :: cosrot !< cosine of rotation angle for coordinate transformation from model to local
91  logical(LGP), public :: transformed !< whether coordinates have been transformed from model to local
92  logical(LGP), public :: advancing !< whether particle is still being tracked for current time step
93  type(listtype), public, pointer :: history !< history of particle positions (for cycle detection)
94  contains
95  procedure, public :: destroy => destroy_particle
96  procedure, public :: get_model_coords
97  procedure, public :: transform => transform_coords
98  procedure, public :: reset_transform
99  end type particletype
100 
101  !> @brief Structure of arrays to store particles.
103  private
104  ! identity
105  character(len=LENBOUNDNAME), dimension(:), pointer, public, contiguous :: name !< optional particle label
106  integer(I4B), dimension(:), pointer, public, contiguous :: imdl !< index of model particle originated in
107  integer(I4B), dimension(:), pointer, public, contiguous :: iprp !< index of release package the particle originated in
108  integer(I4B), dimension(:), pointer, public, contiguous :: irpt !< index of release point in the particle release package the particle originated in
109  ! options
110  logical(LGP), dimension(:), pointer, public, contiguous :: extend !< whether to extend tracking beyond the end of the simulation
111  logical(LGP), dimension(:), pointer, public, contiguous :: frctrn !< force ternary method
112  integer(I4B), dimension(:), pointer, public, contiguous :: istopweaksink !< weak sink option: 0 = do not stop, 1 = stop
113  integer(I4B), dimension(:), pointer, public, contiguous :: istopzone !< stop zone number
114  integer(I4B), dimension(:), pointer, public, contiguous :: idrymeth !< stop in dry cells
115  integer(I4B), dimension(:), pointer, public, contiguous :: iexmeth !< method for iterative solution of particle exit location and time in generalized Pollock's method
116  integer(I4B), dimension(:), pointer, public, contiguous :: icycwin !< cycle detection window size
117  real(dp), dimension(:), pointer, public, contiguous :: extol !< tolerance for iterative solution of particle exit location and time in generalized Pollock's method
118  ! state
119  integer(I4B), dimension(:, :), pointer, public, contiguous :: itrdomain !< array of indices for domains in the tracking domain hierarchy
120  integer(I4B), dimension(:, :), pointer, public, contiguous :: iboundary !< array of indices for tracking domain boundaries
121  integer(I4B), dimension(:), pointer, public, contiguous :: icu !< cell number (user)
122  integer(I4B), dimension(:), pointer, public, contiguous :: ilay !< layer
123  integer(I4B), dimension(:), pointer, public, contiguous :: izone !< current zone number
124  integer(I4B), dimension(:), pointer, public, contiguous :: istatus !< particle status
125  real(dp), dimension(:), pointer, public, contiguous :: x !< model x coord of particle
126  real(dp), dimension(:), pointer, public, contiguous :: y !< model y coord of particle
127  real(dp), dimension(:), pointer, public, contiguous :: z !< model z coord of particle
128  real(dp), dimension(:), pointer, public, contiguous :: trelease !< particle release time
129  real(dp), dimension(:), pointer, public, contiguous :: tstop !< particle stop time
130  real(dp), dimension(:), pointer, public, contiguous :: ttrack !< current tracking time
131  contains
132  procedure, public :: destroy
133  procedure, public :: num_stored
134  procedure, public :: resize
135  procedure, public :: get
136  procedure, public :: put
137  end type particlestoretype
138 
139 contains
140 
141  !> @brief Create a new particle
142  subroutine create_particle(particle)
143  type(particletype), pointer :: particle !< particle
144  allocate (particle)
145  allocate (particle%history)
146  end subroutine create_particle
147 
148  !> @brief Allocate particle store
149  subroutine create_particle_store(store, np, mempath)
150  type(particlestoretype), pointer :: store !< store
151  integer(I4B), intent(in) :: np !< number of particles
152  character(*), intent(in) :: mempath !< path to memory
153 
154  allocate (store)
155  call mem_allocate(store%imdl, np, 'PLIMDL', mempath)
156  call mem_allocate(store%irpt, np, 'PLIRPT', mempath)
157  call mem_allocate(store%iprp, np, 'PLIPRP', mempath)
158  call mem_allocate(store%name, lenboundname, np, 'PLNAME', mempath)
159  call mem_allocate(store%icu, np, 'PLICU', mempath)
160  call mem_allocate(store%ilay, np, 'PLILAY', mempath)
161  call mem_allocate(store%izone, np, 'PLIZONE', mempath)
162  call mem_allocate(store%istatus, np, 'PLISTATUS', mempath)
163  call mem_allocate(store%x, np, 'PLX', mempath)
164  call mem_allocate(store%y, np, 'PLY', mempath)
165  call mem_allocate(store%z, np, 'PLZ', mempath)
166  call mem_allocate(store%trelease, np, 'PLTRELEASE', mempath)
167  call mem_allocate(store%tstop, np, 'PLTSTOP', mempath)
168  call mem_allocate(store%ttrack, np, 'PLTTRACK', mempath)
169  call mem_allocate(store%istopweaksink, np, 'PLISTOPWEAKSINK', mempath)
170  call mem_allocate(store%istopzone, np, 'PLISTOPZONE', mempath)
171  call mem_allocate(store%idrymeth, np, 'PLIDRYMETH', mempath)
172  call mem_allocate(store%frctrn, np, 'PLFRCTRN', mempath)
173  call mem_allocate(store%iexmeth, np, 'PLIEXMETH', mempath)
174  call mem_allocate(store%extol, np, 'PLEXTOL', mempath)
175  call mem_allocate(store%extend, np, 'PLEXTEND', mempath)
176  call mem_allocate(store%icycwin, np, 'PLICYCWIN', mempath)
177  call mem_allocate(store%itrdomain, np, max_level, 'PLIDOMAIN', mempath)
178  call mem_allocate(store%iboundary, np, max_level, 'PLIBOUNDARY', mempath)
179  end subroutine create_particle_store
180 
181  !> @brief Destroy particle store after use.
182  subroutine destroy(this, mempath)
183  class(particlestoretype), intent(inout) :: this !< store
184  character(*), intent(in) :: mempath !< path to memory
185 
186  call mem_deallocate(this%imdl, 'PLIMDL', mempath)
187  call mem_deallocate(this%iprp, 'PLIPRP', mempath)
188  call mem_deallocate(this%irpt, 'PLIRPT', mempath)
189  call mem_deallocate(this%name, 'PLNAME', mempath)
190  call mem_deallocate(this%icu, 'PLICU', mempath)
191  call mem_deallocate(this%ilay, 'PLILAY', mempath)
192  call mem_deallocate(this%izone, 'PLIZONE', mempath)
193  call mem_deallocate(this%istatus, 'PLISTATUS', mempath)
194  call mem_deallocate(this%x, 'PLX', mempath)
195  call mem_deallocate(this%y, 'PLY', mempath)
196  call mem_deallocate(this%z, 'PLZ', mempath)
197  call mem_deallocate(this%trelease, 'PLTRELEASE', mempath)
198  call mem_deallocate(this%tstop, 'PLTSTOP', mempath)
199  call mem_deallocate(this%ttrack, 'PLTTRACK', mempath)
200  call mem_deallocate(this%istopweaksink, 'PLISTOPWEAKSINK', mempath)
201  call mem_deallocate(this%istopzone, 'PLISTOPZONE', mempath)
202  call mem_deallocate(this%idrymeth, 'PLIDRYMETH', mempath)
203  call mem_deallocate(this%frctrn, 'PLFRCTRN', mempath)
204  call mem_deallocate(this%iexmeth, 'PLIEXMETH', mempath)
205  call mem_deallocate(this%extol, 'PLEXTOL', mempath)
206  call mem_deallocate(this%extend, 'PLEXTEND', mempath)
207  call mem_deallocate(this%icycwin, 'PLICYCWIN', mempath)
208  call mem_deallocate(this%itrdomain, 'PLIDOMAIN', mempath)
209  call mem_deallocate(this%iboundary, 'PLIBOUNDARY', mempath)
210  end subroutine destroy
211 
212  !> @brief Destroy a particle after use.
213  subroutine destroy_particle(particle)
214  class(particletype), intent(inout) :: particle !< particle
215  deallocate (particle%history)
216  end subroutine destroy_particle
217 
218  !> @brief Reallocate particle storage to the given size.
219  subroutine resize(this, np, mempath)
220  ! dummy
221  class(particlestoretype), intent(inout) :: this !< particle store
222  integer(I4B), intent(in) :: np !< number of particles
223  character(*), intent(in) :: mempath !< path to memory
224 
225  ! resize arrays
226  call mem_reallocate(this%imdl, np, 'PLIMDL', mempath)
227  call mem_reallocate(this%iprp, np, 'PLIPRP', mempath)
228  call mem_reallocate(this%irpt, np, 'PLIRPT', mempath)
229  call mem_reallocate(this%name, lenboundname, np, 'PLNAME', mempath)
230  call mem_reallocate(this%icu, np, 'PLICU', mempath)
231  call mem_reallocate(this%ilay, np, 'PLILAY', mempath)
232  call mem_reallocate(this%izone, np, 'PLIZONE', mempath)
233  call mem_reallocate(this%istatus, np, 'PLISTATUS', mempath)
234  call mem_reallocate(this%x, np, 'PLX', mempath)
235  call mem_reallocate(this%y, np, 'PLY', mempath)
236  call mem_reallocate(this%z, np, 'PLZ', mempath)
237  call mem_reallocate(this%trelease, np, 'PLTRELEASE', mempath)
238  call mem_reallocate(this%tstop, np, 'PLTSTOP', mempath)
239  call mem_reallocate(this%ttrack, np, 'PLTTRACK', mempath)
240  call mem_reallocate(this%istopweaksink, np, 'PLISTOPWEAKSINK', mempath)
241  call mem_reallocate(this%istopzone, np, 'PLISTOPZONE', mempath)
242  call mem_reallocate(this%idrymeth, np, 'PLIDRYMETH', mempath)
243  call mem_reallocate(this%frctrn, np, 'PLFRCTRN', mempath)
244  call mem_reallocate(this%iexmeth, np, 'PLIEXMETH', mempath)
245  call mem_reallocate(this%extol, np, 'PLEXTOL', mempath)
246  call mem_reallocate(this%extend, np, 'PLEXTEND', mempath)
247  call mem_reallocate(this%icycwin, np, 'PLICYCWIN', mempath)
248  call mem_reallocate(this%itrdomain, np, max_level, 'PLIDOMAIN', mempath)
249  call mem_reallocate(this%iboundary, np, max_level, 'PLIBOUNDARY', mempath)
250  end subroutine resize
251 
252  !> @brief Load a particle from the particle store.
253  !!
254  !! This routine is used to initialize a particle for tracking.
255  !! The advancing flag and coordinate transformation are reset.
256  !<
257  subroutine get(this, particle, imdl, iprp, ip)
258  class(particlestoretype), intent(inout) :: this !< particle store
259  class(particletype), intent(inout) :: particle !< particle
260  integer(I4B), intent(in) :: imdl !< index of model particle originated in
261  integer(I4B), intent(in) :: iprp !< index of particle release package particle originated in
262  integer(I4B), intent(in) :: ip !< index into the particle list
263 
264  call particle%reset_transform()
265  call particle%history%Clear()
266  particle%imdl = imdl
267  particle%iprp = iprp
268  particle%irpt = this%irpt(ip)
269  particle%ip = ip
270  particle%name = this%name(ip)
271  particle%istopweaksink = this%istopweaksink(ip)
272  particle%istopzone = this%istopzone(ip)
273  particle%idrymeth = this%idrymeth(ip)
274  particle%icu = this%icu(ip)
275  particle%ilay = this%ilay(ip)
276  particle%izone = this%izone(ip)
277  particle%istatus = this%istatus(ip)
278  particle%x = this%x(ip)
279  particle%y = this%y(ip)
280  particle%z = this%z(ip)
281  particle%trelease = this%trelease(ip)
282  particle%tstop = this%tstop(ip)
283  particle%ttrack = this%ttrack(ip)
284  particle%advancing = .true.
285  particle%itrdomain(1:max_level) = &
286  this%itrdomain(ip, 1:max_level)
287  particle%itrdomain(1) = imdl
288  particle%iboundary(1:max_level) = &
289  this%iboundary(ip, 1:max_level)
290  particle%frctrn = this%frctrn(ip)
291  particle%iexmeth = this%iexmeth(ip)
292  particle%extol = this%extol(ip)
293  particle%extend = this%extend(ip)
294  particle%icycwin = this%icycwin(ip)
295  end subroutine get
296 
297  !> @brief Save a particle's state to the particle store.
298  subroutine put(this, particle, ip)
299  class(particlestoretype), intent(inout) :: this !< particle storage
300  class(particletype), intent(in) :: particle !< particle
301  integer(I4B), intent(in) :: ip !< particle index
302 
303  this%imdl(ip) = particle%imdl
304  this%iprp(ip) = particle%iprp
305  this%irpt(ip) = particle%irpt
306  this%name(ip) = particle%name
307  this%istopweaksink(ip) = particle%istopweaksink
308  this%istopzone(ip) = particle%istopzone
309  this%idrymeth(ip) = particle%idrymeth
310  this%icu(ip) = particle%icu
311  this%ilay(ip) = particle%ilay
312  this%izone(ip) = particle%izone
313  this%istatus(ip) = particle%istatus
314  this%x(ip) = particle%x
315  this%y(ip) = particle%y
316  this%z(ip) = particle%z
317  this%trelease(ip) = particle%trelease
318  this%tstop(ip) = particle%tstop
319  this%ttrack(ip) = particle%ttrack
320  this%itrdomain( &
321  ip, &
322  1:max_level) = &
323  particle%itrdomain(1:max_level)
324  this%iboundary( &
325  ip, &
326  1:max_level) = &
327  particle%iboundary(1:max_level)
328  this%frctrn(ip) = particle%frctrn
329  this%iexmeth(ip) = particle%iexmeth
330  this%extol(ip) = particle%extol
331  this%extend(ip) = particle%extend
332  this%icycwin(ip) = particle%icycwin
333  end subroutine put
334 
335  !> @brief Transform particle coordinates.
336  !!
337  !! Apply a translation and/or rotation to particle coordinates.
338  !! No rescaling. It's also possible to invert a transformation.
339  !! Be sure to reset the transformation after using it.
340  !<
341  subroutine transform_coords(this, xorigin, yorigin, zorigin, &
342  sinrot, cosrot, invert)
343  use geomutilmodule, only: transform, compose
344  class(particletype), intent(inout) :: this !< particle
345  real(DP), intent(in), optional :: xorigin !< x coordinate of origin
346  real(DP), intent(in), optional :: yorigin !< y coordinate of origin
347  real(DP), intent(in), optional :: zorigin !< z coordinate of origin
348  real(DP), intent(in), optional :: sinrot !< sine of rotation angle
349  real(DP), intent(in), optional :: cosrot !< cosine of rotation angle
350  logical(LGP), intent(in), optional :: invert !< whether to invert
351 
352  call transform(this%x, this%y, this%z, &
353  this%x, this%y, this%z, &
354  xorigin, yorigin, zorigin, &
355  sinrot, cosrot, invert)
356 
357  call compose(this%xorigin, this%yorigin, this%zorigin, &
358  this%sinrot, this%cosrot, &
359  xorigin, yorigin, zorigin, &
360  sinrot, cosrot, invert)
361 
362  this%transformed = .true.
363  end subroutine transform_coords
364 
365  !> @brief Reset particle coordinate transformation properties.
366  subroutine reset_transform(this)
367  class(particletype), intent(inout) :: this !< particle
368 
369  this%xorigin = dzero
370  this%yorigin = dzero
371  this%zorigin = dzero
372  this%sinrot = dzero
373  this%cosrot = done
374  this%cosrot = done
375  this%transformed = .false.
376  end subroutine reset_transform
377 
378  !> @brief Return the particle's model coordinates,
379  !! inverting any applied transformation if needed.
380  !! The particle's state is not altered.
381  subroutine get_model_coords(this, x, y, z)
382  use geomutilmodule, only: transform, compose
383  class(particletype), intent(in) :: this !< particle
384  real(DP), intent(out) :: x !< x coordinate
385  real(DP), intent(out) :: y !< y coordinate
386  real(DP), intent(out) :: z !< z coordinate
387 
388  if (this%transformed) then
389  call transform(this%x, this%y, this%z, x, y, z, &
390  this%xorigin, this%yorigin, this%zorigin, &
391  this%sinrot, this%cosrot, invert=.true.)
392  else
393  x = this%x
394  y = this%y
395  z = this%z
396  end if
397  end subroutine get_model_coords
398 
399  !> @brief Return the number of particles.
400  integer function num_stored(this) result(n)
401  class(particlestoretype) :: this
402  n = size(this%imdl)
403  end function num_stored
404 
405 end module particlemodule
This module contains simulation constants.
Definition: Constants.f90:9
integer(i4b), parameter lenboundname
maximum length of a bound name
Definition: Constants.f90:36
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
integer(i4b), parameter lenmempath
maximum length of the memory path
Definition: Constants.f90:27
real(dp), parameter done
real constant 1
Definition: Constants.f90:76
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
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
This module defines variable data types.
Definition: kind.f90:8
subroutine get(this, particle, imdl, iprp, ip)
Load a particle from the particle store.
Definition: Particle.f90:258
subroutine resize(this, np, mempath)
Reallocate particle storage to the given size.
Definition: Particle.f90:220
subroutine create_particle_store(store, np, mempath)
Allocate particle store.
Definition: Particle.f90:150
subroutine get_model_coords(this, x, y, z)
Return the particle's model coordinates, inverting any applied transformation if needed....
Definition: Particle.f90:382
integer function num_stored(this)
Return the number of particles.
Definition: Particle.f90:401
subroutine reset_transform(this)
Reset particle coordinate transformation properties.
Definition: Particle.f90:367
@ term_weaksink
terminated in a weak sink cell
Definition: Particle.f90:31
@ term_timeout
terminated at stop time or end of simulation
Definition: Particle.f90:37
@ term_inactive
terminated in an inactive cell
Definition: Particle.f90:34
@ term_no_exits
terminated in a cell with no exit face
Definition: Particle.f90:32
@ term_stopzone
terminated in a cell with a stop zone number
Definition: Particle.f90:33
@ term_no_exits_sub
terminated in a subcell with no exit face
Definition: Particle.f90:36
@ term_unreleased
terminated permanently unreleased
Definition: Particle.f90:35
@ term_boundary
terminated at a boundary face
Definition: Particle.f90:30
subroutine transform_coords(this, xorigin, yorigin, zorigin, sinrot, cosrot, invert)
Transform particle coordinates.
Definition: Particle.f90:343
subroutine destroy(this, mempath)
Destroy particle store after use.
Definition: Particle.f90:183
subroutine create_particle(particle)
Create a new particle.
Definition: Particle.f90:143
subroutine put(this, particle, ip)
Save a particle's state to the particle store.
Definition: Particle.f90:299
integer, parameter max_level
Tracking "levels" defined in method modules. Currently only 3 used.
Definition: Particle.f90:12
subroutine destroy_particle(particle)
Destroy a particle after use.
Definition: Particle.f90:214
A generic heterogeneous doubly-linked list.
Definition: List.f90:14
Structure of arrays to store particles.
Definition: Particle.f90:102
Particle tracked by the PRT model.
Definition: Particle.f90:56