MODFLOW 6  version 6.5.0.dev2
MODFLOW 6 Code Documentation
cellutilmodule Module Reference

Functions/Subroutines

subroutine, public cell_poly_to_rect (poly, rect)
 Convert CellPoly representation to CellRect if possible. More...
 
subroutine, public cell_poly_to_quad (poly, quad)
 Convert CellPoly representation to CellRectQuad if possible. More...
 

Function/Subroutine Documentation

◆ cell_poly_to_quad()

subroutine, public cellutilmodule::cell_poly_to_quad ( type(cellpolytype), intent(in), pointer  poly,
type(cellrectquadtype), intent(inout), pointer  quad 
)

Definition at line 114 of file CellUtil.f90.

116  use cellpolymodule, only: cellpolytype
117  use mathutilmodule, only: mod_offset
118  ! -- dummy
119  type(CellPolyType), intent(in), pointer :: poly
120  type(CellRectQuadType), intent(inout), pointer :: quad
121  ! -- local
122  integer :: i, irv, isc
123  double precision :: qhalf, qdisttopbot, q1, q2, q4
124 
125  call create_cell_rect_quad(quad)
126  call quad%init_from(poly%defn)
127  ! kluge note: no check whether conversion is possible; assumes it is
128  ! -- Translate and rotate the rect-quad cell into local coordinates with
129  ! -- x varying from 0 to dx and y varying from 0 to dy. Choose the "south-
130  ! -- west" rectangle vertex to be the local origin so that the rotation
131  ! -- angle is zero if the cell already aligns with the model x and y
132  ! -- coordinates.
133  quad%irvOrigin = quad%get_irectvertSW() ! kluge note: no need to pass all that stuff in call below -- set internally in CellRectQuad
134  call quad%get_rectDimensionsRotation( &
135  quad%irvOrigin, quad%xOrigin, &
136  quad%yOrigin, quad%zOrigin, &
137  quad%dx, quad%dy, &
138  quad%dz, quad%sinrot, &
139  quad%cosrot)
140 
141  ! -- Set the external and internal face flows used for subcells
142  do i = 0, 3
143  irv = mod_offset(i + quad%irvOrigin, 4, 1)
144  isc = mod_offset(i + 3, 4, 1)
145  if (.not. quad%face_is_refined(irv)) then
146  qhalf = 5d-1 * quad%get_rectflow(1, irv)
147  quad%qextl2(isc) = qhalf
148  isc = mod_offset(isc + 1, 4, 1)
149  quad%qextl1(isc) = qhalf
150  else
151  quad%qextl2(isc) = quad%get_rectflow(1, irv)
152  isc = mod_offset(isc + 1, 4, 1)
153  quad%qextl1(isc) = quad%get_rectflow(2, irv)
154  end if
155  end do
156  qdisttopbot = 2.5d-1 * (quad%defn%get_distflow() &
157  + quad%defn%get_botflow() &
158  + quad%defn%get_topflow())
159  q1 = qdisttopbot + quad%qextl1(1) + quad%qextl2(1)
160  q2 = qdisttopbot + quad%qextl1(2) + quad%qextl2(2)
161  q4 = qdisttopbot + quad%qextl1(4) + quad%qextl2(4)
162  quad%qintl(1) = -5d-1 * (q1 + 5d-1 * (q2 - q4))
163  quad%qintl(2) = quad%qintl(1) + q1
164  quad%qintl(3) = quad%qintl(2) + q2
165  quad%qintl(4) = quad%qintl(1) - q4
166  quad%qintl(5) = quad%qintl(1)
subroutine, public create_cell_rect_quad(cell)
Create a new rectangular-quad cell.
Here is the call graph for this function:
Here is the caller graph for this function:

◆ cell_poly_to_rect()

subroutine, public cellutilmodule::cell_poly_to_rect ( type(cellpolytype), intent(in), pointer  poly,
type(cellrecttype), intent(inout), pointer  rect 
)

Definition at line 12 of file CellUtil.f90.

13  use constantsmodule, only: done
15  use cellpolymodule, only: cellpolytype
16  use celldefnmodule, only: celldefntype
17  ! -- dummy
18  type(CellPolyType), intent(in), pointer :: poly
19  type(CellRectType), intent(inout), pointer :: rect
20  ! -- local
21  type(CellDefnType), pointer :: defn
22  integer :: ipv, ipv1, ipv2, ipv3, ipv4
23  integer, dimension(4) :: ipvnxt = (/2, 3, 4, 1/)
24  double precision :: x1, y1, x2, y2, x4, y4
25  double precision :: dx2, dy2, dx4, dy4, areax, areay, areaz
26  double precision :: xOrigin, yOrigin, zOrigin, dx, dy, dz, sinrot, cosrot
27  double precision :: factor, term
28 
29  call create_cell_rect(rect)
30  defn => poly%defn
31  ! -- kluge note: no check whether conversion is possible; assumes it is
32 
33  ! -- Translate and rotate the rectangular cell into local coordinates
34  ! -- with x varying from 0 to dx and y varying from 0 to dy. Choose the
35  ! -- "south-west" vertex to be the local origin so that the rotation
36  ! -- angle is zero if the cell already aligns with the model x and y
37  ! -- coordinates. The "southwest" vertex is found by finding the vertex
38  ! -- formed either by (1) an edge over which x decreases (going
39  ! -- clockwise) followed by an edge over which x does not increase, or
40  ! -- by (2) an edge over which y does not decrease (again going
41  ! -- clockwise) followed by an edge over which y increases. In the end,
42  ! -- ipv1 is the local vertex number (within the cell, taking a value
43  ! -- of 1, 2, 3, or 4) of the southwest vertex, and ipv2, ipv3, and
44  ! -- ipv4 are the local vertex numbers of the remaining three vertices
45  ! -- going clockwise.
46  do ipv = 1, 4
47  ipv4 = ipv
48  ipv1 = ipvnxt(ipv4)
49  x4 = defn%polyvert(1, ipv4)
50  y4 = defn%polyvert(2, ipv4)
51  x1 = defn%polyvert(1, ipv1)
52  y1 = defn%polyvert(2, ipv1)
53  if (x1 .lt. x4) then
54  ipv2 = ipvnxt(ipv1)
55  x2 = defn%polyvert(1, ipv2)
56  if (x2 .le. x1) then
57  y2 = defn%polyvert(2, ipv2)
58  exit
59  end if
60  else if (y1 .ge. y4) then
61  ipv2 = ipvnxt(ipv1)
62  y2 = defn%polyvert(2, ipv2)
63  if (y2 .gt. y1) then
64  x2 = defn%polyvert(1, ipv2)
65  exit
66  end if
67  end if
68  end do
69  ipv3 = ipvnxt(ipv2)
70 
71  ! -- Compute upper bounds on the local coordinates (the rectangular
72  ! -- dimensions of the cell) and the sine and cosine of the rotation
73  ! -- angle, and store local origin information
74  xorigin = x1
75  yorigin = y1
76  zorigin = defn%bot
77  dx2 = x2 - xorigin
78  dy2 = y2 - yorigin
79  dx4 = x4 - xorigin
80  dy4 = y4 - yorigin
81  dx = dsqrt(dx4 * dx4 + dy4 * dy4)
82  dy = dsqrt(dx2 * dx2 + dy2 * dy2)
83  dz = defn%top - zorigin ! todo: need to account for partial saturation
84  sinrot = dy4 / dx
85  cosrot = dx4 / dx
86  rect%defn = poly%defn
87  rect%dx = dx
88  rect%dy = dy
89  rect%dz = dz
90  rect%sinrot = sinrot
91  rect%cosrot = cosrot
92  rect%xOrigin = xorigin
93  rect%yOrigin = yorigin
94  rect%zOrigin = zorigin
95  rect%ipvOrigin = ipv1
96 
97  ! -- Compute (unscaled) cell edge velocities from face flows
98  areax = dx * dz
99  areay = dy * dz
100  areaz = dx * dy
101  factor = done / (defn%retfactor * defn%porosity)
102  term = factor / areax
103  rect%vx1 = defn%faceflow(ipv1) * term
104  rect%vx2 = -defn%faceflow(ipv3) * term
105  term = factor / areay
106  rect%vy1 = defn%faceflow(ipv4) * term
107  rect%vy2 = -defn%faceflow(ipv2) * term
108  term = factor / areaz
109  rect%vz1 = defn%faceflow(6) * term
110  rect%vz2 = -defn%faceflow(7) * term
subroutine, public create_cell_rect(cell)
Create a new rectangular cell.
Definition: CellRect.f90:40
This module contains simulation constants.
Definition: Constants.f90:9
real(dp), parameter done
real constant 1
Definition: Constants.f90:75
Base grid cell definition.
Definition: CellDefn.f90:10
Here is the call graph for this function:
Here is the caller graph for this function: