MODFLOW 6  version 6.7.0.dev1
USGS Modular Hydrologic Model
MethodDis.f90
Go to the documentation of this file.
2 
3  use kindmodule, only: dp, i4b, lgp
4  use constantsmodule, only: done, dzero
5  use methodmodule, only: methodtype
7  use cellmodule
10  use particlemodule
11  use prtfmimodule, only: prtfmitype
12  use dismodule, only: distype
13  use geomutilmodule, only: get_ijk, get_jk
14  use mathutilmodule, only: is_close
15  implicit none
16 
17  private
18  public :: methoddistype
19  public :: create_method_dis
20 
21  type, extends(methodtype) :: methoddistype
22  private
23  contains
24  procedure, public :: apply => apply_dis !< apply the DIS tracking method
25  procedure, public :: deallocate !< deallocate arrays and scalars
26  procedure, public :: load => load_dis !< load the method
27  procedure, public :: pass => pass_dis !< pass the particle to the next domain
28  procedure :: get_top !< get cell top elevation
29  procedure :: update_flowja !< load intercell mass flows
30  procedure :: load_particle !< load particle properties
31  procedure :: load_properties !< load cell properties
32  procedure :: load_neighbors !< load cell face neighbors
33  procedure :: load_flows !< load cell face flows
34  procedure :: load_boundary_flows_to_defn !< load boundary flows to the cell definition
35  procedure :: load_face_flows_to_defn !< load face flows to the cell definition
36  procedure :: load_celldefn !< load cell definition from the grid
37  procedure :: load_cell !< load cell geometry and flows
38  end type methoddistype
39 
40 contains
41 
42  !> @brief Create a new structured grid (DIS) tracking method
43  subroutine create_method_dis(method)
44  ! dummy
45  type(methoddistype), pointer :: method
46  ! local
47  type(cellrecttype), pointer :: cell
48 
49  allocate (method)
50  allocate (method%name)
51  call create_cell_rect(cell)
52  method%cell => cell
53  method%name = "dis"
54  method%delegates = .true.
55  end subroutine create_method_dis
56 
57  !> @brief Destructor the tracking method
58  subroutine deallocate (this)
59  class(methoddistype), intent(inout) :: this
60  deallocate (this%name)
61  end subroutine deallocate
62 
63  subroutine load_cell(this, ic, cell)
64  ! dummy
65  class(methoddistype), intent(inout) :: this
66  integer(I4B), intent(in) :: ic
67  type(cellrecttype), intent(inout) :: cell
68  ! local
69  integer(I4B) :: icu
70  integer(I4B) :: irow
71  integer(I4B) :: jcol
72  integer(I4B) :: klay
73  real(DP) :: areax
74  real(DP) :: areay
75  real(DP) :: areaz
76  real(DP) :: dx
77  real(DP) :: dy
78  real(DP) :: dz
79  real(DP) :: factor
80  real(DP) :: term
81 
82  select type (dis => this%fmi%dis)
83  type is (distype)
84  icu = dis%get_nodeuser(ic)
85 
86  call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, &
87  irow, jcol, klay)
88 
89  dx = dis%delr(jcol)
90  dy = dis%delc(irow)
91  dz = cell%defn%top - cell%defn%bot
92 
93  cell%dx = dx
94  cell%dy = dy
95  cell%dz = dz
96  cell%sinrot = dzero
97  cell%cosrot = done
98  cell%xOrigin = cell%defn%polyvert(1, 1)
99  cell%yOrigin = cell%defn%polyvert(2, 1)
100  cell%zOrigin = cell%defn%bot
101  cell%ipvOrigin = 1
102 
103  factor = done / cell%defn%retfactor
104  factor = factor / cell%defn%porosity
105 
106  areaz = dx * dy
107  term = factor / areaz
108 
109  cell%vz1 = cell%defn%faceflow(6) * term
110  cell%vz2 = -cell%defn%faceflow(7) * term
111 
112  if (this%fmi%ibdgwfsat0(ic) == 0) then
113  cell%vx1 = dzero
114  cell%vx2 = dzero
115  cell%vy1 = dzero
116  cell%vy2 = dzero
117  cell%vz1 = dzero
118  cell%vz2 = dzero
119  return
120  end if
121 
122  areax = dy * dz
123  term = factor / areax
124  cell%vx1 = cell%defn%faceflow(1) * term
125  cell%vx2 = -cell%defn%faceflow(3) * term
126 
127  areay = dx * dz
128  term = factor / areay
129  cell%vy1 = cell%defn%faceflow(4) * term
130  cell%vy2 = -cell%defn%faceflow(2) * term
131 
132  end select
133  end subroutine load_cell
134 
135  !> @brief Load the cell geometry and tracking method
136  subroutine load_dis(this, particle, next_level, submethod)
137  ! dummy
138  class(methoddistype), intent(inout) :: this
139  type(particletype), pointer, intent(inout) :: particle
140  integer(I4B), intent(in) :: next_level
141  class(methodtype), pointer, intent(inout) :: submethod
142  ! local
143  integer(I4B) :: ic
144 
145  select type (cell => this%cell)
146  type is (cellrecttype)
147  ic = particle%idomain(next_level)
148  call this%load_celldefn(ic, cell%defn)
149  call this%load_cell(ic, cell)
150  if (this%fmi%ibdgwfsat0(ic) == 0) then
151  call method_cell_ptb%init( &
152  fmi=this%fmi, &
153  cell=this%cell, &
154  trackctl=this%trackctl, &
155  tracktimes=this%tracktimes)
156  submethod => method_cell_ptb
157  else
158  call method_cell_plck%init( &
159  fmi=this%fmi, &
160  cell=this%cell, &
161  trackctl=this%trackctl, &
162  tracktimes=this%tracktimes)
163  submethod => method_cell_plck
164  end if
165  end select
166  end subroutine load_dis
167 
168  !> @brief Load cell properties into the particle, including
169  ! the z coordinate, entry face, and node and layer numbers.
170  subroutine load_particle(this, cell, particle)
171  use particlemodule, only: term_boundary
172  ! dummy
173  class(methoddistype), intent(inout) :: this
174  type(cellrecttype), pointer, intent(inout) :: cell
175  type(particletype), pointer, intent(inout) :: particle
176  ! local
177  integer(I4B) :: ic
178  integer(I4B) :: icu
179  integer(I4B) :: ilay
180  integer(I4B) :: irow
181  integer(I4B) :: icol
182  integer(I4B) :: inface
183  integer(I4B) :: idiag
184  integer(I4B) :: inbr
185  integer(I4B) :: ipos
186  real(DP) :: z
187  real(DP) :: zrel
188  real(DP) :: topfrom
189  real(DP) :: botfrom
190  real(DP) :: top
191  real(DP) :: bot
192  real(DP) :: sat
193 
194  select type (dis => this%fmi%dis)
195  type is (distype)
196  ! compute reduced/user node numbers and layer
197  inface = particle%iboundary(2)
198  inbr = cell%defn%facenbr(inface)
199  idiag = this%fmi%dis%con%ia(cell%defn%icell)
200  ipos = idiag + inbr
201  ic = dis%con%ja(ipos)
202  icu = dis%get_nodeuser(ic)
203  call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, &
204  irow, icol, ilay)
205 
206  ! if returning to a cell through the bottom
207  ! face after previously leaving it through
208  ! that same face, we've entered a cycle
209  ! as can occur e.g. in wells. terminate
210  ! in the previous cell.
211  if (ic == particle%icp .and. inface == 7 .and. ilay < particle%ilay) then
212  particle%advancing = .false.
213  particle%idomain(2) = particle%icp
214  particle%istatus = term_boundary
215  particle%izone = particle%izp
216  call this%save(particle, reason=3)
217  return
218  else
219  particle%icp = particle%idomain(2)
220  particle%izp = particle%izone
221  end if
222 
223  ! update node numbers and layer
224  particle%idomain(2) = ic
225  particle%icu = icu
226  particle%ilay = ilay
227 
228  ! map/set particle entry face
229  if (inface .eq. 1) then
230  inface = 3
231  else if (inface .eq. 2) then
232  inface = 4
233  else if (inface .eq. 3) then
234  inface = 1
235  else if (inface .eq. 4) then
236  inface = 2
237  else if (inface .eq. 6) then
238  inface = 7
239  else if (inface .eq. 7) then
240  inface = 6
241  end if
242  particle%iboundary(2) = inface
243 
244  ! map z between cells
245  z = particle%z
246  if (inface < 5) then
247  topfrom = cell%defn%top
248  botfrom = cell%defn%bot
249  zrel = (z - botfrom) / (topfrom - botfrom)
250  top = dis%top(ic)
251  bot = dis%bot(ic)
252  sat = this%fmi%gwfsat(ic)
253  z = bot + zrel * sat * (top - bot)
254  end if
255  particle%z = z
256  end select
257  end subroutine load_particle
258 
259  !> @brief Update cell-cell flows of particle mass.
260  !! Every particle is currently assigned unit mass.
261  subroutine update_flowja(this, cell, particle)
262  ! dummy
263  class(methoddistype), intent(inout) :: this
264  type(cellrecttype), pointer, intent(inout) :: cell
265  type(particletype), pointer, intent(inout) :: particle
266  ! local
267  integer(I4B) :: idiag
268  integer(I4B) :: inbr
269  integer(I4B) :: inface
270  integer(I4B) :: ipos
271 
272  inface = particle%iboundary(2)
273  inbr = cell%defn%facenbr(inface)
274  idiag = this%fmi%dis%con%ia(cell%defn%icell)
275  ipos = idiag + inbr
276 
277  ! leaving old cell
278  this%flowja(ipos) = this%flowja(ipos) - done
279 
280  ! entering new cell
281  this%flowja(this%fmi%dis%con%isym(ipos)) &
282  = this%flowja(this%fmi%dis%con%isym(ipos)) + done
283  end subroutine update_flowja
284 
285  !> @brief Pass a particle to the next cell, if there is one
286  subroutine pass_dis(this, particle)
287  use particlemodule, only: term_boundary
288  ! dummy
289  class(methoddistype), intent(inout) :: this
290  type(particletype), pointer, intent(inout) :: particle
291  ! local
292  type(cellrecttype), pointer :: cell
293 
294  select type (c => this%cell)
295  type is (cellrecttype)
296  cell => c
297  ! If the entry face has no neighbors it's a
298  ! boundary face, so terminate the particle.
299  ! todo AMP: reconsider when multiple models supported
300  if (cell%defn%facenbr(particle%iboundary(2)) .eq. 0) then
301  particle%istatus = term_boundary
302  particle%advancing = .false.
303  call this%save(particle, reason=3)
304  else
305  ! Update old to new cell properties
306  call this%load_particle(cell, particle)
307  if (.not. particle%advancing) return
308 
309  ! Update intercell mass flows
310  call this%update_flowja(cell, particle)
311  end if
312  end select
313  end subroutine pass_dis
314 
315  !> @brief Apply the structured tracking method to a particle.
316  subroutine apply_dis(this, particle, tmax)
317  class(methoddistype), intent(inout) :: this
318  type(particletype), pointer, intent(inout) :: particle
319  real(DP), intent(in) :: tmax
320 
321  call this%track(particle, 1, tmax)
322  end subroutine apply_dis
323 
324  !> @brief Returns a top elevation based on index iatop
325  function get_top(this, iatop) result(top)
326  class(methoddistype), intent(inout) :: this
327  integer, intent(in) :: iatop
328  double precision :: top
329 
330  if (iatop .lt. 0) then
331  top = this%fmi%dis%top(-iatop)
332  else
333  top = this%fmi%dis%bot(iatop)
334  end if
335  end function get_top
336 
337  !> @brief Loads cell definition from the grid
338  subroutine load_celldefn(this, ic, defn)
339  ! dummy
340  class(methoddistype), intent(inout) :: this
341  integer(I4B), intent(in) :: ic
342  type(celldefntype), pointer, intent(inout) :: defn
343 
344  ! Load basic cell properties
345  call this%load_properties(ic, defn)
346 
347  ! Load cell polygon vertices
348  call this%fmi%dis%get_polyverts( &
349  defn%icell, &
350  defn%polyvert, &
351  closed=.true.)
352  call this%load_neighbors(defn)
353 
354  ! Load 180 degree face indicators
355  defn%ispv180(1:defn%npolyverts + 1) = .false.
356 
357  call this%load_flows(defn)
358 
359  end subroutine load_celldefn
360 
361  subroutine load_properties(this, ic, defn)
362  ! dummy
363  class(methoddistype), intent(inout) :: this
364  integer(I4B), intent(in) :: ic
365  type(celldefntype), pointer, intent(inout) :: defn
366  ! local
367  integer(I4B) :: irow, icol, ilay, icu
368 
369  defn%icell = ic
370  defn%npolyverts = 4 ! rectangular cell always has 4 vertices
371  defn%iatop = get_iatop(this%fmi%dis%get_ncpl(), &
372  this%fmi%dis%get_nodeuser(ic))
373  defn%top = this%fmi%dis%bot(ic) + &
374  this%fmi%gwfsat(ic) * &
375  (this%fmi%dis%top(ic) - this%fmi%dis%bot(ic))
376  defn%bot = this%fmi%dis%bot(ic)
377  defn%sat = this%fmi%gwfsat(ic)
378  defn%porosity = this%porosity(ic)
379  defn%retfactor = this%retfactor(ic)
380  select type (dis => this%fmi%dis)
381  type is (distype)
382  icu = dis%get_nodeuser(ic)
383  call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, irow, icol, ilay)
384  defn%ilay = ilay
385  end select
386  defn%izone = this%izone(ic)
387  defn%can_be_rect = .true.
388  defn%can_be_quad = .false.
389  end subroutine load_properties
390 
391  !> @brief Loads face neighbors to cell definition from the grid.
392  !! Assumes cell index and number of vertices are already loaded.
393  subroutine load_neighbors(this, defn)
394  ! dummy
395  class(methoddistype), intent(inout) :: this
396  type(celldefntype), pointer, intent(inout) :: defn
397  ! local
398  integer(I4B) :: ic1
399  integer(I4B) :: ic2
400  integer(I4B) :: icu1
401  integer(I4B) :: icu2
402  integer(I4B) :: j1
403  integer(I4B) :: iloc
404  integer(I4B) :: ipos
405  integer(I4B) :: irow1
406  integer(I4B) :: irow2
407  integer(I4B) :: jcol1
408  integer(I4B) :: jcol2
409  integer(I4B) :: klay1
410  integer(I4B) :: klay2
411  integer(I4B) :: iedgeface
412 
413  select type (dis => this%fmi%dis)
414  type is (distype)
415  ! Load face neighbors
416  defn%facenbr = 0
417  ic1 = defn%icell
418  icu1 = dis%get_nodeuser(ic1)
419  call get_ijk(icu1, dis%nrow, dis%ncol, dis%nlay, &
420  irow1, jcol1, klay1)
421  call get_jk(icu1, dis%get_ncpl(), dis%nlay, j1, klay1)
422  do iloc = 1, dis%con%ia(ic1 + 1) - dis%con%ia(ic1) - 1
423  ipos = dis%con%ia(ic1) + iloc
424  ! mask could become relevant if PRT uses interface model
425  if (dis%con%mask(ipos) == 0) cycle
426  ic2 = dis%con%ja(ipos)
427  icu2 = dis%get_nodeuser(ic2)
428  call get_ijk(icu2, dis%nrow, dis%ncol, dis%nlay, &
429  irow2, jcol2, klay2)
430  if (klay2 == klay1) then
431  ! Edge (polygon) face neighbor
432  if (irow2 > irow1) then
433  ! Neighbor to the S
434  iedgeface = 4
435  else if (jcol2 > jcol1) then
436  ! Neighbor to the E
437  iedgeface = 3
438  else if (irow2 < irow1) then
439  ! Neighbor to the N
440  iedgeface = 2
441  else
442  ! Neighbor to the W
443  iedgeface = 1
444  end if
445  defn%facenbr(iedgeface) = int(iloc, 1)
446  else if (klay2 > klay1) then
447  ! Bottom face neighbor
448  defn%facenbr(defn%npolyverts + 2) = int(iloc, 1)
449  else
450  ! Top face neighbor
451  defn%facenbr(defn%npolyverts + 3) = int(iloc, 1)
452  end if
453  end do
454  end select
455  ! List of edge (polygon) faces wraps around
456  defn%facenbr(defn%npolyverts + 1) = defn%facenbr(1)
457  end subroutine load_neighbors
458 
459  !> @brief Load flows into the cell definition.
460  !! These include face, boundary and net distributed flows.
461  !! Assumes cell index and number of vertices are already loaded.
462  subroutine load_flows(this, defn)
463  class(methoddistype), intent(inout) :: this
464  type(celldefntype), intent(inout) :: defn
465 
466  ! Load face flows, including boundary flows. As with cell verts,
467  ! the face flow array wraps around. Top and bottom flows make up
468  ! the last two elements, respectively, for size npolyverts + 3.
469  ! If there is no flow through any face, set a no-exit-face flag.
470  defn%faceflow = dzero
471  defn%inoexitface = 1
472  call this%load_boundary_flows_to_defn(defn)
473  call this%load_face_flows_to_defn(defn)
474 
475  ! Add up net distributed flow
476  defn%distflow = this%fmi%SourceFlows(defn%icell) + &
477  this%fmi%SinkFlows(defn%icell) + &
478  this%fmi%StorageFlows(defn%icell)
479 
480  ! Set weak sink flag
481  if (this%fmi%SinkFlows(defn%icell) .ne. dzero) then
482  defn%iweaksink = 1
483  else
484  defn%iweaksink = 0
485  end if
486  end subroutine load_flows
487 
488  subroutine load_face_flows_to_defn(this, defn)
489  ! dummy
490  class(methoddistype), intent(inout) :: this
491  type(celldefntype), intent(inout) :: defn
492  ! local
493  integer(I4B) :: m, n, nfaces
494  real(DP) :: q
495 
496  nfaces = defn%npolyverts + 3
497  do m = 1, nfaces
498  n = defn%facenbr(m)
499  if (n > 0) then
500  q = this%fmi%gwfflowja(this%fmi%dis%con%ia(defn%icell) + n)
501  defn%faceflow(m) = defn%faceflow(m) + q
502  end if
503  if (defn%faceflow(m) < dzero) defn%inoexitface = 0
504  end do
505  end subroutine load_face_flows_to_defn
506 
507  !> @brief Add boundary flows to the cell definition faceflow array.
508  !! Assumes cell index and number of vertices are already loaded.
509  subroutine load_boundary_flows_to_defn(this, defn)
510  ! dummy
511  class(methoddistype), intent(inout) :: this
512  type(celldefntype), intent(inout) :: defn
513  ! local
514  integer(I4B) :: ioffset
515 
516  ioffset = (defn%icell - 1) * max_poly_cells
517  defn%faceflow(1) = defn%faceflow(1) + &
518  this%fmi%BoundaryFlows(ioffset + 1)
519  defn%faceflow(2) = defn%faceflow(2) + &
520  this%fmi%BoundaryFlows(ioffset + 2)
521  defn%faceflow(3) = defn%faceflow(3) + &
522  this%fmi%BoundaryFlows(ioffset + 3)
523  defn%faceflow(4) = defn%faceflow(4) + &
524  this%fmi%BoundaryFlows(ioffset + 4)
525  defn%faceflow(5) = defn%faceflow(1)
526  defn%faceflow(6) = defn%faceflow(6) + &
527  this%fmi%BoundaryFlows(ioffset + 9)
528  defn%faceflow(7) = defn%faceflow(7) + &
529  this%fmi%BoundaryFlows(ioffset + 10)
530  end subroutine load_boundary_flows_to_defn
531 
532 end module methoddismodule
integer(i4b) function, public get_iatop(ncpl, icu)
Get the index corresponding to top elevation of a cell in the grid. This is -1 if the cell is in the ...
Definition: CellDefn.f90:54
integer(i4b), parameter, public max_poly_cells
Definition: Cell.f90:9
subroutine, public create_cell_rect(cell)
Create a new rectangular cell.
Definition: CellRect.f90:40
This module contains simulation constants.
Definition: Constants.f90:9
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
real(dp), parameter done
real constant 1
Definition: Constants.f90:76
Definition: Dis.f90:1
subroutine, public get_ijk(nodenumber, nrow, ncol, nlay, irow, icol, ilay)
Get row, column and layer indices from node number and grid dimensions. If nodenumber is invalid,...
Definition: GeomUtil.f90:100
subroutine, public get_jk(nodenumber, ncpl, nlay, icpl, ilay)
Get layer index and within-layer node index from node number and grid dimensions. If nodenumber is in...
Definition: GeomUtil.f90:128
This module defines variable data types.
Definition: kind.f90:8
pure logical function, public is_close(a, b, rtol, atol, symmetric)
Check if a real value is approximately equal to another.
Definition: MathUtil.f90:46
Cell-level tracking methods.
type(methodcellpollocktype), pointer, public method_cell_plck
type(methodcellpasstobottype), pointer, public method_cell_ptb
subroutine load_flows(this, defn)
Load flows into the cell definition. These include face, boundary and net distributed flows....
Definition: MethodDis.f90:463
subroutine load_boundary_flows_to_defn(this, defn)
Add boundary flows to the cell definition faceflow array. Assumes cell index and number of vertices a...
Definition: MethodDis.f90:510
double precision function get_top(this, iatop)
Returns a top elevation based on index iatop.
Definition: MethodDis.f90:326
subroutine load_particle(this, cell, particle)
Load cell properties into the particle, including.
Definition: MethodDis.f90:171
subroutine load_properties(this, ic, defn)
Definition: MethodDis.f90:362
subroutine update_flowja(this, cell, particle)
Update cell-cell flows of particle mass. Every particle is currently assigned unit mass.
Definition: MethodDis.f90:262
subroutine, public create_method_dis(method)
Create a new structured grid (DIS) tracking method.
Definition: MethodDis.f90:44
subroutine load_celldefn(this, ic, defn)
Loads cell definition from the grid.
Definition: MethodDis.f90:339
subroutine apply_dis(this, particle, tmax)
Apply the structured tracking method to a particle.
Definition: MethodDis.f90:317
subroutine load_dis(this, particle, next_level, submethod)
Load the cell geometry and tracking method.
Definition: MethodDis.f90:137
subroutine load_cell(this, ic, cell)
Definition: MethodDis.f90:64
subroutine load_face_flows_to_defn(this, defn)
Definition: MethodDis.f90:489
subroutine load_neighbors(this, defn)
Loads face neighbors to cell definition from the grid. Assumes cell index and number of vertices are ...
Definition: MethodDis.f90:394
subroutine pass_dis(this, particle)
Pass a particle to the next cell, if there is one.
Definition: MethodDis.f90:287
Particle tracking strategies.
Definition: Method.f90:2
@ term_boundary
terminated at a boundary face
Definition: Particle.f90:33
Base grid cell definition.
Definition: CellDefn.f90:10
Structured grid discretization.
Definition: Dis.f90:23
Base type for particle tracking methods.
Definition: Method.f90:31
Particle tracked by the PRT model.
Definition: Particle.f90:78