14 real(dp),
intent(in) :: x, a, b
15 between = ((x >= a .and. x <= b) .or. (x <= a .and. x >= b))
25 real(dp),
intent(in) :: x
26 real(dp),
intent(in) :: y
27 real(dp),
allocatable,
intent(in) :: poly(:, :)
29 integer(I4B) :: i, ii, num_verts
30 real(dp) :: xa, xb, ya, yb, c = 0.0_dp
33 num_verts =
size(poly, 2)
34 xa = poly(1, num_verts)
35 ya = poly(2, num_verts)
37 do i = 0, num_verts - 1
38 ii = mod(i, num_verts) + 1
42 if ((x == xa .and. y == ya) .or. &
43 (x == xb .and. y == yb))
then
47 else if (ya == yb .and. &
53 else if (
between(y, ya, yb))
then
54 if ((y == ya .and. yb >= ya) .or. &
55 (y == yb .and. ya >= yb))
then
61 c = (xa - x) * (yb - y) - (xb - x) * (ya - y)
66 else if ((ya < yb) .eqv. (c > 0))
then
79 function get_node(ilay, irow, icol, nlay, nrow, ncol)
80 integer(I4B),
intent(in) :: ilay, irow, icol, nlay, nrow, ncol
83 if (nlay > 0 .and. nrow > 0 .and. ncol > 0 .and. &
84 ilay > 0 .and. ilay <= nlay .and. &
85 irow > 0 .and. irow <= nrow .and. &
86 icol > 0 .and. icol <= ncol)
then
88 icol + ncol * (irow - 1) + (ilay - 1) * nrow * ncol
96 subroutine get_ijk(nodenumber, nrow, ncol, nlay, irow, icol, ilay)
98 integer(I4B),
intent(in) :: nodenumber
99 integer(I4B),
intent(in) :: nrow
100 integer(I4B),
intent(in) :: ncol
101 integer(I4B),
intent(in) :: nlay
102 integer(I4B),
intent(out) :: irow
103 integer(I4B),
intent(out) :: icol
104 integer(I4B),
intent(out) :: ilay
106 integer(I4B) :: nodes
109 nodes = nlay * nrow * ncol
110 if (nodenumber < 1 .or. nodenumber > nodes)
then
115 ilay = (nodenumber - 1) / (ncol * nrow) + 1
116 ij = nodenumber - (ilay - 1) * ncol * nrow
117 irow = (ij - 1) / ncol + 1
118 icol = ij - (irow - 1) * ncol
124 subroutine get_jk(nodenumber, ncpl, nlay, icpl, ilay)
126 integer(I4B),
intent(in) :: nodenumber
127 integer(I4B),
intent(in) :: ncpl
128 integer(I4B),
intent(in) :: nlay
129 integer(I4B),
intent(out) :: icpl
130 integer(I4B),
intent(out) :: ilay
132 integer(I4B) :: nodes
135 if (nodenumber < 1 .or. nodenumber > nodes)
then
139 ilay = (nodenumber - 1) / ncpl + 1
140 icpl = nodenumber - (ilay - 1) * ncpl
145 pure function skew(v, s, invert)
result(res)
147 real(dp),
intent(in) :: v(2)
148 real(dp),
intent(in) :: s(3)
149 logical(LGP),
intent(in),
optional :: invert
152 logical(LGP) :: linvert
153 real(dp) :: sxx, sxy, syy
156 if (
present(invert))
then
165 if (.not. linvert)
then
166 res(1) = sxx * v(1) + sxy * v(2)
170 res(1) = (v(1) - sxy * res(2)) / sxx
177 xorigin, yorigin, zorigin, &
181 real(dp) :: xin, yin, zin
182 real(dp) :: xout, yout, zout
183 real(dp),
optional :: xorigin, yorigin, zorigin
184 real(dp),
optional :: sinrot, cosrot
185 logical(LGP),
optional :: invert
187 logical(LGP) :: ltranslate, lrotate, linvert
189 real(dp) :: lxorigin, lyorigin, lzorigin
190 real(dp) :: lsinrot, lcosrot
193 call defaults(lxorigin, lyorigin, lzorigin, &
194 lsinrot, lcosrot, linvert, &
195 ltranslate, lrotate, &
196 xorigin, yorigin, zorigin, &
197 sinrot, cosrot, invert)
200 if (.not. linvert)
then
203 xout = xin - lxorigin
204 yout = yin - lyorigin
205 zout = zin - lzorigin
214 xout = x * lcosrot + y * lsinrot
215 yout = -x * lsinrot + y * lcosrot
220 x = xin * lcosrot - yin * lsinrot
221 y = xin * lsinrot + yin * lcosrot
229 zout = zin + lzorigin
235 subroutine compose(xorigin, yorigin, zorigin, &
237 xorigin_new, yorigin_new, zorigin_new, &
238 sinrot_new, cosrot_new, &
241 real(dp) :: xorigin, yorigin, zorigin
242 real(dp) :: sinrot, cosrot
243 real(dp),
optional :: xorigin_new, yorigin_new, zorigin_new
244 real(dp),
optional :: sinrot_new, cosrot_new
245 logical(LGP),
optional :: invert
247 logical(LGP) :: ltranslate, lrotate, linvert
248 real(dp) :: xorigin_add, yorigin_add, zorigin_add
249 real(dp) :: sinrot_add, cosrot_add
250 real(dp) :: x0, y0, z0, s0, c0
253 call defaults(xorigin_add, yorigin_add, zorigin_add, &
254 sinrot_add, cosrot_add, linvert, &
255 ltranslate, lrotate, &
256 xorigin_new, yorigin_new, zorigin_new, &
257 sinrot_new, cosrot_new, invert)
267 if (.not. linvert)
then
274 call transform(xorigin_add, yorigin_add, zorigin_add, &
275 xorigin, yorigin, zorigin, &
276 x0, y0, z0, s0, c0, .true.)
282 sinrot = cosrot_add * s0 + sinrot_add * c0
283 cosrot = cosrot_add * c0 - sinrot_add * s0
293 call transform(-xorigin_add, -yorigin_add, zorigin_add, &
294 x0, y0, z0, xorigin, yorigin, zorigin, &
295 -sinrot_add, cosrot_add, .true.)
297 xorigin = c0 * x0 - s0 * y0
298 yorigin = s0 * x0 + c0 * y0
305 sinrot = cosrot_add * s0 - sinrot_add * c0
306 cosrot = cosrot_add * c0 + sinrot_add * s0
314 invert, translate, rotate, &
315 xorigin_opt, yorigin_opt, zorigin_opt, &
316 sinrot_opt, cosrot_opt, invert_opt)
318 real(DP) :: xorigin, yorigin, zorigin
319 real(DP) :: sinrot, cosrot
320 logical(LGP) :: invert, translate, rotate
321 real(DP),
optional :: xorigin_opt, yorigin_opt, zorigin_opt
322 real(DP),
optional :: sinrot_opt, cosrot_opt
323 logical(LGP),
optional :: invert_opt
327 if (
present(xorigin_opt))
then
328 xorigin = xorigin_opt
332 if (
present(yorigin_opt))
then
333 yorigin = yorigin_opt
337 if (
present(zorigin_opt))
then
338 zorigin = zorigin_opt
344 if (
present(sinrot_opt))
then
346 if (
present(cosrot_opt))
then
351 cosrot = dsqrt(
done - sinrot * sinrot)
354 else if (
present(cosrot_opt))
then
358 sinrot = dsqrt(
done - cosrot * cosrot)
362 if (
present(invert_opt)) invert = invert_opt
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 transform(xin, yin, zin, xout, yout, zout, xorigin, yorigin, zorigin, sinrot, cosrot, invert)
Apply a 3D translation and optional 2D rotation to coordinates.
integer(i4b) function, public get_node(ilay, irow, icol, nlay, nrow, ncol)
Get node number, given layer, row, and column indices for a structured grid. If any argument is inval...
subroutine defaults(xorigin, yorigin, zorigin, sinrot, cosrot, invert, translate, rotate, xorigin_opt, yorigin_opt, zorigin_opt, sinrot_opt, cosrot_opt, invert_opt)
Process arguments and set defaults. Internal use only.
logical function, public between(x, a, b)
Check if a value is between two other values (inclusive).
subroutine, public compose(xorigin, yorigin, zorigin, sinrot, cosrot, xorigin_new, yorigin_new, zorigin_new, sinrot_new, cosrot_new, invert)
Apply a 3D translation and 2D rotation to an existing transformation.
logical function, public point_in_polygon(x, y, poly)
Check if a point is within a polygon.
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,...
pure real(dp) function, dimension(2), public skew(v, s, invert)
Skew a 2D vector along the x-axis.
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.