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