MODFLOW 6 version 6.9.0.dev0
USGS Modular Hydrologic Model
Loading...
Searching...
No Matches
Data Types | Functions/Subroutines
disvgeom Module Reference

Data Types

type  disvgeomtype

Functions/Subroutines

subroutine init (this, nlay, ncpl, nodes, top_grid, bot_grid, iavert, javert, vertex_grid, cellxy_grid, nodereduced, nodeuser)
 Initialize.
subroutine set_kj (this, k, j)
 Set node IDs.
subroutine set_nodered (this, nodered)
 Set reduced node number.
subroutine cell_setup (this)
 Set top and bottom elevations of grid cell.
subroutine cprops (this, cell2, hwva, cl1, cl2, ax, ihc)
subroutine edge_normal (this, cell2, xcomp, ycomp)
 Return the x and y components of an outward normal facing vector.
subroutine connection_vector (this, cell2, nozee, satn, satm, xcomp, ycomp, zcomp, conlen)
 Return the x y and z components of a unit vector that points from from the center of this to the center of cell2, and the straight-line connection length.
logical function shares_edge (this, cell2)
 Return true if this shares a horizontal edge with cell2.
subroutine, public shared_edge (ivlist1, ivlist2, ivert1, ivert2)
 Find two common vertices shared by cell1 and cell2.
real(dp) function get_area (this)
 Calculate and return the area of the cell.
real(dp) function anglex (x1, y1, x2, y2)
 Calculate the angle that the x-axis makes with a line that is normal to the two points.
real(dp) function distance (x1, y1, x2, y2)
 Calculate distance between two points.
real(dp) function distance_normal (x0, y0, x1, y1, x2, y2)
 Calculate normal distance from point (x0, y0) to line defined by two points, (x1, y1), (x2, y2).
subroutine line_unit_normal (x0, y0, x1, y1, xcomp, ycomp)
 Calculate the normal vector components (xcomp and ycomp) for a line defined by two points, (x0, y0), (x1, y1).
subroutine, public line_unit_vector (x0, y0, z0, x1, y1, z1, xcomp, ycomp, zcomp, vmag)
 Calculate the vector components (xcomp, ycomp, and zcomp) for a line defined by two points, (x0, y0, z0), (x1, y1, z1).

Function/Subroutine Documentation

◆ anglex()

real(dp) function disvgeom::anglex ( real(dp), intent(in) x1,
real(dp), intent(in) y1,
real(dp), intent(in) x2,
real(dp), intent(in) y2 )
private

This assumes that vertices are numbered clockwise so that the angle is for the normal outward of cell n.

Definition at line 394 of file DisvGeom.f90.

395 ! -- modules
396 use constantsmodule, only: dzero, dtwo, dpi
397 ! -- dummy
398 real(DP), intent(in) :: x1
399 real(DP), intent(in) :: x2
400 real(DP), intent(in) :: y1
401 real(DP), intent(in) :: y2
402 ! -- return
403 real(DP) :: ax
404 ! -- local
405 real(DP) :: dx
406 real(DP) :: dy
407 !
408 dx = x2 - x1
409 dy = y2 - y1
410 ax = atan2(dx, -dy)
411 if (ax < dzero) ax = dtwo * dpi + ax
constantsmodule
This module contains simulation constants.
Definition Constants.f90:9
constantsmodule::dpi
real(dp), parameter dpi
real constant
Definition Constants.f90:128
constantsmodule::dzero
real(dp), parameter dzero
real constant zero
Definition Constants.f90:65
constantsmodule::dtwo
real(dp), parameter dtwo
real constant 2
Definition Constants.f90:79
Here is the caller graph for this function:

◆ cell_setup()

subroutine disvgeom::cell_setup ( class(disvgeomtype) this)
private

Definition at line 128 of file DisvGeom.f90.

129 ! -- dummy
130 class(DisvGeomType) :: this
131 !
132 this%top = this%top_grid(this%nodered)
133 this%bot = this%bot_grid(this%nodered)

◆ connection_vector()

subroutine disvgeom::connection_vector ( class(disvgeomtype) this,
type(disvgeomtype) cell2,
logical, intent(in) nozee,
real(dp), intent(in) satn,
real(dp), intent(in) satm,
real(dp), intent(out) xcomp,
real(dp), intent(out) ycomp,
real(dp), intent(out) zcomp,
real(dp), intent(out) conlen )
private

Definition at line 238 of file DisvGeom.f90.

240 ! -- module
241 use constantsmodule, only: dzero, dhalf, done
242 ! -- dummy
243 class(DisvGeomType) :: this
244 type(DisvGeomType) :: cell2
245 logical, intent(in) :: nozee
246 real(DP), intent(in) :: satn
247 real(DP), intent(in) :: satm
248 real(DP), intent(out) :: xcomp
249 real(DP), intent(out) :: ycomp
250 real(DP), intent(out) :: zcomp
251 real(DP), intent(out) :: conlen
252 ! -- local
253 real(DP) :: x1, y1, z1, x2, y2, z2
254 !
255 x1 = this%cellxy_grid(1, this%j)
256 y1 = this%cellxy_grid(2, this%j)
257 x2 = this%cellxy_grid(1, cell2%j)
258 y2 = this%cellxy_grid(2, cell2%j)
259 if (nozee) then
260 z1 = dzero
261 z2 = dzero
262 else
263 z1 = this%bot + dhalf * satn * (this%top - this%bot)
264 z2 = cell2%bot + dhalf * satm * (cell2%top - cell2%bot)
265 end if
266 !
267 call line_unit_vector(x1, y1, z1, x2, y2, z2, xcomp, ycomp, zcomp, &
268 conlen)
constantsmodule::dhalf
real(dp), parameter dhalf
real constant 1/2
Definition Constants.f90:68
constantsmodule::done
real(dp), parameter done
real constant 1
Definition Constants.f90:76
Here is the call graph for this function:

◆ cprops()

subroutine disvgeom::cprops ( class(disvgeomtype) this,
type(disvgeomtype) cell2,
real(dp), intent(out) hwva,
real(dp), intent(out) cl1,
real(dp), intent(out) cl2,
real(dp), intent(out) ax,
integer(i4b), intent(out) ihc )
private

Definition at line 136 of file DisvGeom.f90.

137 ! -- module
138 use constantsmodule, only: dzero, dhalf, done
139 ! -- dummy
140 class(DisvGeomType) :: this
141 type(DisvGeomType) :: cell2
142 real(DP), intent(out) :: hwva
143 real(DP), intent(out) :: cl1
144 real(DP), intent(out) :: cl2
145 real(DP), intent(out) :: ax
146 integer(I4B), intent(out) :: ihc
147 ! -- local
148 integer(I4B) :: ivert1, ivert2
149 integer(I4B) :: istart1, istart2, istop1, istop2
150 real(DP) :: x0, y0, x1, y1, x2, y2
151 !
152 if (this%j == cell2%j) then
153 !
154 ! -- Cells share same j index, so must be a vertical connection
155 ihc = 0
156 hwva = this%get_area()
157 cl1 = dhalf * (this%top - this%bot)
158 cl2 = dhalf * (cell2%top - cell2%bot)
159 ax = dzero
160 else
161 !
162 ! -- Must be horizontal connection
163 ihc = 1
164 istart1 = this%iavert(this%j)
165 istop1 = this%iavert(this%j + 1) - 1
166 istart2 = cell2%iavert(cell2%j)
167 istop2 = this%iavert(cell2%j + 1) - 1
168 call shared_edge(this%javert(istart1:istop1), &
169 this%javert(istart2:istop2), &
170 ivert1, ivert2)
171 if (ivert1 == 0 .or. ivert2 == 0) then
172 !
173 ! -- Cells do not share an edge
174 hwva = dzero
175 cl1 = done
176 cl2 = done
177 else
178 x1 = this%vertex_grid(1, ivert1)
179 y1 = this%vertex_grid(2, ivert1)
180 x2 = this%vertex_grid(1, ivert2)
181 y2 = this%vertex_grid(2, ivert2)
182 hwva = distance(x1, y1, x2, y2)
183 !
184 ! -- cl1
185 x0 = this%cellxy_grid(1, this%j)
186 y0 = this%cellxy_grid(2, this%j)
187 cl1 = distance_normal(x0, y0, x1, y1, x2, y2)
188 !
189 ! -- cl2
190 x0 = this%cellxy_grid(1, cell2%j)
191 y0 = this%cellxy_grid(2, cell2%j)
192 cl2 = distance_normal(x0, y0, x1, y1, x2, y2)
193 !
194 ! -- anglex
195 x1 = this%vertex_grid(1, ivert1)
196 y1 = this%vertex_grid(2, ivert1)
197 x2 = this%vertex_grid(1, ivert2)
198 y2 = this%vertex_grid(2, ivert2)
199 ax = anglex(x1, y1, x2, y2)
200 end if
201 end if
Here is the call graph for this function:

◆ distance()

real(dp) function disvgeom::distance ( real(dp), intent(in) x1,
real(dp), intent(in) y1,
real(dp), intent(in) x2,
real(dp), intent(in) y2 )
private

Definition at line 416 of file DisvGeom.f90.

417 ! -- dummy
418 real(DP), intent(in) :: x1
419 real(DP), intent(in) :: x2
420 real(DP), intent(in) :: y1
421 real(DP), intent(in) :: y2
422 ! -- return
423 real(DP) :: d
424 !
425 d = (x1 - x2)**2 + (y1 - y2)**2
426 d = sqrt(d)
Here is the caller graph for this function:

◆ distance_normal()

real(dp) function disvgeom::distance_normal ( real(dp), intent(in) x0,
real(dp), intent(in) y0,
real(dp), intent(in) x1,
real(dp), intent(in) y1,
real(dp), intent(in) x2,
real(dp), intent(in) y2 )
private

Definition at line 432 of file DisvGeom.f90.

433 ! -- dummy
434 real(DP), intent(in) :: x0
435 real(DP), intent(in) :: y0
436 real(DP), intent(in) :: x1
437 real(DP), intent(in) :: y1
438 real(DP), intent(in) :: x2
439 real(DP), intent(in) :: y2
440 ! -- return
441 real(DP) :: d
442 !
443 d = abs((x2 - x1) * (y1 - y0) - (x1 - x0) * (y2 - y1))
444 d = d / distance(x1, y1, x2, y2)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ edge_normal()

subroutine disvgeom::edge_normal ( class(disvgeomtype) this,
type(disvgeomtype) cell2,
real(dp), intent(out) xcomp,
real(dp), intent(out) ycomp )
private

Definition at line 206 of file DisvGeom.f90.

207 ! -- module
208 use constantsmodule, only: dzero, dhalf, done
209 ! -- dummy
210 class(DisvGeomType) :: this
211 type(DisvGeomType) :: cell2
212 real(DP), intent(out) :: xcomp
213 real(DP), intent(out) :: ycomp
214 ! -- local
215 integer(I4B) :: ivert1, ivert2
216 integer(I4B) :: istart1, istart2, istop1, istop2
217 real(DP) :: x1, y1, x2, y2
218 !
219 istart1 = this%iavert(this%j)
220 istop1 = this%iavert(this%j + 1) - 1
221 istart2 = cell2%iavert(cell2%j)
222 istop2 = this%iavert(cell2%j + 1) - 1
223 call shared_edge(this%javert(istart1:istop1), &
224 this%javert(istart2:istop2), &
225 ivert1, ivert2)
226 x1 = this%vertex_grid(1, ivert1)
227 y1 = this%vertex_grid(2, ivert1)
228 x2 = this%vertex_grid(1, ivert2)
229 y2 = this%vertex_grid(2, ivert2)
230 !
231 call line_unit_normal(x1, y1, x2, y2, xcomp, ycomp)
Here is the call graph for this function:

◆ get_area()

real(dp) function disvgeom::get_area ( class(disvgeomtype) this)
private

a = 1/2 *[(x1*y2 + x2*y3 + x3*y4 + ... + xn*y1) - (x2*y1 + x3*y2 + x4*y3 + ... + x1*yn)]

Definition at line 339 of file DisvGeom.f90.

340 ! -- module
341 use constantsmodule, only: dzero, dhalf
342 ! -- dummy
343 class(DisvGeomType) :: this
344 ! -- return
345 real(DP) :: area
346 ! -- local
347 integer(I4B) :: ivert
348 integer(I4B) :: nvert
349 integer(I4B) :: icount
350 integer(I4B) :: iv1
351 real(DP) :: x
352 real(DP) :: y
353 real(DP) :: x1
354 real(DP) :: y1
355 !
356 area = dzero
357 nvert = this%iavert(this%j + 1) - this%iavert(this%j)
358 icount = 1
359 iv1 = this%javert(this%iavert(this%j))
360 x1 = this%vertex_grid(1, iv1)
361 y1 = this%vertex_grid(2, iv1)
362 do ivert = this%iavert(this%j), this%iavert(this%j + 1) - 1
363 x = this%vertex_grid(1, this%javert(ivert))
364 if (icount < nvert) then
365 y = this%vertex_grid(2, this%javert(ivert + 1))
366 else
367 y = this%vertex_grid(2, this%javert(this%iavert(this%j)))
368 end if
369 area = area + (x - x1) * (y - y1)
370 icount = icount + 1
371 end do
372 !
373 icount = 1
374 do ivert = this%iavert(this%j), this%iavert(this%j + 1) - 1
375 y = this%vertex_grid(2, this%javert(ivert))
376 if (icount < nvert) then
377 x = this%vertex_grid(1, this%javert(ivert + 1))
378 else
379 x = this%vertex_grid(1, this%javert(this%iavert(this%j)))
380 end if
381 area = area - (x - x1) * (y - y1)
382 icount = icount + 1
383 end do
384 !
385 area = abs(area) * dhalf
Here is the call graph for this function:

◆ init()

subroutine disvgeom::init ( class(disvgeomtype) this,
integer(i4b), intent(in) nlay,
integer(i4b), intent(in) ncpl,
integer(i4b), intent(in) nodes,
real(dp), dimension(nodes), target top_grid,
real(dp), dimension(nodes), target bot_grid,
integer(i4b), dimension(:), target iavert,
integer(i4b), dimension(:), target javert,
real(dp), dimension(:, :), target vertex_grid,
real(dp), dimension(:, :), target cellxy_grid,
integer(i4b), dimension(ncpl, nlay), target nodereduced,
integer(i4b), dimension(nodes), target nodeuser )

Definition at line 50 of file DisvGeom.f90.

52 ! -- dummy
53 class(DisvGeomType) :: this
54 integer(I4B), intent(in) :: nlay
55 integer(I4B), intent(in) :: ncpl
56 integer(I4B), intent(in) :: nodes
57 real(DP), dimension(nodes), target :: top_grid
58 real(DP), dimension(nodes), target :: bot_grid
59 integer(I4B), dimension(:), target :: iavert
60 integer(I4B), dimension(:), target :: javert
61 real(DP), dimension(:, :), target :: vertex_grid
62 real(DP), dimension(:, :), target :: cellxy_grid
63 integer(I4B), dimension(ncpl, nlay), target :: nodereduced
64 integer(I4B), dimension(nodes), target :: nodeuser
65 ! -- local
66 integer(I4B) :: nodesuser
67 !
68 this%nlay = nlay
69 this%ncpl = ncpl
70 this%nodes = nodes
71 this%top_grid => top_grid
72 this%bot_grid => bot_grid
73 this%iavert => iavert
74 this%javert => javert
75 this%vertex_grid => vertex_grid
76 this%cellxy_grid => cellxy_grid
77 this%nodereduced => nodereduced
78 this%nodeuser => nodeuser
79 nodesuser = ncpl * nlay
80 !
81 if (nodes < nodesuser) then
82 this%reduced = .true.
83 else
84 this%reduced = .false.
85 end if

◆ line_unit_normal()

subroutine disvgeom::line_unit_normal ( real(dp), intent(in) x0,
real(dp), intent(in) y0,
real(dp), intent(in) x1,
real(dp), intent(in) y1,
real(dp), intent(out) xcomp,
real(dp), intent(out) ycomp )
private

Definition at line 450 of file DisvGeom.f90.

451 ! -- dummy
452 real(DP), intent(in) :: x0
453 real(DP), intent(in) :: y0
454 real(DP), intent(in) :: x1
455 real(DP), intent(in) :: y1
456 real(DP), intent(out) :: xcomp
457 real(DP), intent(out) :: ycomp
458 ! -- local
459 real(DP) :: dx, dy, vmag
460 !
461 dx = x1 - x0
462 dy = y1 - y0
463 vmag = sqrt(dx**2 + dy**2)
464 xcomp = -dy / vmag
465 ycomp = dx / vmag
Here is the caller graph for this function:

◆ line_unit_vector()

subroutine, public disvgeom::line_unit_vector ( real(dp), intent(in) x0,
real(dp), intent(in) y0,
real(dp), intent(in) z0,
real(dp), intent(in) x1,
real(dp), intent(in) y1,
real(dp), intent(in) z1,
real(dp), intent(out) xcomp,
real(dp), intent(out) ycomp,
real(dp), intent(out) zcomp,
real(dp) vmag )

Also return the magnitude of the original vector, vmag.

Definition at line 473 of file DisvGeom.f90.

475 ! -- dummy
476 real(DP), intent(in) :: x0
477 real(DP), intent(in) :: y0
478 real(DP), intent(in) :: z0
479 real(DP), intent(in) :: x1
480 real(DP), intent(in) :: y1
481 real(DP), intent(in) :: z1
482 real(DP), intent(out) :: xcomp
483 real(DP), intent(out) :: ycomp
484 real(DP), intent(out) :: zcomp
485 real(DP) :: vmag
486 ! -- local
487 real(DP) :: dx, dy, dz
488 !
489 dx = x1 - x0
490 dy = y1 - y0
491 dz = z1 - z0
492 vmag = sqrt(dx**2 + dy**2 + dz**2)
493 xcomp = dx / vmag
494 ycomp = dy / vmag
495 zcomp = dz / vmag
Here is the caller graph for this function:

◆ set_kj()

subroutine disvgeom::set_kj ( class(disvgeomtype) this,
integer(i4b), intent(in) k,
integer(i4b), intent(in) j )
private

Definition at line 90 of file DisvGeom.f90.

91 ! -- dummy
92 class(DisvGeomType) :: this
93 integer(I4B), intent(in) :: k
94 integer(I4B), intent(in) :: j
95 !
96 this%k = k
97 this%j = j
98 this%nodeusr = get_node(k, 1, j, this%nlay, 1, this%ncpl)
99 if (this%reduced) then
100 this%nodered = this%nodereduced(k, j)
101 else
102 this%nodered = this%nodeusr
103 end if
104 call this%cell_setup()
Here is the call graph for this function:
Here is the caller graph for this function:

◆ set_nodered()

subroutine disvgeom::set_nodered ( class(disvgeomtype) this,
integer(i4b), intent(in) nodered )
private

Definition at line 109 of file DisvGeom.f90.

110 ! -- dummy
111 class(DisvGeomType) :: this
112 integer(I4B), intent(in) :: nodered
113 !
114 this%nodered = nodered
115 !
116 if (this%reduced) then
117 this%nodeusr = this%nodeuser(nodered)
118 else
119 this%nodeusr = nodered
120 end if
121 !
122 call get_jk(this%nodeusr, this%ncpl, this%nlay, this%j, this%k)
123 call this%cell_setup()
Here is the call graph for this function:
Here is the caller graph for this function:

◆ shared_edge()

subroutine, public disvgeom::shared_edge ( integer(i4b), dimension(:) ivlist1,
integer(i4b), dimension(:) ivlist2,
integer(i4b), intent(out) ivert1,
integer(i4b), intent(out) ivert2 )

Return 0 if there are no shared edges. Proceed forward through ivlist1 and backward through ivlist2 as a clockwise face in cell1 must correspond to a counter clockwise face in cell2.

Definition at line 302 of file DisvGeom.f90.

303 ! -- dummy
304 integer(I4B), dimension(:) :: ivlist1
305 integer(I4B), dimension(:) :: ivlist2
306 integer(I4B), intent(out) :: ivert1
307 integer(I4B), intent(out) :: ivert2
308 ! -- local
309 integer(I4B) :: nv1
310 integer(I4B) :: nv2
311 integer(I4B) :: il1
312 integer(I4B) :: il2
313 logical :: found
314 !
315 found = .false.
316 nv1 = size(ivlist1)
317 nv2 = size(ivlist2)
318 ivert1 = 0
319 ivert2 = 0
320 outerloop: do il1 = 1, nv1 - 1
321 do il2 = nv2, 2, -1
322 if (ivlist1(il1) == ivlist2(il2) .and. &
323 ivlist1(il1 + 1) == ivlist2(il2 - 1)) then
324 found = .true.
325 ivert1 = ivlist1(il1)
326 ivert2 = ivlist1(il1 + 1)
327 exit outerloop
328 end if
329 end do
330 if (found) exit
331 end do outerloop
Here is the caller graph for this function:

◆ shares_edge()

logical function disvgeom::shares_edge ( class(disvgeomtype) this,
type(disvgeomtype) cell2 )
private

Definition at line 273 of file DisvGeom.f90.

274 ! -- dummy
275 class(DisvGeomType) :: this
276 type(DisvGeomType) :: cell2
277 ! -- return
278 logical l
279 ! -- local
280 integer(I4B) :: istart1, istop1, istart2, istop2
281 integer(I4B) :: ivert1, ivert2
282 !
283 istart1 = this%iavert(this%j)
284 istop1 = this%iavert(this%j + 1) - 1
285 istart2 = cell2%iavert(cell2%j)
286 istop2 = this%iavert(cell2%j + 1) - 1
287 call shared_edge(this%javert(istart1:istop1), &
288 this%javert(istart2:istop2), &
289 ivert1, ivert2)
290 l = .true.
291 if (ivert1 == 0 .or. ivert2 == 0) then
292 l = .false.
293 end if
Here is the call graph for this function: