27 procedure,
public :: pass =>
pass_dis
50 allocate (method%name)
54 method%delegates = .true.
60 deallocate (this%name)
66 integer(I4B),
intent(in) :: ic
82 select type (dis => this%fmi%dis)
84 icu = dis%get_nodeuser(ic)
86 call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, &
91 dz = cell%defn%top - cell%defn%bot
98 cell%xOrigin = cell%defn%polyvert(1, 1)
99 cell%yOrigin = cell%defn%polyvert(2, 1)
100 cell%zOrigin = cell%defn%bot
103 factor =
done / cell%defn%retfactor
104 factor = factor / cell%defn%porosity
107 term = factor / areaz
109 cell%vz1 = cell%defn%faceflow(6) * term
110 cell%vz2 = -cell%defn%faceflow(7) * term
112 if (this%fmi%ibdgwfsat0(ic) == 0)
then
123 term = factor / areax
124 cell%vx1 = cell%defn%faceflow(1) * term
125 cell%vx2 = -cell%defn%faceflow(3) * term
128 term = factor / areay
129 cell%vy1 = cell%defn%faceflow(4) * term
130 cell%vy2 = -cell%defn%faceflow(2) * term
136 subroutine load_dis(this, particle, next_level, submethod)
140 integer(I4B),
intent(in) :: next_level
141 class(
methodtype),
pointer,
intent(inout) :: submethod
145 select type (cell => this%cell)
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
154 trackctl=this%trackctl, &
155 tracktimes=this%tracktimes)
161 trackctl=this%trackctl, &
162 tracktimes=this%tracktimes)
182 integer(I4B) :: inface
183 integer(I4B) :: idiag
194 select type (dis => this%fmi%dis)
197 inface = particle%iboundary(2)
198 inbr = cell%defn%facenbr(inface)
199 idiag = this%fmi%dis%con%ia(cell%defn%icell)
201 ic = dis%con%ja(ipos)
202 icu = dis%get_nodeuser(ic)
203 call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, &
211 if (ic == particle%icp .and. inface == 7 .and. ilay < particle%ilay)
then
212 particle%advancing = .false.
213 particle%idomain(2) = particle%icp
215 particle%izone = particle%izp
216 call this%save(particle, reason=3)
219 particle%icp = particle%idomain(2)
220 particle%izp = particle%izone
224 particle%idomain(2) = ic
229 if (inface .eq. 1)
then
231 else if (inface .eq. 2)
then
233 else if (inface .eq. 3)
then
235 else if (inface .eq. 4)
then
237 else if (inface .eq. 6)
then
239 else if (inface .eq. 7)
then
242 particle%iboundary(2) = inface
247 topfrom = cell%defn%top
248 botfrom = cell%defn%bot
249 zrel = (z - botfrom) / (topfrom - botfrom)
252 sat = this%fmi%gwfsat(ic)
253 z = bot + zrel * sat * (top - bot)
267 integer(I4B) :: idiag
269 integer(I4B) :: inface
272 inface = particle%iboundary(2)
273 inbr = cell%defn%facenbr(inface)
274 idiag = this%fmi%dis%con%ia(cell%defn%icell)
278 this%flowja(ipos) = this%flowja(ipos) -
done
281 this%flowja(this%fmi%dis%con%isym(ipos)) &
282 = this%flowja(this%fmi%dis%con%isym(ipos)) +
done
294 select type (c => this%cell)
300 if (cell%defn%facenbr(particle%iboundary(2)) .eq. 0)
then
302 particle%advancing = .false.
303 call this%save(particle, reason=3)
306 call this%load_particle(cell, particle)
307 if (.not. particle%advancing)
return
310 call this%update_flowja(cell, particle)
319 real(DP),
intent(in) :: tmax
321 call this%track(particle, 1, tmax)
327 integer,
intent(in) :: iatop
328 double precision :: top
330 if (iatop .lt. 0)
then
331 top = this%fmi%dis%top(-iatop)
333 top = this%fmi%dis%bot(iatop)
341 integer(I4B),
intent(in) :: ic
345 call this%load_properties(ic, defn)
348 call this%fmi%dis%get_polyverts( &
352 call this%load_neighbors(defn)
355 defn%ispv180(1:defn%npolyverts + 1) = .false.
357 call this%load_flows(defn)
364 integer(I4B),
intent(in) :: ic
367 integer(I4B) :: irow, icol, ilay, icu
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)
382 icu = dis%get_nodeuser(ic)
383 call get_ijk(icu, dis%nrow, dis%ncol, dis%nlay, irow, icol, ilay)
386 defn%izone = this%izone(ic)
387 defn%can_be_rect = .true.
388 defn%can_be_quad = .false.
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
413 select type (dis => this%fmi%dis)
418 icu1 = dis%get_nodeuser(ic1)
419 call get_ijk(icu1, dis%nrow, dis%ncol, dis%nlay, &
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
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, &
430 if (klay2 == klay1)
then
432 if (irow2 > irow1)
then
435 else if (jcol2 > jcol1)
then
438 else if (irow2 < irow1)
then
445 defn%facenbr(iedgeface) = int(iloc, 1)
446 else if (klay2 > klay1)
then
448 defn%facenbr(defn%npolyverts + 2) = int(iloc, 1)
451 defn%facenbr(defn%npolyverts + 3) = int(iloc, 1)
456 defn%facenbr(defn%npolyverts + 1) = defn%facenbr(1)
470 defn%faceflow =
dzero
472 call this%load_boundary_flows_to_defn(defn)
473 call this%load_face_flows_to_defn(defn)
476 defn%distflow = this%fmi%SourceFlows(defn%icell) + &
477 this%fmi%SinkFlows(defn%icell) + &
478 this%fmi%StorageFlows(defn%icell)
481 if (this%fmi%SinkFlows(defn%icell) .ne.
dzero)
then
493 integer(I4B) :: m, n, nfaces
496 nfaces = defn%npolyverts + 3
500 q = this%fmi%gwfflowja(this%fmi%dis%con%ia(defn%icell) + n)
501 defn%faceflow(m) = defn%faceflow(m) + q
503 if (defn%faceflow(m) <
dzero) defn%inoexitface = 0
514 integer(I4B) :: ioffset
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)
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_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.
@ 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.