28 procedure,
public :: pass =>
pass_dis
51 allocate (method%name)
55 method%delegates = .true.
61 deallocate (this%name)
67 integer(I4B),
intent(in) :: ic
83 select type (dis => this%fmi%dis)
85 icu = dis%get_nodeuser(ic)
87 call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, &
92 dz = cell%defn%top - cell%defn%bot
99 cell%xOrigin = cell%defn%polyvert(1, 1)
100 cell%yOrigin = cell%defn%polyvert(2, 1)
101 cell%zOrigin = cell%defn%bot
104 factor =
done / cell%defn%retfactor
105 factor = factor / cell%defn%porosity
108 term = factor / areaz
110 cell%vz1 = cell%defn%faceflow(6) * term
111 cell%vz2 = -cell%defn%faceflow(7) * term
113 if (this%fmi%ibdgwfsat0(ic) == 0)
then
124 term = factor / areax
125 cell%vx1 = cell%defn%faceflow(1) * term
126 cell%vx2 = -cell%defn%faceflow(3) * term
129 term = factor / areay
130 cell%vy1 = cell%defn%faceflow(4) * term
131 cell%vy2 = -cell%defn%faceflow(2) * term
137 subroutine load_dis(this, particle, next_level, submethod)
141 integer(I4B),
intent(in) :: next_level
142 class(
methodtype),
pointer,
intent(inout) :: submethod
146 select type (cell => this%cell)
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
155 events=this%events, &
156 tracktimes=this%tracktimes)
162 events=this%events, &
163 tracktimes=this%tracktimes)
184 integer(I4B) :: inface
185 integer(I4B) :: idiag
196 select type (dis => this%fmi%dis)
200 inbr = cell%defn%facenbr(inface)
201 idiag = this%fmi%dis%con%ia(cell%defn%icell)
203 ic = dis%con%ja(ipos)
204 icu = dis%get_nodeuser(ic)
205 call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, &
213 if (ic == particle%icp .and. inface == 7 .and. ilay < particle%ilay)
then
215 particle%izone = particle%izp
216 call this%terminate(particle, &
221 particle%izp = particle%izone
230 if (inface .eq. 1)
then
232 else if (inface .eq. 2)
then
234 else if (inface .eq. 3)
then
236 else if (inface .eq. 4)
then
238 else if (inface .eq. 6)
then
240 else if (inface .eq. 7)
then
248 topfrom = cell%defn%top
249 botfrom = cell%defn%bot
250 zrel = (z - botfrom) / (topfrom - botfrom)
253 sat = this%fmi%gwfsat(ic)
254 z = bot + zrel * sat * (top - bot)
268 integer(I4B) :: idiag
270 integer(I4B) :: inface
274 inbr = cell%defn%facenbr(inface)
275 idiag = this%fmi%dis%con%ia(cell%defn%icell)
279 this%flowja(ipos) = this%flowja(ipos) -
done
282 this%flowja(this%fmi%dis%con%isym(ipos)) &
283 = this%flowja(this%fmi%dis%con%isym(ipos)) +
done
296 select type (c => this%cell)
302 if (cell%defn%facenbr(particle%iboundary(
level_feature)) .eq. 0)
then
303 call this%terminate(particle, &
307 call this%load_particle(cell, particle)
308 if (.not. particle%advancing)
return
311 call this%update_flowja(cell, particle)
320 real(DP),
intent(in) :: tmax
322 call this%track(particle, 1, tmax)
328 integer,
intent(in) :: iatop
329 double precision :: top
331 if (iatop .lt. 0)
then
332 top = this%fmi%dis%top(-iatop)
334 top = this%fmi%dis%bot(iatop)
342 integer(I4B),
intent(in) :: ic
346 call this%load_properties(ic, defn)
349 call this%fmi%dis%get_polyverts( &
353 call this%load_neighbors(defn)
356 defn%ispv180(1:defn%npolyverts + 1) = .false.
358 call this%load_flows(defn)
365 integer(I4B),
intent(in) :: ic
368 integer(I4B) :: irow, icol, ilay, icu
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)
383 icu = dis%get_nodeuser(ic)
384 call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, irow, icol, ilay)
387 defn%izone = this%izone(ic)
388 defn%can_be_rect = .true.
389 defn%can_be_quad = .false.
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
414 select type (dis => this%fmi%dis)
419 icu1 = dis%get_nodeuser(ic1)
420 call get_ijk(icu1, dis%nrow, dis%ncol, dis%nlay, &
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
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, &
431 if (klay2 == klay1)
then
433 if (irow2 > irow1)
then
436 else if (jcol2 > jcol1)
then
439 else if (irow2 < irow1)
then
446 defn%facenbr(iedgeface) = int(iloc, 1)
447 else if (klay2 > klay1)
then
449 defn%facenbr(defn%npolyverts + 2) = int(iloc, 1)
452 defn%facenbr(defn%npolyverts + 3) = int(iloc, 1)
457 defn%facenbr(defn%npolyverts + 1) = defn%facenbr(1)
471 defn%faceflow =
dzero
473 call this%load_boundary_flows_to_defn(defn)
474 call this%load_face_flows_to_defn(defn)
477 defn%distflow = this%fmi%SourceFlows(defn%icell) + &
478 this%fmi%SinkFlows(defn%icell) + &
479 this%fmi%StorageFlows(defn%icell)
482 if (this%fmi%SinkFlows(defn%icell) .ne.
dzero)
then
494 integer(I4B) :: m, n, nfaces
497 nfaces = defn%npolyverts + 3
501 q = this%fmi%gwfflowja(this%fmi%dis%con%ia(defn%icell) + n)
502 defn%faceflow(m) = defn%faceflow(m) + q
504 if (defn%faceflow(m) <
dzero) defn%inoexitface = 0
515 integer(I4B) :: max_faces
516 integer(I4B) :: ioffset
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)
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 ...
subroutine, public create_cell_rect(cell)
Create a new rectangular cell.
This module contains simulation constants.
real(dp), parameter dzero
real constant zero
real(dp), parameter done
real constant 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,...
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.
pure logical function, public is_close(a, b, rtol, atol, symmetric)
Check if a real value is approximately equal to another.
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....
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...
double precision function get_top(this, iatop)
Returns a top elevation based on index iatop.
subroutine load_particle(this, cell, particle)
Load cell properties into the particle, including.
subroutine load_properties(this, ic, defn)
subroutine update_flowja(this, cell, particle)
Update cell-cell flows of particle mass. Every particle is currently assigned unit mass.
subroutine, public create_method_dis(method)
Create a new structured grid (DIS) tracking method.
subroutine load_celldefn(this, ic, defn)
Loads cell definition from the grid.
subroutine apply_dis(this, particle, tmax)
Apply the structured tracking method to a particle.
subroutine load_dis(this, particle, next_level, submethod)
Load the cell geometry and tracking method.
subroutine load_cell(this, ic, cell)
subroutine load_face_flows_to_defn(this, defn)
subroutine load_neighbors(this, defn)
Loads face neighbors to cell definition from the grid. Assumes cell index and number of vertices are ...
subroutine pass_dis(this, particle)
Pass a particle to the next cell, if there is one.
Particle tracking strategies.
@, public terminate
particle terminated
@ term_boundary
terminated at a boundary face
Base grid cell definition.
Structured grid discretization.
Base type for particle tracking methods.
Particle tracked by the PRT model.