MODFLOW 6  version 6.6.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 :: icp !< previous cell number (reduced)
93  integer(I4B), dimension(:), pointer, public, contiguous :: icu !< cell number (user)
94  integer(I4B), dimension(:), pointer, public, contiguous :: ilay !< layer
95  integer(I4B), dimension(:), pointer, public, contiguous :: izone !< current zone number
96  integer(I4B), dimension(:), pointer, public, contiguous :: izp !< previous zone number
97  integer(I4B), dimension(:), pointer, public, contiguous :: istatus !< particle status
98  real(dp), dimension(:), pointer, public, contiguous :: x !< model x coord of particle
99  real(dp), dimension(:), pointer, public, contiguous :: y !< model y coord of particle
100  real(dp), dimension(:), pointer, public, contiguous :: z !< model z coord of particle
101  real(dp), dimension(:), pointer, public, contiguous :: trelease !< particle release time
102  real(dp), dimension(:), pointer, public, contiguous :: tstop !< particle stop time
103  real(dp), dimension(:), pointer, public, contiguous :: ttrack !< current tracking time
104  integer(I4B), dimension(:), pointer, public, contiguous :: ifrctrn !< force ternary method
105  integer(I4B), dimension(:), pointer, public, contiguous :: iexmeth !< method for iterative solution of particle exit location and time in generalized Pollock's method
106  real(dp), dimension(:), pointer, public, contiguous :: extol !< tolerance for iterative solution of particle exit location and time in generalized Pollock's method
107  integer(LGP), dimension(:), pointer, public, contiguous :: extend !< whether to extend tracking beyond the end of the simulation
108  contains
109  procedure, public :: deallocate
110  procedure, public :: num_stored
111  procedure, public :: resize
112  procedure, public :: save_particle
113  end type particlestoretype
114 
115 contains
116 
117  !> @brief Create a new particle
118  subroutine create_particle(particle)
119  type(particletype), pointer :: particle !< particle
120  allocate (particle)
121  allocate (particle%idomain(levelmax))
122  allocate (particle%iboundary(levelmax))
123  end subroutine create_particle
124 
125  !> @brief Create a new particle store
126  subroutine allocate_particle_store(this, np, mempath)
127  type(particlestoretype), pointer :: this !< store
128  integer(I4B), intent(in) :: np !< number of particles
129  character(*), intent(in) :: mempath !< path to memory
130 
131  allocate (this)
132  call mem_allocate(this%imdl, np, 'PLIMDL', mempath)
133  call mem_allocate(this%irpt, np, 'PLIRPT', mempath)
134  call mem_allocate(this%iprp, np, 'PLIPRP', mempath)
135  call mem_allocate(this%name, lenboundname, np, 'PLNAME', mempath)
136  call mem_allocate(this%icp, np, 'PLICP', mempath)
137  call mem_allocate(this%icu, np, 'PLICU', mempath)
138  call mem_allocate(this%ilay, np, 'PLILAY', mempath)
139  call mem_allocate(this%izone, np, 'PLIZONE', mempath)
140  call mem_allocate(this%izp, np, 'PLIZP', mempath)
141  call mem_allocate(this%istatus, np, 'PLISTATUS', mempath)
142  call mem_allocate(this%x, np, 'PLX', mempath)
143  call mem_allocate(this%y, np, 'PLY', mempath)
144  call mem_allocate(this%z, np, 'PLZ', mempath)
145  call mem_allocate(this%trelease, np, 'PLTRELEASE', mempath)
146  call mem_allocate(this%tstop, np, 'PLTSTOP', mempath)
147  call mem_allocate(this%ttrack, np, 'PLTTRACK', mempath)
148  call mem_allocate(this%istopweaksink, np, 'PLISTOPWEAKSINK', mempath)
149  call mem_allocate(this%istopzone, np, 'PLISTOPZONE', mempath)
150  call mem_allocate(this%idrymeth, np, 'PLIDRYMETH', mempath)
151  call mem_allocate(this%ifrctrn, np, 'PLIFRCTRN', mempath)
152  call mem_allocate(this%iexmeth, np, 'PLIEXMETH', mempath)
153  call mem_allocate(this%extol, np, 'PLEXTOL', mempath)
154  call mem_allocate(this%extend, np, 'PLIEXTEND', mempath)
155  call mem_allocate(this%idomain, np, levelmax, 'PLIDOMAIN', mempath)
156  call mem_allocate(this%iboundary, np, levelmax, 'PLIBOUNDARY', mempath)
157  end subroutine allocate_particle_store
158 
159  !> @brief Deallocate particle arrays
160  subroutine deallocate (this, mempath)
161  class(particlestoretype), intent(inout) :: this !< store
162  character(*), intent(in) :: mempath !< path to memory
163 
164  call mem_deallocate(this%imdl, 'PLIMDL', mempath)
165  call mem_deallocate(this%iprp, 'PLIPRP', mempath)
166  call mem_deallocate(this%irpt, 'PLIRPT', mempath)
167  call mem_deallocate(this%name, 'PLNAME', mempath)
168  call mem_deallocate(this%icp, 'PLICP', mempath)
169  call mem_deallocate(this%icu, 'PLICU', mempath)
170  call mem_deallocate(this%ilay, 'PLILAY', mempath)
171  call mem_deallocate(this%izone, 'PLIZONE', mempath)
172  call mem_deallocate(this%izp, 'PLIZP', mempath)
173  call mem_deallocate(this%istatus, 'PLISTATUS', mempath)
174  call mem_deallocate(this%x, 'PLX', mempath)
175  call mem_deallocate(this%y, 'PLY', mempath)
176  call mem_deallocate(this%z, 'PLZ', mempath)
177  call mem_deallocate(this%trelease, 'PLTRELEASE', mempath)
178  call mem_deallocate(this%tstop, 'PLTSTOP', mempath)
179  call mem_deallocate(this%ttrack, 'PLTTRACK', mempath)
180  call mem_deallocate(this%istopweaksink, 'PLISTOPWEAKSINK', mempath)
181  call mem_deallocate(this%istopzone, 'PLISTOPZONE', mempath)
182  call mem_deallocate(this%idrymeth, 'PLIDRYMETH', mempath)
183  call mem_deallocate(this%ifrctrn, 'PLIFRCTRN', mempath)
184  call mem_deallocate(this%iexmeth, 'PLIEXMETH', mempath)
185  call mem_deallocate(this%extol, 'PLEXTOL', mempath)
186  call mem_deallocate(this%extend, 'PLIEXTEND', mempath)
187  call mem_deallocate(this%idomain, 'PLIDOMAIN', mempath)
188  call mem_deallocate(this%iboundary, 'PLIBOUNDARY', mempath)
189  end subroutine deallocate
190 
191  !> @brief Reallocate particle arrays
192  subroutine resize(this, np, mempath)
193  ! dummy
194  class(particlestoretype), intent(inout) :: this !< particle store
195  integer(I4B), intent(in) :: np !< number of particles
196  character(*), intent(in) :: mempath !< path to memory
197 
198  ! resize arrays
199  call mem_reallocate(this%imdl, np, 'PLIMDL', mempath)
200  call mem_reallocate(this%iprp, np, 'PLIPRP', mempath)
201  call mem_reallocate(this%irpt, np, 'PLIRPT', mempath)
202  call mem_reallocate(this%name, lenboundname, np, 'PLNAME', mempath)
203  call mem_reallocate(this%icp, np, 'PLICP', mempath)
204  call mem_reallocate(this%icu, np, 'PLICU', mempath)
205  call mem_reallocate(this%ilay, np, 'PLILAY', mempath)
206  call mem_reallocate(this%izone, np, 'PLIZONE', mempath)
207  call mem_reallocate(this%izp, np, 'PLIZP', mempath)
208  call mem_reallocate(this%istatus, np, 'PLISTATUS', mempath)
209  call mem_reallocate(this%x, np, 'PLX', mempath)
210  call mem_reallocate(this%y, np, 'PLY', mempath)
211  call mem_reallocate(this%z, np, 'PLZ', mempath)
212  call mem_reallocate(this%trelease, np, 'PLTRELEASE', mempath)
213  call mem_reallocate(this%tstop, np, 'PLTSTOP', mempath)
214  call mem_reallocate(this%ttrack, np, 'PLTTRACK', mempath)
215  call mem_reallocate(this%istopweaksink, np, 'PLISTOPWEAKSINK', mempath)
216  call mem_reallocate(this%istopzone, np, 'PLISTOPZONE', mempath)
217  call mem_reallocate(this%idrymeth, np, 'PLIDRYMETH', mempath)
218  call mem_reallocate(this%ifrctrn, np, 'PLIFRCTRN', mempath)
219  call mem_reallocate(this%iexmeth, np, 'PLIEXMETH', mempath)
220  call mem_reallocate(this%extol, np, 'PLEXTOL', mempath)
221  call mem_reallocate(this%extend, np, 'PLIEXTEND', mempath)
222  call mem_reallocate(this%idomain, np, levelmax, 'PLIDOMAIN', mempath)
223  call mem_reallocate(this%iboundary, np, levelmax, 'PLIBOUNDARY', mempath)
224  end subroutine resize
225 
226  !> @brief Load a particle from the particle store.
227  !!
228  !! This routine is used to initialize a particle for tracking.
229  !! The advancing flag and coordinate transformation are reset.
230  !<
231  subroutine load_particle(this, store, imdl, iprp, ip)
232  class(particletype), intent(inout) :: this !< particle
233  type(particlestoretype), intent(in) :: store !< particle storage
234  integer(I4B), intent(in) :: imdl !< index of model particle originated in
235  integer(I4B), intent(in) :: iprp !< index of particle release package particle originated in
236  integer(I4B), intent(in) :: ip !< index into the particle list
237 
238  call this%reset_transform()
239  this%imdl = imdl
240  this%iprp = iprp
241  this%irpt = store%irpt(ip)
242  this%ip = ip
243  this%name = store%name(ip)
244  this%istopweaksink = store%istopweaksink(ip)
245  this%istopzone = store%istopzone(ip)
246  this%idrymeth = store%idrymeth(ip)
247  this%icp = store%icp(ip)
248  this%icu = store%icu(ip)
249  this%ilay = store%ilay(ip)
250  this%izone = store%izone(ip)
251  this%izp = store%izp(ip)
252  this%istatus = store%istatus(ip)
253  this%x = store%x(ip)
254  this%y = store%y(ip)
255  this%z = store%z(ip)
256  this%trelease = store%trelease(ip)
257  this%tstop = store%tstop(ip)
258  this%ttrack = store%ttrack(ip)
259  this%advancing = .true.
260  this%idomain(1:levelmax) = &
261  store%idomain(ip, 1:levelmax)
262  this%idomain(1) = imdl
263  this%iboundary(1:levelmax) = &
264  store%iboundary(ip, 1:levelmax)
265  this%ifrctrn = store%ifrctrn(ip)
266  this%iexmeth = store%iexmeth(ip)
267  this%extol = store%extol(ip)
268  this%iextend = store%extend(ip)
269  end subroutine load_particle
270 
271  !> @brief Save a particle's state to the particle store
272  subroutine save_particle(this, particle, ip)
273  class(particlestoretype), intent(inout) :: this !< particle storage
274  type(particletype), intent(in) :: particle !< particle
275  integer(I4B), intent(in) :: ip !< particle index
276 
277  this%imdl(ip) = particle%imdl
278  this%iprp(ip) = particle%iprp
279  this%irpt(ip) = particle%irpt
280  this%name(ip) = particle%name
281  this%istopweaksink(ip) = particle%istopweaksink
282  this%istopzone(ip) = particle%istopzone
283  this%idrymeth(ip) = particle%idrymeth
284  this%icp(ip) = particle%icp
285  this%icu(ip) = particle%icu
286  this%ilay(ip) = particle%ilay
287  this%izone(ip) = particle%izone
288  this%izp(ip) = particle%izp
289  this%istatus(ip) = particle%istatus
290  this%x(ip) = particle%x
291  this%y(ip) = particle%y
292  this%z(ip) = particle%z
293  this%trelease(ip) = particle%trelease
294  this%tstop(ip) = particle%tstop
295  this%ttrack(ip) = particle%ttrack
296  this%idomain( &
297  ip, &
298  1:levelmax) = &
299  particle%idomain(1:levelmax)
300  this%iboundary( &
301  ip, &
302  1:levelmax) = &
303  particle%iboundary(1:levelmax)
304  this%ifrctrn(ip) = particle%ifrctrn
305  this%iexmeth(ip) = particle%iexmeth
306  this%extol(ip) = particle%extol
307  this%extend(ip) = particle%iextend
308  end subroutine save_particle
309 
310  !> @brief Transform particle coordinates.
311  !!
312  !! Apply a translation and/or rotation to particle coordinates.
313  !! No rescaling. It's also possible to invert a transformation.
314  !! Be sure to reset the transformation after using it.
315  !<
316  subroutine transform_coords(this, xorigin, yorigin, zorigin, &
317  sinrot, cosrot, invert)
318  use geomutilmodule, only: transform, compose
319  class(particletype), intent(inout) :: this !< particle
320  real(DP), intent(in), optional :: xorigin !< x coordinate of origin
321  real(DP), intent(in), optional :: yorigin !< y coordinate of origin
322  real(DP), intent(in), optional :: zorigin !< z coordinate of origin
323  real(DP), intent(in), optional :: sinrot !< sine of rotation angle
324  real(DP), intent(in), optional :: cosrot !< cosine of rotation angle
325  logical(LGP), intent(in), optional :: invert !< whether to invert
326 
327  call transform(this%x, this%y, this%z, &
328  this%x, this%y, this%z, &
329  xorigin, yorigin, zorigin, &
330  sinrot, cosrot, invert)
331 
332  call compose(this%xorigin, this%yorigin, this%zorigin, &
333  this%sinrot, this%cosrot, &
334  xorigin, yorigin, zorigin, &
335  sinrot, cosrot, invert)
336 
337  this%transformed = .true.
338  end subroutine transform_coords
339 
340  subroutine reset_transform(this)
341  class(particletype), intent(inout) :: this !< particle
342 
343  this%xorigin = dzero
344  this%yorigin = dzero
345  this%zorigin = dzero
346  this%sinrot = dzero
347  this%cosrot = done
348  this%cosrot = done
349  this%transformed = .false.
350  end subroutine reset_transform
351 
352  !> @brief Return the particle's model coordinates.
353  subroutine get_model_coords(this, x, y, z)
354  use geomutilmodule, only: transform, compose
355  class(particletype), intent(inout) :: this !< particle
356  real(DP), intent(out) :: x !< x coordinate
357  real(DP), intent(out) :: y !< y coordinate
358  real(DP), intent(out) :: z !< z coordinate
359 
360  if (this%transformed) then
361  ! Untransform coordinates
362  call transform(this%x, this%y, this%z, x, y, z, &
363  this%xorigin, this%yorigin, this%zorigin, &
364  this%sinrot, this%cosrot, invert=.true.)
365  else
366  x = this%x
367  y = this%y
368  z = this%z
369  end if
370  end subroutine get_model_coords
371 
372  integer function num_stored(this) result(n)
373  class(particlestoretype) :: this
374  n = size(this%imdl)
375  end function num_stored
376 
377 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:232
particlemodule::allocate_particle_store
subroutine, public allocate_particle_store(this, np, mempath)
Create a new particle store.
Definition: Particle.f90:127
particlemodule::resize
subroutine resize(this, np, mempath)
Reallocate particle arrays.
Definition: Particle.f90:193
particlemodule::get_model_coords
subroutine get_model_coords(this, x, y, z)
Return the particle's model coordinates.
Definition: Particle.f90:354
particlemodule::num_stored
integer function num_stored(this)
Definition: Particle.f90:373
particlemodule::save_particle
subroutine save_particle(this, particle, ip)
Save a particle's state to the particle store.
Definition: Particle.f90:273
particlemodule::reset_transform
subroutine reset_transform(this)
Definition: Particle.f90:341
particlemodule::transform_coords
subroutine transform_coords(this, xorigin, yorigin, zorigin, sinrot, cosrot, invert)
Transform particle coordinates.
Definition: Particle.f90:318
particlemodule::levelmax
integer, parameter, public levelmax
Definition: Particle.f90:14
particlemodule::deallocate
subroutine deallocate(this, mempath)
Deallocate particle arrays.
Definition: Particle.f90:161
particlemodule::create_particle
subroutine, public create_particle(particle)
Create a new particle.
Definition: Particle.f90:119
particlemodule::particlestoretype
Structure of arrays to store particles.
Definition: Particle.f90:78
particlemodule::particletype
Particle tracked by the PRT model.
Definition: Particle.f90:32