MODFLOW 6  version 6.8.0.dev0
USGS Modular Hydrologic Model
methoddismodule Module Reference

Data Types

type  methoddistype
 

Functions/Subroutines

subroutine, public create_method_dis (method)
 Create a new structured grid (DIS) tracking method. More...
 
subroutine deallocate (this)
 Destructor the tracking method. More...
 
subroutine load_cell (this, ic, cell)
 
subroutine load_dis (this, particle, next_level, submethod)
 Load the cell geometry and tracking method. More...
 
subroutine load_particle (this, cell, particle)
 Load cell properties into the particle, including. More...
 
subroutine update_flowja (this, cell, particle)
 Update cell-cell flows of particle mass. Every particle is currently assigned unit mass. More...
 
subroutine pass_dis (this, particle)
 Pass a particle to the next cell, if there is one. More...
 
subroutine apply_dis (this, particle, tmax)
 Apply the structured tracking method to a particle. More...
 
double precision function get_top (this, iatop)
 Returns a top elevation based on index iatop. More...
 
subroutine load_cell_defn (this, ic, defn)
 Loads cell definition from the grid. More...
 
subroutine load_cell_properties (this, ic, defn)
 
subroutine load_cell_neighbors (this, defn)
 Loads face neighbors to cell definition from the grid. Assumes cell index and number of vertices are already loaded. More...
 
subroutine load_cell_flows (this, defn)
 Load flows into the cell definition. These include face, boundary and net distributed flows. Assumes cell index and number of vertices are already loaded. More...
 
subroutine load_cell_face_flows (this, defn)
 
subroutine load_cell_boundary_flows (this, defn)
 Add boundary flows to the cell definition faceflow array. Assumes cell index and number of vertices are already loaded. More...
 

Function/Subroutine Documentation

◆ apply_dis()

subroutine methoddismodule::apply_dis ( class(methoddistype), intent(inout)  this,
type(particletype), intent(inout), pointer  particle,
real(dp), intent(in)  tmax 
)

Definition at line 314 of file MethodDis.f90.

315  class(MethodDisType), intent(inout) :: this
316  type(ParticleType), pointer, intent(inout) :: particle
317  real(DP), intent(in) :: tmax
318 
319  call this%track(particle, 1, tmax)

◆ create_method_dis()

subroutine, public methoddismodule::create_method_dis ( type(methoddistype), pointer  method)

Definition at line 49 of file MethodDis.f90.

50  ! dummy
51  type(MethodDisType), pointer :: method
52  ! local
53  type(CellRectType), pointer :: cell
54 
55  allocate (method)
56  allocate (method%name)
57  call create_cell_rect(cell)
58  method%cell => cell
59  method%name = "dis"
60  method%delegates = .true.
61 
62  ! Create cell methods this method can delegate to
63  call create_method_cell_pollock(method%method_cell_plck)
64  call create_method_cell_ptb(method%method_cell_ptb)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ deallocate()

subroutine methoddismodule::deallocate ( class(methoddistype), intent(inout)  this)
private

Definition at line 68 of file MethodDis.f90.

69  class(MethodDisType), intent(inout) :: this
70  deallocate (this%name)
71 
72  ! Deallocate owned cell methods
73  call this%method_cell_plck%deallocate()
74  deallocate (this%method_cell_plck)
75  call this%method_cell_ptb%deallocate()
76  deallocate (this%method_cell_ptb)

◆ get_top()

double precision function methoddismodule::get_top ( class(methoddistype), intent(inout)  this,
integer, intent(in)  iatop 
)
private

Definition at line 323 of file MethodDis.f90.

324  class(MethodDisType), intent(inout) :: this
325  integer, intent(in) :: iatop
326  double precision :: top
327 
328  if (iatop .lt. 0) then
329  top = this%fmi%dis%top(-iatop)
330  else
331  top = this%fmi%dis%bot(iatop)
332  end if

◆ load_cell()

subroutine methoddismodule::load_cell ( class(methoddistype), intent(inout)  this,
integer(i4b), intent(in)  ic,
type(cellrecttype), intent(inout)  cell 
)
private

Definition at line 79 of file MethodDis.f90.

80  ! dummy
81  class(MethodDisType), intent(inout) :: this
82  integer(I4B), intent(in) :: ic
83  type(CellRectType), intent(inout) :: cell
84  ! local
85  integer(I4B) :: icu
86  integer(I4B) :: irow
87  integer(I4B) :: jcol
88  integer(I4B) :: klay
89  real(DP) :: areax
90  real(DP) :: areay
91  real(DP) :: areaz
92  real(DP) :: dx
93  real(DP) :: dy
94  real(DP) :: dz
95  real(DP) :: factor
96  real(DP) :: term
97 
98  select type (dis => this%fmi%dis)
99  type is (distype)
100  icu = dis%get_nodeuser(ic)
101 
102  call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, &
103  irow, jcol, klay)
104 
105  dx = dis%delr(jcol)
106  dy = dis%delc(irow)
107  dz = cell%defn%top - cell%defn%bot
108 
109  cell%dx = dx
110  cell%dy = dy
111  cell%dz = dz
112  cell%sinrot = dzero
113  cell%cosrot = done
114  cell%xOrigin = cell%defn%polyvert(1, 1)
115  cell%yOrigin = cell%defn%polyvert(2, 1)
116  cell%zOrigin = cell%defn%bot
117  cell%ipvOrigin = 1
118 
119  factor = done / cell%defn%retfactor
120  factor = factor / cell%defn%porosity
121 
122  areaz = dx * dy
123  term = factor / areaz
124 
125  cell%vz1 = cell%defn%faceflow(6) * term
126  cell%vz2 = -cell%defn%faceflow(7) * term
127 
128  if (this%fmi%ibdgwfsat0(ic) == 0) then
129  cell%vx1 = dzero
130  cell%vx2 = dzero
131  cell%vy1 = dzero
132  cell%vy2 = dzero
133  cell%vz1 = dzero
134  cell%vz2 = dzero
135  return
136  end if
137 
138  areax = dy * dz
139  term = factor / areax
140  cell%vx1 = cell%defn%faceflow(1) * term
141  cell%vx2 = -cell%defn%faceflow(3) * term
142 
143  areay = dx * dz
144  term = factor / areay
145  cell%vy1 = cell%defn%faceflow(4) * term
146  cell%vy2 = -cell%defn%faceflow(2) * term
147 
148  end select
Here is the call graph for this function:

◆ load_cell_boundary_flows()

subroutine methoddismodule::load_cell_boundary_flows ( class(methoddistype), intent(inout)  this,
type(celldefntype), intent(inout), pointer  defn 
)
private

Definition at line 503 of file MethodDis.f90.

504  ! dummy
505  class(MethodDisType), intent(inout) :: this
506  type(CellDefnType), pointer, intent(inout) :: defn
507  ! local
508 
509  defn%faceflow(1) = defn%faceflow(1) + &
510  this%fmi%BoundaryFlows(defn%icell, 1)
511  defn%faceflow(2) = defn%faceflow(2) + &
512  this%fmi%BoundaryFlows(defn%icell, 2)
513  defn%faceflow(3) = defn%faceflow(3) + &
514  this%fmi%BoundaryFlows(defn%icell, 3)
515  defn%faceflow(4) = defn%faceflow(4) + &
516  this%fmi%BoundaryFlows(defn%icell, 4)
517  defn%faceflow(5) = defn%faceflow(1)
518  defn%faceflow(6) = defn%faceflow(6) + &
519  this%fmi%BoundaryFlows(defn%icell, this%fmi%max_faces - 1)
520  defn%faceflow(7) = defn%faceflow(7) + &
521  this%fmi%BoundaryFlows(defn%icell, this%fmi%max_faces)

◆ load_cell_defn()

subroutine methoddismodule::load_cell_defn ( class(methoddistype), intent(inout)  this,
integer(i4b), intent(in)  ic,
type(celldefntype), intent(inout), pointer  defn 
)
private

Definition at line 336 of file MethodDis.f90.

337  ! dummy
338  class(MethodDisType), intent(inout) :: this
339  integer(I4B), intent(in) :: ic
340  type(CellDefnType), pointer, intent(inout) :: defn
341 
342  call this%load_cell_properties(ic, defn)
343  call this%fmi%dis%get_polyverts( &
344  defn%icell, &
345  defn%polyvert, &
346  closed=.true.)
347  call this%load_cell_neighbors(defn)
348  call this%load_cell_saturation_status(defn)
349  defn%ispv180(1:defn%npolyverts + 1) = .false.
350  call this%load_cell_flows(defn)
351 

◆ load_cell_face_flows()

subroutine methoddismodule::load_cell_face_flows ( class(methoddistype), intent(inout)  this,
type(celldefntype), intent(inout), pointer  defn 
)
private

Definition at line 483 of file MethodDis.f90.

484  ! dummy
485  class(MethodDisType), intent(inout) :: this
486  type(CellDefnType), pointer, intent(inout) :: defn
487  ! local
488  integer(I4B) :: m, n, nfaces
489  real(DP) :: q
490 
491  nfaces = defn%npolyverts + 3
492  do m = 1, nfaces
493  n = defn%facenbr(m)
494  if (n > 0) then
495  q = this%fmi%gwfflowja(this%fmi%dis%con%ia(defn%icell) + n)
496  defn%faceflow(m) = defn%faceflow(m) + q
497  end if
498  end do

◆ load_cell_flows()

subroutine methoddismodule::load_cell_flows ( class(methoddistype), intent(inout)  this,
type(celldefntype), intent(inout), pointer  defn 
)
private

Definition at line 456 of file MethodDis.f90.

457  class(MethodDisType), intent(inout) :: this
458  type(CellDefnType), pointer, intent(inout) :: defn
459 
460  ! Load face flows, including boundary flows. As with cell verts,
461  ! the face flow array wraps around. Top and bottom flows make up
462  ! the last two elements, respectively, for size npolyverts + 3.
463  ! If there is no flow through any face, set a no-exit-face flag.
464  defn%faceflow = dzero
465  call this%load_cell_boundary_flows(defn)
466  call this%load_cell_face_flows(defn)
467  call this%cap_cell_wt_flow(defn)
468  call this%load_cell_no_exit_face(defn)
469 
470  ! Add up net distributed flow
471  defn%distflow = this%fmi%SourceFlows(defn%icell) + &
472  this%fmi%SinkFlows(defn%icell) + &
473  this%fmi%StorageFlows(defn%icell)
474 
475  ! Set weak sink flag
476  if (this%fmi%SinkFlows(defn%icell) .ne. dzero) then
477  defn%iweaksink = 1
478  else
479  defn%iweaksink = 0
480  end if

◆ load_cell_neighbors()

subroutine methoddismodule::load_cell_neighbors ( class(methoddistype), intent(inout)  this,
type(celldefntype), intent(inout), pointer  defn 
)
private

Definition at line 387 of file MethodDis.f90.

388  ! dummy
389  class(MethodDisType), intent(inout) :: this
390  type(CellDefnType), pointer, intent(inout) :: defn
391  ! local
392  integer(I4B) :: ic1
393  integer(I4B) :: ic2
394  integer(I4B) :: icu1
395  integer(I4B) :: icu2
396  integer(I4B) :: j1
397  integer(I4B) :: iloc
398  integer(I4B) :: ipos
399  integer(I4B) :: irow1
400  integer(I4B) :: irow2
401  integer(I4B) :: jcol1
402  integer(I4B) :: jcol2
403  integer(I4B) :: klay1
404  integer(I4B) :: klay2
405  integer(I4B) :: iedgeface
406 
407  select type (dis => this%fmi%dis)
408  type is (distype)
409  ! Load face neighbors
410  defn%facenbr = 0
411  ic1 = defn%icell
412  icu1 = dis%get_nodeuser(ic1)
413  call get_ijk(icu1, dis%nrow, dis%ncol, dis%nlay, &
414  irow1, jcol1, klay1)
415  call get_jk(icu1, dis%get_ncpl(), dis%nlay, j1, klay1)
416  do iloc = 1, dis%con%ia(ic1 + 1) - dis%con%ia(ic1) - 1
417  ipos = dis%con%ia(ic1) + iloc
418  ! mask could become relevant if PRT uses interface model
419  if (dis%con%mask(ipos) == 0) cycle
420  ic2 = dis%con%ja(ipos)
421  icu2 = dis%get_nodeuser(ic2)
422  call get_ijk(icu2, dis%nrow, dis%ncol, dis%nlay, &
423  irow2, jcol2, klay2)
424  if (klay2 == klay1) then
425  ! Edge (polygon) face neighbor
426  if (irow2 > irow1) then
427  ! Neighbor to the S
428  iedgeface = 4
429  else if (jcol2 > jcol1) then
430  ! Neighbor to the E
431  iedgeface = 3
432  else if (irow2 < irow1) then
433  ! Neighbor to the N
434  iedgeface = 2
435  else
436  ! Neighbor to the W
437  iedgeface = 1
438  end if
439  defn%facenbr(iedgeface) = int(iloc, 1)
440  else if (klay2 > klay1) then
441  ! Bottom face neighbor
442  defn%facenbr(defn%npolyverts + 2) = int(iloc, 1)
443  else
444  ! Top face neighbor
445  defn%facenbr(defn%npolyverts + 3) = int(iloc, 1)
446  end if
447  end do
448  end select
449  ! List of edge (polygon) faces wraps around
450  defn%facenbr(defn%npolyverts + 1) = defn%facenbr(1)
Here is the call graph for this function:

◆ load_cell_properties()

subroutine methoddismodule::load_cell_properties ( class(methoddistype), intent(inout)  this,
integer(i4b), intent(in)  ic,
type(celldefntype), intent(inout), pointer  defn 
)
private

Definition at line 354 of file MethodDis.f90.

355  ! dummy
356  class(MethodDisType), intent(inout) :: this
357  integer(I4B), intent(in) :: ic
358  type(CellDefnType), pointer, intent(inout) :: defn
359  ! local
360  integer(I4B) :: irow, icol, ilay, icu
361 
362  defn%icell = ic
363  defn%npolyverts = 4 ! rectangular cell always has 4 vertices
364  defn%iatop = get_iatop(this%fmi%dis%get_ncpl(), &
365  this%fmi%dis%get_nodeuser(ic))
366  defn%top = this%fmi%dis%bot(ic) + &
367  this%fmi%gwfsat(ic) * &
368  (this%fmi%dis%top(ic) - this%fmi%dis%bot(ic))
369  defn%bot = this%fmi%dis%bot(ic)
370  defn%sat = this%fmi%gwfsat(ic)
371  defn%porosity = this%porosity(ic)
372  defn%retfactor = this%retfactor(ic)
373  select type (dis => this%fmi%dis)
374  type is (distype)
375  icu = dis%get_nodeuser(ic)
376  call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, irow, icol, ilay)
377  defn%ilay = ilay
378  end select
379  defn%izone = this%izone(ic)
380  defn%can_be_rect = .true.
381  defn%can_be_quad = .false.
382 
Here is the call graph for this function:

◆ load_dis()

subroutine methoddismodule::load_dis ( class(methoddistype), intent(inout)  this,
type(particletype), intent(inout), pointer  particle,
integer(i4b), intent(in)  next_level,
class(methodtype), intent(inout), pointer  submethod 
)
private

Definition at line 152 of file MethodDis.f90.

153  ! dummy
154  class(MethodDisType), intent(inout) :: this
155  type(ParticleType), pointer, intent(inout) :: particle
156  integer(I4B), intent(in) :: next_level
157  class(MethodType), pointer, intent(inout) :: submethod
158  ! local
159  integer(I4B) :: ic
160 
161  select type (cell => this%cell)
162  type is (cellrecttype)
163  ic = particle%itrdomain(next_level)
164  call this%load_cell_defn(ic, cell%defn)
165  call this%load_cell(ic, cell)
166  if (this%fmi%ibdgwfsat0(ic) == 0) then
167  call this%method_cell_ptb%init( &
168  fmi=this%fmi, &
169  cell=this%cell, &
170  events=this%events, &
171  tracktimes=this%tracktimes)
172  submethod => this%method_cell_ptb
173  else
174  call this%method_cell_plck%init( &
175  fmi=this%fmi, &
176  cell=this%cell, &
177  events=this%events, &
178  tracktimes=this%tracktimes)
179  submethod => this%method_cell_plck
180  end if
181  end select

◆ load_particle()

subroutine methoddismodule::load_particle ( class(methoddistype), intent(inout)  this,
type(cellrecttype), intent(inout), pointer  cell,
type(particletype), intent(inout), pointer  particle 
)
private

Definition at line 186 of file MethodDis.f90.

187  use particlemodule, only: term_boundary
188  use particleeventmodule, only: terminate
189  ! dummy
190  class(MethodDisType), intent(inout) :: this
191  type(CellRectType), pointer, intent(inout) :: cell
192  type(ParticleType), pointer, intent(inout) :: particle
193  ! local
194  integer(I4B) :: ic
195  integer(I4B) :: icu
196  integer(I4B) :: ilay
197  integer(I4B) :: irow
198  integer(I4B) :: icol
199  integer(I4B) :: inface
200  integer(I4B) :: idiag
201  integer(I4B) :: inbr
202  integer(I4B) :: ipos
203  real(DP) :: z
204  real(DP) :: zrel
205  real(DP) :: topfrom
206  real(DP) :: botfrom
207  real(DP) :: top
208  real(DP) :: bot
209  real(DP) :: sat
210 
211  select type (dis => this%fmi%dis)
212  type is (distype)
213  ! compute reduced/user node numbers and layer
214  inface = particle%iboundary(level_feature)
215  inbr = cell%defn%facenbr(inface)
216  idiag = this%fmi%dis%con%ia(cell%defn%icell)
217  ipos = idiag + inbr
218  ic = dis%con%ja(ipos)
219  icu = dis%get_nodeuser(ic)
220  call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, irow, icol, ilay)
221 
222  ! update node numbers and layer
223  particle%itrdomain(level_feature) = ic
224  particle%icu = icu
225  particle%ilay = ilay
226 
227  ! map/set particle entry face
228  if (inface .eq. 1) then
229  inface = 3
230  else if (inface .eq. 2) then
231  inface = 4
232  else if (inface .eq. 3) then
233  inface = 1
234  else if (inface .eq. 4) then
235  inface = 2
236  else if (inface .eq. 6) then
237  inface = 7
238  else if (inface .eq. 7) then
239  inface = 6
240  end if
241  particle%iboundary(level_feature) = inface
242 
243  ! map z between cells
244  z = particle%z
245  if (inface < 5) then
246  topfrom = cell%defn%top
247  botfrom = cell%defn%bot
248  zrel = (z - botfrom) / (topfrom - botfrom)
249  top = dis%top(ic)
250  bot = dis%bot(ic)
251  sat = this%fmi%gwfsat(ic)
252  z = bot + zrel * sat * (top - bot)
253  end if
254  particle%z = z
255  end select
@, public terminate
particle terminated
@ term_boundary
terminated at a boundary face
Definition: Particle.f90:31
Here is the call graph for this function:

◆ pass_dis()

subroutine methoddismodule::pass_dis ( class(methoddistype), intent(inout)  this,
type(particletype), intent(inout), pointer  particle 
)
private

Definition at line 285 of file MethodDis.f90.

286  use particlemodule, only: term_boundary
287  use particleeventmodule, only: terminate
288  ! dummy
289  class(MethodDisType), intent(inout) :: this
290  type(ParticleType), pointer, intent(inout) :: particle
291  ! local
292  type(CellRectType), pointer :: cell
293  integer(I4B) :: iface
294  logical(LGP) :: no_neighbors
295 
296  select type (c => this%cell)
297  type is (cellrecttype)
298  cell => c
299  iface = particle%iboundary(level_feature)
300  no_neighbors = cell%defn%facenbr(iface) == 0
301 
302  ! todo AMP: reconsider when multiple models supported
303  if (no_neighbors) then
304  call this%terminate(particle, status=term_boundary)
305  return
306  end if
307 
308  call this%load_particle(cell, particle)
309  call this%update_flowja(cell, particle)
310  end select

◆ update_flowja()

subroutine methoddismodule::update_flowja ( class(methoddistype), intent(inout)  this,
type(cellrecttype), intent(inout), pointer  cell,
type(particletype), intent(inout), pointer  particle 
)

Definition at line 260 of file MethodDis.f90.

261  ! dummy
262  class(MethodDisType), intent(inout) :: this
263  type(CellRectType), pointer, intent(inout) :: cell
264  type(ParticleType), pointer, intent(inout) :: particle
265  ! local
266  integer(I4B) :: idiag
267  integer(I4B) :: inbr
268  integer(I4B) :: inface
269  integer(I4B) :: ipos
270 
271  inface = particle%iboundary(level_feature)
272  inbr = cell%defn%facenbr(inface)
273  idiag = this%fmi%dis%con%ia(cell%defn%icell)
274  ipos = idiag + inbr
275 
276  ! leaving old cell
277  this%flowja(ipos) = this%flowja(ipos) - done
278 
279  ! entering new cell
280  this%flowja(this%fmi%dis%con%isym(ipos)) &
281  = this%flowja(this%fmi%dis%con%isym(ipos)) + done