55 allocate (method%name)
59 method%delegates = .true.
66 deallocate (this%name)
70 subroutine load_disv(this, particle, next_level, submethod)
78 integer(I4B),
intent(in) :: next_level
79 class(
methodtype),
pointer,
intent(inout) :: submethod
86 select type (cell => this%cell)
88 ic = particle%idomain(next_level)
89 call this%load_cell_defn(ic, cell%defn)
90 if (this%fmi%ibdgwfsat0(ic) == 0)
then
96 trackctl=this%trackctl, &
97 tracktimes=this%tracktimes)
99 else if (particle%ifrctrn > 0)
then
104 trackctl=this%trackctl, &
105 tracktimes=this%tracktimes)
107 else if (cell%defn%can_be_rect)
then
115 trackctl=this%trackctl, &
116 tracktimes=this%tracktimes)
118 else if (cell%defn%can_be_quad)
then
126 trackctl=this%trackctl, &
127 tracktimes=this%tracktimes)
134 trackctl=this%trackctl, &
135 tracktimes=this%tracktimes)
150 integer(I4B) :: inface
155 integer(I4B) :: idiag
160 select type (dis => this%fmi%dis)
162 inface = particle%iboundary(2)
163 idiag = dis%con%ia(cell%defn%icell)
164 inbr = cell%defn%facenbr(inface)
166 ic = dis%con%ja(ipos)
167 icu = dis%get_nodeuser(ic)
168 call get_jk(icu, dis%ncpl, dis%nlay, icpl, ilay)
175 if (ic == particle%icp .and. inface == 7 .and. ilay < particle%ilay)
then
176 particle%advancing = .false.
177 particle%idomain(2) = particle%icp
179 particle%izone = particle%izp
180 call this%save(particle, reason=3)
183 particle%icp = particle%idomain(2)
184 particle%izp = particle%izone
187 particle%idomain(2) = ic
192 call this%map_neighbor(cell%defn, inface, z)
194 particle%iboundary(2) = inface
195 particle%idomain(3:) = 0
196 particle%iboundary(3:) = 0
209 integer(I4B) :: idiag
212 idiag = this%fmi%dis%con%ia(cell%defn%icell)
213 inbr = cell%defn%facenbr(particle%iboundary(2))
217 this%flowja(ipos) = this%flowja(ipos) -
done
220 this%flowja(this%fmi%dis%con%isym(ipos)) &
221 = this%flowja(this%fmi%dis%con%isym(ipos)) +
done
234 select type (c => this%cell)
240 if (cell%defn%facenbr(particle%iboundary(2)) .eq. 0)
then
242 particle%advancing = .false.
243 call this%save(particle, reason=3)
247 call this%load_particle(cell, particle)
248 if (.not. particle%advancing)
return
251 call this%update_flowja(cell, particle)
261 integer(I4B),
intent(inout) :: inface
262 double precision,
intent(inout) :: z
265 integer(I4B) :: npolyvertsin
267 integer(I4B) :: npolyverts
269 integer(I4B) :: inbrnbr
280 inbr = defn%facenbr(inface)
281 if (inbr .eq. 0)
then
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)
292 npolyvertsin = defn%npolyverts
293 npolyverts = this%neighbor%npolyverts
294 if (inface .eq. npolyvertsin + 2)
then
296 inface = npolyverts + 3
297 else if (inface .eq. npolyvertsin + 3)
then
299 inface = npolyverts + 2
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
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)
326 real(DP),
intent(in) :: tmax
328 call this%track(particle, 1, tmax)
335 integer(I4B),
intent(in) :: ic
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)
349 integer(I4B),
intent(in) :: ic
355 integer(I4B) :: icu, icpl, ilay
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)
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)
372 icu = dis%get_nodeuser(ic)
373 call get_jk(icu, dis%ncpl, dis%nlay, icpl, ilay)
383 call this%fmi%dis%get_polyverts( &
387 defn%npolyverts =
size(defn%polyvert, dim=2) - 1
407 integer(I4B) :: istart1
408 integer(I4B) :: istart2
409 integer(I4B) :: istop1
410 integer(I4B) :: istop2
411 integer(I4B) :: isharedface
413 integer(I4B) :: nfaces
414 integer(I4B) :: nslots
417 nfaces = defn%npolyverts + 3
418 nslots =
size(defn%facenbr)
419 if (nslots < nfaces)
call expandarray(defn%facenbr, nfaces - nslots)
421 select type (dis => this%fmi%dis)
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
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
441 dis%javert(istart2:istop2), &
443 if (isharedface /= 0)
then
445 defn%facenbr(isharedface) = int(iloc, 1)
449 defn%facenbr(defn%npolyverts + 2) = int(iloc, 1)
450 else if (k2 < k1)
then
452 defn%facenbr(defn%npolyverts + 3) = int(iloc, 1)
454 call pstop(1,
"k2 should be <> k1, since no shared edge face")
460 defn%facenbr(defn%npolyverts + 1) = defn%facenbr(1)
471 integer(I4B) :: nfaces
472 integer(I4B) :: nslots
475 nfaces = defn%npolyverts + 3
476 nslots =
size(defn%faceflow)
477 if (nslots < nfaces)
call expandarray(defn%faceflow, nfaces - nslots)
483 defn%faceflow =
dzero
485 call this%load_boundary_flows_to_defn_poly(defn)
486 call this%load_face_flows_to_defn_poly(defn)
489 defn%distflow = this%fmi%SourceFlows(defn%icell) + &
490 this%fmi%SinkFlows(defn%icell) + &
491 this%fmi%StorageFlows(defn%icell)
494 if (this%fmi%SinkFlows(defn%icell) .ne.
dzero)
then
506 integer(I4B) :: m, n, nfaces
509 nfaces = defn%npolyverts + 3
513 q = this%fmi%gwfflowja(this%fmi%dis%con%ia(defn%icell) + n)
514 defn%faceflow(m) = defn%faceflow(m) + q
516 if (defn%faceflow(m) <
dzero) defn%inoexitface = 0
527 integer(I4B) :: ioffset
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)
557 integer(I4B) :: ioffset
561 integer(I4B) :: mdiff
563 integer(I4B) :: irectvert(5)
565 ioffset = (defn%icell - 1) * 10
571 else if (n .eq. 4)
then
576 qbf = this%fmi%BoundaryFlows(ioffset + nbf)
578 do m = 1, defn%npolyverts
579 if (.not. defn%ispv180(m))
then
584 irectvert(5) = irectvert(1)
586 m2 = irectvert(n + 1)
587 if (m2 .lt. m1) m2 = m2 + defn%npolyverts
589 if (mdiff .eq. 1)
then
591 defn%faceflow(m1) = defn%faceflow(m1) + qbf
595 defn%faceflow(m1) = defn%faceflow(m1) + qbf
596 defn%faceflow(m1 + 1) = defn%faceflow(m1 + 1) + qbf
600 m = defn%npolyverts + 1
601 defn%faceflow(m) = defn%faceflow(1)
604 defn%faceflow(m) = defn%faceflow(m) + &
605 this%fmi%BoundaryFlows(ioffset + 9)
608 defn%faceflow(m) = defn%faceflow(m) + &
609 this%fmi%BoundaryFlows(ioffset + 10)
621 integer(I4B) :: npolyverts
622 integer(I4B) :: ioffset
626 npolyverts = defn%npolyverts
629 do iv = 1, npolyverts
630 defn%faceflow(iv) = &
631 defn%faceflow(iv) + &
632 this%fmi%BoundaryFlows(ioffset + iv)
634 defn%faceflow(npolyverts + 1) = defn%faceflow(1)
635 defn%faceflow(npolyverts + 2) = &
636 defn%faceflow(npolyverts + 2) + &
638 defn%faceflow(npolyverts + 3) = &
639 defn%faceflow(npolyverts + 3) + &
653 integer(I4B) :: npolyverts
659 integer(I4B) :: num90
660 integer(I4B) :: num180
671 s0mag, s2x, s2y, s2mag, sinang, cosang, dotprod
675 npolyverts = defn%npolyverts
677 if (
size(defn%ispv180) < npolyverts + 3) &
680 defn%ispv180(1:npolyverts + 1) = .false.
681 defn%can_be_rect = .false.
682 defn%can_be_quad = .false.
694 else if (m1 .eq. npolyverts)
then
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)
709 s0mag = dsqrt(s0x * s0x + s0y * s0y)
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
720 defn%ispv180(m) = .true.
723 if (dabs(cosang) .lt. epsang) num90 = num90 + 1
729 defn%ispv180(npolyverts + 1) = defn%ispv180(1)
732 if (num90 .eq. 4)
then
733 if (num180 .eq. 0)
then
734 defn%can_be_rect = .true.
736 defn%can_be_quad = .true.
subroutine, public create_defn(cellDefn)
Create a new cell definition object.
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 ...
integer(i4b), parameter, public max_poly_cells
subroutine, public create_cell_poly(cell)
Create a new polygonal cell.
subroutine, public cell_poly_to_rect(poly, rect)
Convert CellPoly representation to CellRect. Assumes the conversion is possible.
subroutine, public cell_poly_to_quad(poly, quad)
Convert CellPoly representation to CellRectQuad. Assumes the conversion is possible.
This module contains simulation constants.
real(dp), parameter dzero
real constant zero
real(dp), parameter done
real constant 1
subroutine pstop(status, message)
Stop the program, optionally specifying an error status code.
subroutine, public shared_face(iverts1, iverts2, iface)
Find the lateral face shared by two cells.
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...
This module defines variable data types.
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.
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 ...
subroutine load_neighbors(this, defn)
Loads face neighbors to cell definition from the grid Assumes cell index and number of vertices are a...
subroutine map_neighbor(this, defn, inface, z)
Map location on cell face to shared cell face of neighbor.
subroutine load_indicators(this, defn)
Load 180-degree vertex indicator array and set flags indicating how cell can be represented....
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...
subroutine update_flowja(this, cell, particle)
subroutine load_disv(this, particle, next_level, submethod)
Load the cell and the tracking method.
subroutine load_face_flows_to_defn_poly(this, defn)
subroutine load_polygon(this, defn)
subroutine load_particle(this, cell, particle)
subroutine load_properties(this, ic, defn)
Loads cell properties to cell definition from the grid.
subroutine load_flows(this, defn)
Load flows into the cell definition. These include face, boundary and net distributed flows....
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...
subroutine load_cell_defn(this, ic, defn)
Load cell definition from the grid.
subroutine apply_disv(this, particle, tmax)
Apply the DISV tracking method to a particle.
subroutine pass_disv(this, particle)
Pass a particle to the next cell, if there is one.
Particle tracking strategies.
@ term_boundary
terminated at a boundary face
Base grid cell definition.
Base type for grid cells of a concrete type. Contains a cell-definition which is information shared b...
Vertex grid discretization.
Base type for particle tracking methods.
Particle tracked by the PRT model.