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