MODFLOW 6  version 6.7.0.dev1
USGS Modular Hydrologic Model
MethodDisv.f90
Go to the documentation of this file.
2 
3  use kindmodule, only: dp, i4b, lgp
4  use errorutilmodule, only: pstop
5  use constantsmodule, only: done, dzero
6  use methodmodule, only: methodtype
8  use cellmodule, only: max_poly_cells
10  use cellpolymodule
11  use particlemodule
12  use prtfmimodule, only: prtfmitype
13  use disvmodule, only: disvtype
16  implicit none
17 
18  private
19  public :: methoddisvtype
20  public :: create_method_disv
21 
22  type, extends(methodtype) :: methoddisvtype
23  private
24  type(celldefntype), pointer :: neighbor => null() !< ptr to a neighbor defn
25  contains
26  procedure, public :: apply => apply_disv !< apply the DISV tracking method
27  procedure, public :: deallocate !< deallocate arrays and scalars
28  procedure, public :: load => load_disv !< load the cell method
29  procedure, public :: load_cell_defn !< load cell definition from the grid
30  procedure, public :: pass => pass_disv !< pass the particle to the next cell
31  procedure :: map_neighbor !< map a location on the cell face to the shared face of a neighbor
32  procedure :: update_flowja !< update intercell mass flows
33  procedure :: load_particle !< load particle properties
34  procedure :: load_properties !< load cell properties
35  procedure :: load_polygon !< load cell polygon
36  procedure :: load_neighbors !< load cell face neighbors
37  procedure :: load_indicators !< load cell 180-degree vertex indicator
38  procedure :: load_flows !< load the cell's flows
39  procedure :: load_boundary_flows_to_defn_rect !< load boundary flows to a rectangular cell definition
40  procedure :: load_boundary_flows_to_defn_rect_quad !< load boundary flows to a rectangular-quad cell definition
41  procedure :: load_boundary_flows_to_defn_poly !< load boundary flows to a polygonal cell definition
42  procedure :: load_face_flows_to_defn_poly !< load face flows to a polygonal cell definition
43  end type methoddisvtype
44 
45 contains
46 
47  !> @brief Create a new vertex grid (DISV) tracking method
48  subroutine create_method_disv(method)
49  ! dummy
50  type(methoddisvtype), pointer :: method
51  ! local
52  type(cellpolytype), pointer :: cell
53 
54  allocate (method)
55  allocate (method%name)
56  call create_cell_poly(cell)
57  method%cell => cell
58  method%name = "disv"
59  method%delegates = .true.
60  call create_defn(method%neighbor)
61  end subroutine create_method_disv
62 
63  !> @brief Destroy the tracking method
64  subroutine deallocate (this)
65  class(methoddisvtype), intent(inout) :: this
66  deallocate (this%name)
67  end subroutine deallocate
68 
69  !> @brief Load the cell and the tracking method
70  subroutine load_disv(this, particle, next_level, submethod)
71  use cellmodule
72  use cellrectmodule
74  use cellutilmodule
75  ! dummy
76  class(methoddisvtype), intent(inout) :: this
77  type(particletype), pointer, intent(inout) :: particle
78  integer(I4B), intent(in) :: next_level
79  class(methodtype), pointer, intent(inout) :: submethod
80  ! local
81  integer(I4B) :: ic
82  class(celltype), pointer :: base
83  type(cellrecttype), pointer :: rect
84  type(cellrectquadtype), pointer :: quad
85 
86  select type (cell => this%cell)
87  type is (cellpolytype)
88  ic = particle%idomain(next_level)
89  call this%load_cell_defn(ic, cell%defn)
90  if (this%fmi%ibdgwfsat0(ic) == 0) then
91  ! Cell is active but dry, so select and initialize pass-to-bottom
92  ! cell method and set cell method pointer
93  call method_cell_ptb%init( &
94  fmi=this%fmi, &
95  cell=this%cell, &
96  trackctl=this%trackctl, &
97  tracktimes=this%tracktimes)
98  submethod => method_cell_ptb
99  else if (particle%ifrctrn > 0) then
100  ! Force the ternary method
101  call method_cell_tern%init( &
102  fmi=this%fmi, &
103  cell=this%cell, &
104  trackctl=this%trackctl, &
105  tracktimes=this%tracktimes)
106  submethod => method_cell_tern
107  else if (cell%defn%can_be_rect) then
108  ! Cell is a rectangle, convert it to a rectangular cell type and
109  ! initialize Pollock's method
110  call cell_poly_to_rect(cell, rect)
111  base => rect
112  call method_cell_plck%init( &
113  fmi=this%fmi, &
114  cell=base, &
115  trackctl=this%trackctl, &
116  tracktimes=this%tracktimes)
117  submethod => method_cell_plck
118  else if (cell%defn%can_be_quad) then
119  ! Cell is quad-refined, convert to a quad rect cell type and
120  ! initialize the corresponding method
121  call cell_poly_to_quad(cell, quad)
122  base => quad
123  call method_cell_quad%init( &
124  fmi=this%fmi, &
125  cell=base, &
126  trackctl=this%trackctl, &
127  tracktimes=this%tracktimes)
128  submethod => method_cell_quad
129  else
130  ! Default to the ternary method
131  call method_cell_tern%init( &
132  fmi=this%fmi, &
133  cell=this%cell, &
134  trackctl=this%trackctl, &
135  tracktimes=this%tracktimes)
136  submethod => method_cell_tern
137  end if
138  end select
139  end subroutine load_disv
140 
141  subroutine load_particle(this, cell, particle)
142  ! modules
143  use disvmodule, only: disvtype
144  use particlemodule, only: term_boundary
145  ! dummy
146  class(methoddisvtype), intent(inout) :: this
147  type(cellpolytype), pointer, intent(inout) :: cell
148  type(particletype), pointer, intent(inout) :: particle
149  ! local
150  integer(I4B) :: inface
151  integer(I4B) :: ipos
152  integer(I4B) :: ic
153  integer(I4B) :: icu
154  integer(I4B) :: inbr
155  integer(I4B) :: idiag
156  integer(I4B) :: icpl
157  integer(I4B) :: ilay
158  real(DP) :: z
159 
160  select type (dis => this%fmi%dis)
161  type is (disvtype)
162  inface = particle%iboundary(2)
163  idiag = dis%con%ia(cell%defn%icell)
164  inbr = cell%defn%facenbr(inface)
165  ipos = idiag + inbr
166  ic = dis%con%ja(ipos)
167  icu = dis%get_nodeuser(ic)
168  call get_jk(icu, dis%ncpl, dis%nlay, icpl, ilay)
169 
170  ! if returning to a cell through the bottom
171  ! face after previously leaving it through
172  ! that same face, we've entered a cycle
173  ! as can occur e.g. in wells. terminate
174  ! in the previous cell.
175  if (ic == particle%icp .and. inface == 7 .and. ilay < particle%ilay) then
176  particle%advancing = .false.
177  particle%idomain(2) = particle%icp
178  particle%istatus = term_boundary
179  particle%izone = particle%izp
180  call this%save(particle, reason=3)
181  return
182  else
183  particle%icp = particle%idomain(2)
184  particle%izp = particle%izone
185  end if
186 
187  particle%idomain(2) = ic
188  particle%icu = icu
189  particle%ilay = ilay
190 
191  z = particle%z
192  call this%map_neighbor(cell%defn, inface, z)
193 
194  particle%iboundary(2) = inface
195  particle%idomain(3:) = 0
196  particle%iboundary(3:) = 0
197  particle%z = z
198  end select
199 
200  end subroutine
201 
202  subroutine update_flowja(this, cell, particle)
203  ! dummy
204  class(methoddisvtype), intent(inout) :: this
205  type(cellpolytype), pointer, intent(inout) :: cell
206  type(particletype), pointer, intent(inout) :: particle
207  ! local
208  integer(I4B) :: inbr
209  integer(I4B) :: idiag
210  integer(I4B) :: ipos
211 
212  idiag = this%fmi%dis%con%ia(cell%defn%icell)
213  inbr = cell%defn%facenbr(particle%iboundary(2))
214  ipos = idiag + inbr
215 
216  ! leaving old cell
217  this%flowja(ipos) = this%flowja(ipos) - done
218 
219  ! entering new cell
220  this%flowja(this%fmi%dis%con%isym(ipos)) &
221  = this%flowja(this%fmi%dis%con%isym(ipos)) + done
222 
223  end subroutine update_flowja
224 
225  !> @brief Pass a particle to the next cell, if there is one
226  subroutine pass_disv(this, particle)
227  use particlemodule, only: term_boundary
228  ! dummy
229  class(methoddisvtype), intent(inout) :: this
230  type(particletype), pointer, intent(inout) :: particle
231  ! local
232  type(cellpolytype), pointer :: cell
233 
234  select type (c => this%cell)
235  type is (cellpolytype)
236  cell => c
237  ! If the entry face has no neighbors it's a
238  ! boundary face, so terminate the particle.
239  ! todo AMP: reconsider when multiple models supported
240  if (cell%defn%facenbr(particle%iboundary(2)) .eq. 0) then
241  particle%istatus = term_boundary
242  particle%advancing = .false.
243  call this%save(particle, reason=3)
244  else
245  ! Otherwise, load cell properties into the
246  ! particle. It may be marked to terminate.
247  call this%load_particle(cell, particle)
248  if (.not. particle%advancing) return
249 
250  ! Update intercell mass flows
251  call this%update_flowja(cell, particle)
252  end if
253  end select
254  end subroutine pass_disv
255 
256  !> @brief Map location on cell face to shared cell face of neighbor
257  subroutine map_neighbor(this, defn, inface, z)
258  ! dummy
259  class(methoddisvtype), intent(inout) :: this
260  type(celldefntype), pointer, intent(inout) :: defn
261  integer(I4B), intent(inout) :: inface
262  double precision, intent(inout) :: z
263  ! local
264  integer(I4B) :: icin
265  integer(I4B) :: npolyvertsin
266  integer(I4B) :: ic
267  integer(I4B) :: npolyverts
268  integer(I4B) :: inbr
269  integer(I4B) :: inbrnbr
270  integer(I4B) :: j
271  integer(I4B) :: m
272  real(DP) :: zrel
273  real(DP) :: topfrom
274  real(DP) :: botfrom
275  real(DP) :: top
276  real(DP) :: bot
277  real(DP) :: sat
278 
279  ! Map to shared cell face of neighbor
280  inbr = defn%facenbr(inface)
281  if (inbr .eq. 0) then
282  ! Exterior face; no neighbor to map to
283  ! todo AMP: reconsider when multiple models allowed
284  inface = -1
285  else
286  ! Load definition for neighbor cell (neighbor with shared face)
287  icin = defn%icell
288  j = this%fmi%dis%con%ia(icin)
289  ic = this%fmi%dis%con%ja(j + inbr)
290  call this%load_cell_defn(ic, this%neighbor)
291 
292  npolyvertsin = defn%npolyverts
293  npolyverts = this%neighbor%npolyverts
294  if (inface .eq. npolyvertsin + 2) then
295  ! Exits through bot, enters through top
296  inface = npolyverts + 3
297  else if (inface .eq. npolyvertsin + 3) then
298  ! Exits through top, enters through bot
299  inface = npolyverts + 2
300  else
301  ! Exits and enters through shared polygon face
302  j = this%fmi%dis%con%ia(ic)
303  do m = 1, npolyverts + 3
304  inbrnbr = this%neighbor%facenbr(m)
305  if (this%fmi%dis%con%ja(j + inbrnbr) .eq. icin) then
306  inface = m
307  exit
308  end if
309  end do
310  ! Map z between cells
311  topfrom = defn%top
312  botfrom = defn%bot
313  zrel = (z - botfrom) / (topfrom - botfrom)
314  top = this%fmi%dis%top(ic)
315  bot = this%fmi%dis%bot(ic)
316  sat = this%fmi%gwfsat(ic)
317  z = bot + zrel * sat * (top - bot)
318  end if
319  end if
320  end subroutine map_neighbor
321 
322  !> @brief Apply the DISV tracking method to a particle.
323  subroutine apply_disv(this, particle, tmax)
324  class(methoddisvtype), intent(inout) :: this
325  type(particletype), pointer, intent(inout) :: particle
326  real(DP), intent(in) :: tmax
327 
328  call this%track(particle, 1, tmax)
329  end subroutine apply_disv
330 
331  !> @brief Load cell definition from the grid
332  subroutine load_cell_defn(this, ic, defn)
333  ! dummy
334  class(methoddisvtype), intent(inout) :: this
335  integer(I4B), intent(in) :: ic
336  type(celldefntype), pointer, intent(inout) :: defn
337 
338  call this%load_properties(ic, defn)
339  call this%load_polygon(defn)
340  call this%load_neighbors(defn)
341  call this%load_indicators(defn)
342  call this%load_flows(defn)
343  end subroutine load_cell_defn
344 
345  !> @brief Loads cell properties to cell definition from the grid.
346  subroutine load_properties(this, ic, defn)
347  ! dummy
348  class(methoddisvtype), intent(inout) :: this
349  integer(I4B), intent(in) :: ic
350  type(celldefntype), pointer, intent(inout) :: defn
351  ! local
352  real(DP) :: top
353  real(DP) :: bot
354  real(DP) :: sat
355  integer(I4B) :: icu, icpl, ilay
356 
357  defn%icell = ic
358  defn%iatop = get_iatop(this%fmi%dis%get_ncpl(), &
359  this%fmi%dis%get_nodeuser(ic))
360  top = this%fmi%dis%top(ic)
361  bot = this%fmi%dis%bot(ic)
362  sat = this%fmi%gwfsat(ic)
363  top = bot + sat * (top - bot)
364  defn%top = top
365  defn%bot = bot
366  defn%sat = sat
367  defn%porosity = this%porosity(ic)
368  defn%retfactor = this%retfactor(ic)
369  defn%izone = this%izone(ic)
370  select type (dis => this%fmi%dis)
371  type is (disvtype)
372  icu = dis%get_nodeuser(ic)
373  call get_jk(icu, dis%ncpl, dis%nlay, icpl, ilay)
374  defn%ilay = ilay
375  end select
376  end subroutine load_properties
377 
378  subroutine load_polygon(this, defn)
379  ! dummy
380  class(methoddisvtype), intent(inout) :: this
381  type(celldefntype), pointer, intent(inout) :: defn
382 
383  call this%fmi%dis%get_polyverts( &
384  defn%icell, &
385  defn%polyvert, &
386  closed=.true.)
387  defn%npolyverts = size(defn%polyvert, dim=2) - 1
388  end subroutine load_polygon
389 
390  !> @brief Loads face neighbors to cell definition from the grid
391  !! Assumes cell index and number of vertices are already loaded.
392  subroutine load_neighbors(this, defn)
393  ! dummy
394  class(methoddisvtype), intent(inout) :: this
395  type(celldefntype), pointer, intent(inout) :: defn
396  ! local
397  integer(I4B) :: ic1
398  integer(I4B) :: ic2
399  integer(I4B) :: icu1
400  integer(I4B) :: icu2
401  integer(I4B) :: j1
402  integer(I4B) :: j2
403  integer(I4B) :: k1
404  integer(I4B) :: k2
405  integer(I4B) :: iloc
406  integer(I4B) :: ipos
407  integer(I4B) :: istart1
408  integer(I4B) :: istart2
409  integer(I4B) :: istop1
410  integer(I4B) :: istop2
411  integer(I4B) :: isharedface
412  integer(I4B) :: ncpl
413  integer(I4B) :: nfaces
414  integer(I4B) :: nslots
415 
416  ! expand facenbr array if needed
417  nfaces = defn%npolyverts + 3
418  nslots = size(defn%facenbr)
419  if (nslots < nfaces) call expandarray(defn%facenbr, nfaces - nslots)
420 
421  select type (dis => this%fmi%dis)
422  type is (disvtype)
423  ! Load face neighbors.
424  defn%facenbr = 0
425  ic1 = defn%icell
426  icu1 = dis%get_nodeuser(ic1)
427  ncpl = dis%get_ncpl()
428  call get_jk(icu1, ncpl, dis%nlay, j1, k1)
429  istart1 = dis%iavert(j1)
430  istop1 = dis%iavert(j1 + 1) - 1
431  do iloc = 1, dis%con%ia(ic1 + 1) - dis%con%ia(ic1) - 1
432  ipos = dis%con%ia(ic1) + iloc
433  ! mask could become relevant if PRT uses interface model
434  if (dis%con%mask(ipos) == 0) cycle
435  ic2 = dis%con%ja(ipos)
436  icu2 = dis%get_nodeuser(ic2)
437  call get_jk(icu2, ncpl, dis%nlay, j2, k2)
438  istart2 = dis%iavert(j2)
439  istop2 = dis%iavert(j2 + 1) - 1
440  call shared_face(dis%javert(istart1:istop1), &
441  dis%javert(istart2:istop2), &
442  isharedface)
443  if (isharedface /= 0) then
444  ! Edge (polygon) face neighbor
445  defn%facenbr(isharedface) = int(iloc, 1)
446  else
447  if (k2 > k1) then
448  ! Bottom face neighbor
449  defn%facenbr(defn%npolyverts + 2) = int(iloc, 1)
450  else if (k2 < k1) then
451  ! Top face neighbor
452  defn%facenbr(defn%npolyverts + 3) = int(iloc, 1)
453  else
454  call pstop(1, "k2 should be <> k1, since no shared edge face")
455  end if
456  end if
457  end do
458  end select
459  ! List of edge (polygon) faces wraps around
460  defn%facenbr(defn%npolyverts + 1) = defn%facenbr(1)
461  end subroutine load_neighbors
462 
463  !> @brief Load flows into the cell definition.
464  !! These include face, boundary and net distributed flows.
465  !! Assumes cell index and number of vertices are already loaded.
466  subroutine load_flows(this, defn)
467  ! dummy
468  class(methoddisvtype), intent(inout) :: this
469  type(celldefntype), intent(inout) :: defn
470  ! local
471  integer(I4B) :: nfaces
472  integer(I4B) :: nslots
473 
474  ! expand faceflow array if needed
475  nfaces = defn%npolyverts + 3
476  nslots = size(defn%faceflow)
477  if (nslots < nfaces) call expandarray(defn%faceflow, nfaces - nslots)
478 
479  ! Load face flows, including boundary flows. As with cell verts,
480  ! the face flow array wraps around. Top and bottom flows make up
481  ! the last two elements, respectively, for size npolyverts + 3.
482  ! If there is no flow through any face, set a no-exit-face flag.
483  defn%faceflow = dzero
484  defn%inoexitface = 1
485  call this%load_boundary_flows_to_defn_poly(defn)
486  call this%load_face_flows_to_defn_poly(defn)
487 
488  ! Add up net distributed flow
489  defn%distflow = this%fmi%SourceFlows(defn%icell) + &
490  this%fmi%SinkFlows(defn%icell) + &
491  this%fmi%StorageFlows(defn%icell)
492 
493  ! Set weak sink flag
494  if (this%fmi%SinkFlows(defn%icell) .ne. dzero) then
495  defn%iweaksink = 1
496  else
497  defn%iweaksink = 0
498  end if
499  end subroutine load_flows
500 
501  subroutine load_face_flows_to_defn_poly(this, defn)
502  ! dummy
503  class(methoddisvtype), intent(inout) :: this
504  type(celldefntype), intent(inout) :: defn
505  ! local
506  integer(I4B) :: m, n, nfaces
507  real(DP) :: q
508 
509  nfaces = defn%npolyverts + 3
510  do m = 1, nfaces
511  n = defn%facenbr(m)
512  if (n > 0) then
513  q = this%fmi%gwfflowja(this%fmi%dis%con%ia(defn%icell) + n)
514  defn%faceflow(m) = defn%faceflow(m) + q
515  end if
516  if (defn%faceflow(m) < dzero) defn%inoexitface = 0
517  end do
518  end subroutine load_face_flows_to_defn_poly
519 
520  !> @brief Load boundary flows from the grid into a rectangular cell.
521  !! Assumes cell index and number of vertices are already loaded.
522  subroutine load_boundary_flows_to_defn_rect(this, defn)
523  ! dummy
524  class(methoddisvtype), intent(inout) :: this
525  type(celldefntype), intent(inout) :: defn
526  ! local
527  integer(I4B) :: ioffset
528 
529  ! assignment of BoundaryFlows to faceflow below assumes clockwise
530  ! ordering of faces, with face 1 being the "western" face
531  ioffset = (defn%icell - 1) * 10
532  defn%faceflow(1) = defn%faceflow(1) + &
533  this%fmi%BoundaryFlows(ioffset + 4)
534  defn%faceflow(2) = defn%faceflow(2) + &
535  this%fmi%BoundaryFlows(ioffset + 2)
536  defn%faceflow(3) = defn%faceflow(3) + &
537  this%fmi%BoundaryFlows(ioffset + 3)
538  defn%faceflow(4) = defn%faceflow(4) + &
539  this%fmi%BoundaryFlows(ioffset + 1)
540  defn%faceflow(5) = defn%faceflow(1)
541  defn%faceflow(6) = defn%faceflow(6) + &
542  this%fmi%BoundaryFlows(ioffset + 9)
543  defn%faceflow(7) = defn%faceflow(7) + &
544  this%fmi%BoundaryFlows(ioffset + 10)
545  end subroutine load_boundary_flows_to_defn_rect
546 
547  !> @brief Load boundary flows from the grid into rectangular quadcell.
548  !! Assumes cell index and number of vertices are already loaded.
550  ! dummy
551  class(methoddisvtype), intent(inout) :: this
552  type(celldefntype), intent(inout) :: defn
553  ! local
554  integer(I4B) :: m
555  integer(I4B) :: n
556  integer(I4B) :: nn
557  integer(I4B) :: ioffset
558  integer(I4B) :: nbf
559  integer(I4B) :: m1
560  integer(I4B) :: m2
561  integer(I4B) :: mdiff
562  real(DP) :: qbf
563  integer(I4B) :: irectvert(5)
564 
565  ioffset = (defn%icell - 1) * 10
566 
567  ! Polygon faces in positions 1 through npolyverts
568  do n = 1, 4
569  if (n .eq. 2) then
570  nbf = 4
571  else if (n .eq. 4) then
572  nbf = 1
573  else
574  nbf = n
575  end if
576  qbf = this%fmi%BoundaryFlows(ioffset + nbf)
577  nn = 0
578  do m = 1, defn%npolyverts
579  if (.not. defn%ispv180(m)) then
580  nn = nn + 1
581  irectvert(nn) = m
582  end if
583  end do
584  irectvert(5) = irectvert(1)
585  m1 = irectvert(n)
586  m2 = irectvert(n + 1)
587  if (m2 .lt. m1) m2 = m2 + defn%npolyverts
588  mdiff = m2 - m1
589  if (mdiff .eq. 1) then
590  ! Assign BoundaryFlow to corresponding polygon face
591  defn%faceflow(m1) = defn%faceflow(m1) + qbf
592  else
593  ! Split BoundaryFlow between two faces on quad-refined edge
594  qbf = 5d-1 * qbf
595  defn%faceflow(m1) = defn%faceflow(m1) + qbf
596  defn%faceflow(m1 + 1) = defn%faceflow(m1 + 1) + qbf
597  end if
598  end do
599  ! Wrap around to 1 in position npolyverts+1
600  m = defn%npolyverts + 1
601  defn%faceflow(m) = defn%faceflow(1)
602  ! Bottom in position npolyverts+2
603  m = m + 1
604  defn%faceflow(m) = defn%faceflow(m) + &
605  this%fmi%BoundaryFlows(ioffset + 9)
606  ! Top in position npolyverts+3
607  m = m + 1
608  defn%faceflow(m) = defn%faceflow(m) + &
609  this%fmi%BoundaryFlows(ioffset + 10)
610 
612 
613  !> @brief Load boundary flows from the grid into a polygonal cell.
614  !! Assumes cell index and number of vertices are already loaded.
615  subroutine load_boundary_flows_to_defn_poly(this, defn)
616  ! dummy
617  class(methoddisvtype), intent(inout) :: this
618  type(celldefntype), intent(inout) :: defn
619  ! local
620  integer(I4B) :: ic
621  integer(I4B) :: npolyverts
622  integer(I4B) :: ioffset
623  integer(I4B) :: iv
624 
625  ic = defn%icell
626  npolyverts = defn%npolyverts
627 
628  ioffset = (ic - 1) * max_poly_cells
629  do iv = 1, npolyverts
630  defn%faceflow(iv) = &
631  defn%faceflow(iv) + &
632  this%fmi%BoundaryFlows(ioffset + iv)
633  end do
634  defn%faceflow(npolyverts + 1) = defn%faceflow(1)
635  defn%faceflow(npolyverts + 2) = &
636  defn%faceflow(npolyverts + 2) + &
637  this%fmi%BoundaryFlows(ioffset + max_poly_cells - 1)
638  defn%faceflow(npolyverts + 3) = &
639  defn%faceflow(npolyverts + 3) + &
640  this%fmi%BoundaryFlows(ioffset + max_poly_cells)
641 
642  end subroutine load_boundary_flows_to_defn_poly
643 
644  !> @brief Load 180-degree vertex indicator array and set flags
645  !! indicating how cell can be represented. Assumes cell index
646  !! and number of vertices are already loaded.
647  !<
648  subroutine load_indicators(this, defn)
649  ! dummy
650  class(methoddisvtype), intent(inout) :: this
651  type(celldefntype), pointer, intent(inout) :: defn
652  ! local
653  integer(I4B) :: npolyverts
654  integer(I4B) :: m
655  integer(I4B) :: m0
656  integer(I4B) :: m1
657  integer(I4B) :: m2
658  integer(I4B) :: ic
659  integer(I4B) :: num90
660  integer(I4B) :: num180
661  real(DP) :: x0
662  real(DP) :: y0
663  real(DP) :: x1
664  real(DP) :: y1
665  real(DP) :: x2
666  real(DP) :: y2
667  real(DP) :: epsang
668  real(DP) :: s0x
669  real(DP) :: s0y
670  real(DP) :: &
671  s0mag, s2x, s2y, s2mag, sinang, cosang, dotprod
672  logical(LGP) last180
673 
674  ic = defn%icell
675  npolyverts = defn%npolyverts
676 
677  if (size(defn%ispv180) < npolyverts + 3) &
678  call expandarray(defn%ispv180, npolyverts + 1)
679 
680  defn%ispv180(1:npolyverts + 1) = .false.
681  defn%can_be_rect = .false.
682  defn%can_be_quad = .false.
683  epsang = 1d-5
684  num90 = 0
685  num180 = 0
686  last180 = .false.
687 
688  ! assumes non-self-intersecting polygon
689  do m = 1, npolyverts
690  m1 = m
691  if (m1 .eq. 1) then
692  m0 = npolyverts
693  m2 = 2
694  else if (m1 .eq. npolyverts) then
695  m0 = npolyverts - 1
696  m2 = 1
697  else
698  m0 = m1 - 1
699  m2 = m1 + 1
700  end if
701  x0 = defn%polyvert(1, m0)
702  y0 = defn%polyvert(2, m0)
703  x1 = defn%polyvert(1, m1)
704  y1 = defn%polyvert(2, m1)
705  x2 = defn%polyvert(1, m2)
706  y2 = defn%polyvert(2, m2)
707  s0x = x0 - x1
708  s0y = y0 - y1
709  s0mag = dsqrt(s0x * s0x + s0y * s0y)
710  s2x = x2 - x1
711  s2y = y2 - y1
712  s2mag = dsqrt(s2x * s2x + s2y * s2y)
713  sinang = (s0x * s2y - s0y * s2x) / (s0mag * s2mag)
714  cosang = dsqrt(dabs(done - (sinang * sinang)))
715  if (dabs(sinang) .lt. epsang) then
716  dotprod = s0x * s2x + s0y * s2y
717  if (dotprod .lt. dzero) then
718  num180 = num180 + 1
719  last180 = .true.
720  defn%ispv180(m) = .true.
721  end if
722  else
723  if (dabs(cosang) .lt. epsang) num90 = num90 + 1
724  last180 = .false.
725  end if
726  end do
727 
728  ! List of 180-degree indicators wraps around for convenience
729  defn%ispv180(npolyverts + 1) = defn%ispv180(1)
730 
731  ! Set rect/quad flags
732  if (num90 .eq. 4) then
733  if (num180 .eq. 0) then
734  defn%can_be_rect = .true.
735  else
736  defn%can_be_quad = .true.
737  end if
738  end if
739  end subroutine load_indicators
740 
741 end module methoddisvmodule
subroutine, public create_defn(cellDefn)
Create a new cell definition object.
Definition: CellDefn.f90:42
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_poly(cell)
Create a new polygonal cell.
Definition: CellPoly.f90:20
subroutine, public cell_poly_to_rect(poly, rect)
Convert CellPoly representation to CellRect. Assumes the conversion is possible.
Definition: CellUtil.f90:14
subroutine, public cell_poly_to_quad(poly, quad)
Convert CellPoly representation to CellRectQuad. Assumes the conversion is possible.
Definition: CellUtil.f90:115
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
subroutine pstop(status, message)
Stop the program, optionally specifying an error status code.
Definition: ErrorUtil.f90:24
subroutine, public shared_face(iverts1, iverts2, iface)
Find the lateral face shared by two cells.
Definition: GeomUtil.f90:403
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
Cell-level tracking methods.
type(methodcellpollocktype), pointer, public method_cell_plck
type(methodcellpasstobottype), pointer, public method_cell_ptb
type(methodcellpollockquadtype), pointer, public method_cell_quad
type(methodcellternarytype), pointer, public method_cell_tern
subroutine, public create_method_disv(method)
Create a new vertex grid (DISV) tracking method.
Definition: MethodDisv.f90:49
subroutine load_boundary_flows_to_defn_rect(this, defn)
Load boundary flows from the grid into a rectangular cell. Assumes cell index and number of vertices ...
Definition: MethodDisv.f90:523
subroutine load_neighbors(this, defn)
Loads face neighbors to cell definition from the grid Assumes cell index and number of vertices are a...
Definition: MethodDisv.f90:393
subroutine map_neighbor(this, defn, inface, z)
Map location on cell face to shared cell face of neighbor.
Definition: MethodDisv.f90:258
subroutine load_indicators(this, defn)
Load 180-degree vertex indicator array and set flags indicating how cell can be represented....
Definition: MethodDisv.f90:649
subroutine load_boundary_flows_to_defn_rect_quad(this, defn)
Load boundary flows from the grid into rectangular quadcell. Assumes cell index and number of vertice...
Definition: MethodDisv.f90:550
subroutine update_flowja(this, cell, particle)
Definition: MethodDisv.f90:203
subroutine load_disv(this, particle, next_level, submethod)
Load the cell and the tracking method.
Definition: MethodDisv.f90:71
subroutine load_face_flows_to_defn_poly(this, defn)
Definition: MethodDisv.f90:502
subroutine load_polygon(this, defn)
Definition: MethodDisv.f90:379
subroutine load_particle(this, cell, particle)
Definition: MethodDisv.f90:142
subroutine load_properties(this, ic, defn)
Loads cell properties to cell definition from the grid.
Definition: MethodDisv.f90:347
subroutine load_flows(this, defn)
Load flows into the cell definition. These include face, boundary and net distributed flows....
Definition: MethodDisv.f90:467
subroutine load_boundary_flows_to_defn_poly(this, defn)
Load boundary flows from the grid into a polygonal cell. Assumes cell index and number of vertices ar...
Definition: MethodDisv.f90:616
subroutine load_cell_defn(this, ic, defn)
Load cell definition from the grid.
Definition: MethodDisv.f90:333
subroutine apply_disv(this, particle, tmax)
Apply the DISV tracking method to a particle.
Definition: MethodDisv.f90:324
subroutine pass_disv(this, particle)
Pass a particle to the next cell, if there is one.
Definition: MethodDisv.f90:227
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
Base type for grid cells of a concrete type. Contains a cell-definition which is information shared b...
Definition: Cell.f90:13
Vertex grid discretization.
Definition: Disv.f90:24
Base type for particle tracking methods.
Definition: Method.f90:31
Particle tracked by the PRT model.
Definition: Particle.f90:78