MODFLOW 6  version 6.7.0.dev2
USGS Modular Hydrologic Model
methodcellpollockquadmodule Module Reference

Data Types

type  methodcellpollockquadtype
 

Functions/Subroutines

procedure subroutine, public create_method_cell_quad (method)
 Create a new Pollock quad-refined cell method. More...
 
subroutine deallocate (this)
 Deallocate the Pollock quad-refined cell method. More...
 
subroutine load_mcpq (this, particle, next_level, submethod)
 Load subcell into tracking method. More...
 
subroutine pass_mcpq (this, particle)
 Pass particle to next subcell if there is one, or to the cell face. More...
 
subroutine apply_mcpq (this, particle, tmax)
 Solve the quad-rectangular cell via Pollock's method. More...
 
subroutine load_subcell (this, particle, subcell)
 Load the rectangular subcell from the rectangular cell. More...
 

Function/Subroutine Documentation

◆ apply_mcpq()

subroutine methodcellpollockquadmodule::apply_mcpq ( class(methodcellpollockquadtype), intent(inout)  this,
type(particletype), intent(inout), pointer  particle,
real(dp), intent(in)  tmax 
)
private

Definition at line 202 of file MethodCellPollockQuad.f90.

203  ! dummy
204  class(MethodCellPollockQuadType), intent(inout) :: this
205  type(ParticleType), pointer, intent(inout) :: particle
206  real(DP), intent(in) :: tmax
207  ! local
208  double precision :: xOrigin, yOrigin, zOrigin, sinrot, cosrot
209 
210  select type (cell => this%cell)
211  type is (cellrectquadtype)
212  call this%assess(particle, this%cell%defn, tmax)
213  if (.not. particle%advancing) return
214 
215  ! Transform model coordinates to local cell coordinates
216  ! (translated/rotated but not scaled relative to model)
217  xorigin = cell%xOrigin
218  yorigin = cell%yOrigin
219  zorigin = cell%zOrigin
220  sinrot = cell%sinrot
221  cosrot = cell%cosrot
222  call particle%transform(xorigin, yorigin, zorigin, &
223  sinrot, cosrot)
224 
225  ! Track the particle over the cell
226  call this%track(particle, 2, tmax)
227 
228  ! Transform cell coordinates back to model coordinates
229  call particle%transform(xorigin, yorigin, zorigin, &
230  sinrot, cosrot, invert=.true.)
231  call particle%reset_transform()
232  end select

◆ create_method_cell_quad()

procedure subroutine, public methodcellpollockquadmodule::create_method_cell_quad ( type(methodcellpollockquadtype), pointer  method)

Definition at line 31 of file MethodCellPollockQuad.f90.

32  ! dummy
33  type(MethodCellPollockQuadType), pointer :: method
34  ! local
35  type(CellRectQuadType), pointer :: cell
36  type(SubcellRectType), pointer :: subcell
37 
38  allocate (method)
39  call create_cell_rect_quad(cell)
40  method%cell => cell
41  method%name => method%cell%type
42  method%delegates = .true.
43  call create_subcell_rect(subcell)
44  method%subcell => subcell
Here is the call graph for this function:
Here is the caller graph for this function:

◆ deallocate()

subroutine methodcellpollockquadmodule::deallocate ( class(methodcellpollockquadtype), intent(inout)  this)
private

Definition at line 48 of file MethodCellPollockQuad.f90.

49  class(MethodCellPollockQuadType), intent(inout) :: this
50  deallocate (this%name)

◆ load_mcpq()

subroutine methodcellpollockquadmodule::load_mcpq ( class(methodcellpollockquadtype), intent(inout)  this,
type(particletype), intent(inout), pointer  particle,
integer, intent(in)  next_level,
class(methodtype), intent(inout), pointer  submethod 
)
private

Definition at line 54 of file MethodCellPollockQuad.f90.

55  class(MethodCellPollockQuadType), intent(inout) :: this
56  type(ParticleType), pointer, intent(inout) :: particle
57  integer, intent(in) :: next_level
58  class(MethodType), pointer, intent(inout) :: submethod
59 
60  select type (subcell => this%subcell)
61  type is (subcellrecttype)
62  call this%load_subcell(particle, subcell)
63  end select
64  call method_subcell_plck%init( &
65  fmi=this%fmi, &
66  cell=this%cell, &
67  subcell=this%subcell, &
68  events=this%events, &
69  tracktimes=this%tracktimes)
70  submethod => method_subcell_plck

◆ load_subcell()

subroutine methodcellpollockquadmodule::load_subcell ( class(methodcellpollockquadtype), intent(inout)  this,
type(particletype), intent(inout), pointer  particle,
class(subcellrecttype), intent(inout)  subcell 
)
private

Definition at line 236 of file MethodCellPollockQuad.f90.

237  ! dummy
238  class(MethodCellPollockQuadType), intent(inout) :: this
239  type(ParticleType), pointer, intent(inout) :: particle
240  class(SubcellRectType), intent(inout) :: subcell
241  ! local
242  real(DP) :: dx, dy, dz, areax, areay, areaz
243  real(DP) :: dxprel, dyprel
244  integer(I4B) :: isc, npolyverts, m1, m2
245  real(DP) :: qextl1, qextl2, qintl1, qintl2
246  real(DP) :: factor, term
247 
248  select type (cell => this%cell)
249  type is (cellrectquadtype)
250  factor = done / cell%defn%retfactor
251  factor = factor / cell%defn%porosity
252  npolyverts = cell%defn%npolyverts
253 
254  isc = particle%itrdomain(level_subfeature)
255  ! Subcells 1, 2, 3, and 4 are Pollock's subcells A, B, C, and D,
256  ! respectively
257 
258  dx = cell%dx
259  dy = cell%dy
260  ! If not already known, determine subcell number
261  if (isc .le. 0) then
262  dxprel = particle%x / dx
263  dyprel = particle%y / dy
264 
265  if (dyprel .ge. 5d-1) then
266  if (dxprel .le. 5d-1) then
267  isc = 4
268  else
269  isc = 1
270  end if
271  else
272  if (dxprel .le. 5d-1) then
273  isc = 3
274  else
275  isc = 2
276  end if
277  end if
278 
279  subcell%isubcell = isc
280  particle%itrdomain(level_subfeature) = isc
281  end if
282  dx = 5d-1 * dx
283  dy = 5d-1 * dy
284  dz = cell%defn%top - &
285  cell%defn%bot
286  areax = dy * dz
287  areay = dx * dz
288  areaz = dx * dy
289  qintl1 = cell%qintl(isc)
290  ! qintl list wraps around, so isc+1=5 is ok
291  qintl2 = cell%qintl(isc + 1)
292  qextl1 = cell%qextl1(isc)
293  qextl2 = cell%qextl2(isc)
294 
295  subcell%dx = dx
296  subcell%dy = dy
297  subcell%dz = dz
298  subcell%sinrot = dzero
299  subcell%cosrot = done
300  subcell%zOrigin = dzero
301  select case (isc)
302  case (1)
303  subcell%xOrigin = dx
304  subcell%yOrigin = dy
305  term = factor / areax
306  subcell%vx1 = qintl1 * term
307  subcell%vx2 = -qextl2 * term
308  term = factor / areay
309  subcell%vy1 = -qintl2 * term
310  subcell%vy2 = -qextl1 * term
311  case (2)
312  subcell%xOrigin = dx
313  subcell%yOrigin = dzero
314  term = factor / areax
315  subcell%vx1 = -qintl2 * term
316  subcell%vx2 = -qextl1 * term
317  term = factor / areay
318  subcell%vy1 = qextl2 * term
319  subcell%vy2 = -qintl1 * term
320  case (3)
321  subcell%xOrigin = dzero
322  subcell%yOrigin = dzero
323  term = factor / areax
324  subcell%vx1 = qextl2 * term
325  subcell%vx2 = -qintl1 * term
326  term = factor / areay
327  subcell%vy1 = qextl1 * term
328  subcell%vy2 = qintl2 * term
329  case (4)
330  subcell%xOrigin = dzero
331  subcell%yOrigin = dy
332  term = factor / areax
333  subcell%vx1 = qextl1 * term
334  subcell%vx2 = qintl2 * term
335  term = factor / areay
336  subcell%vy1 = qintl1 * term
337  subcell%vy2 = -qextl2 * term
338  end select
339  m1 = npolyverts + 2
340  m2 = m1 + 1
341  term = factor / areaz
342  subcell%vz1 = 2.5d-1 * cell%defn%faceflow(m1) * term
343  subcell%vz2 = -2.5d-1 * cell%defn%faceflow(m2) * term
344  end select

◆ pass_mcpq()

subroutine methodcellpollockquadmodule::pass_mcpq ( class(methodcellpollockquadtype), intent(inout)  this,
type(particletype), intent(inout), pointer  particle 
)
private

Definition at line 74 of file MethodCellPollockQuad.f90.

75  ! dummy
76  class(MethodCellPollockQuadType), intent(inout) :: this
77  type(ParticleType), pointer, intent(inout) :: particle
78  ! local
79  integer(I4B) :: isc, exitFace, npolyverts, inface, infaceoff
80 
81  select type (cell => this%cell)
82  type is (cellrectquadtype)
83  exitface = particle%iboundary(level_subfeature)
84  isc = particle%itrdomain(level_subfeature)
85  npolyverts = cell%defn%npolyverts
86 
87  ! exitFace uses MODPATH 7 iface convention here
88  select case (exitface)
89  case (0)
90  ! Subcell interior (cell interior)
91  inface = -1
92  case (1)
93  select case (isc)
94  case (1)
95  ! W face, subcell 1 --> E face, subcell 4 (cell interior)
96  particle%itrdomain(level_subfeature) = 4
97  particle%iboundary(level_subfeature) = 2
98  inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
99  case (2)
100  ! W face, subcell 2 --> E face, subcell 3 (cell interior)
101  particle%itrdomain(level_subfeature) = 3
102  particle%iboundary(level_subfeature) = 2
103  inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
104  case (3)
105  ! W face, subcell 3 (cell face)
106  inface = 1 ! want Domain(2) = -Domain(2); Boundary(2) = inface
107  infaceoff = 0
108  case (4)
109  ! W face, subcell 4 (cell face)
110  inface = 2 ! want Domain(2) = -Domain(2); Boundary(2) = inface
111  infaceoff = -1
112  end select
113  case (2)
114  select case (isc)
115  case (1)
116  ! E face, subcell 1 (cell face)
117  inface = 3 ! want Domain(2) = -Domain(2); Boundary(2) = inface
118  infaceoff = 0
119  case (2)
120  ! E face, subcell 2 (cell face)
121  inface = 4 ! want Domain(2) = -Domain(2); Boundary(2) = inface
122  infaceoff = -1
123  case (3)
124  ! E face, subcell 3 --> W face, subcell 2 (cell interior)
125  particle%itrdomain(level_subfeature) = 2
126  particle%iboundary(level_subfeature) = 1
127  inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
128  case (4)
129  ! E face, subcell 4 --> W face subcell 1 (cell interior)
130  particle%itrdomain(level_subfeature) = 1
131  particle%iboundary(level_subfeature) = 1
132  inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
133  end select
134  case (3)
135  select case (isc)
136  case (1)
137  ! S face, subcell 1 --> N face, subcell 2 (cell interior)
138  particle%itrdomain(level_subfeature) = 2
139  particle%iboundary(level_subfeature) = 4
140  inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
141  case (2)
142  ! S face, subcell 2 (cell face)
143  inface = 4 ! want Domain(2) = -Domain(2); Boundary(2) = inface
144  infaceoff = 0
145  case (3)
146  ! S face, subcell 3 (cell face)
147  inface = 1 ! want Domain(2) = -Domain(2); Boundary(2) = inface
148  infaceoff = -1
149  case (4)
150  ! S face, subcell 4 --> N face, subcell 3 (cell interior)
151  particle%itrdomain(level_subfeature) = 3
152  particle%iboundary(level_subfeature) = 4
153  inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
154  end select
155  case (4)
156  select case (isc)
157  case (1)
158  ! N face, subcell 1 (cell face)
159  inface = 3 ! want Domain(2) = -Domain(2); Boundary(2) = inface
160  infaceoff = -1
161  case (2)
162  ! N face, subcell 2 --> S face, subcell 1 (cell interior)
163  particle%itrdomain(level_subfeature) = 1
164  particle%iboundary(level_subfeature) = 3
165  inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
166  case (3)
167  ! N face, subcell 3 --> S face, subcell 4 (cell interior)
168  particle%itrdomain(level_subfeature) = 4
169  particle%iboundary(level_subfeature) = 3
170  inface = 0 ! want Domain(2) unchanged; Boundary(2) = 0
171  case (4)
172  ! N face, subcell 4 (cell face)
173  inface = 2 ! want Domain(2) = -Domain(2); Boundary(2) = inface
174  infaceoff = 0
175  end select
176  case (5)
177  ! Subcell bottom (cell bottom)
178  inface = npolyverts + 2 ! want Domain(2) = -Domain(2); Boundary(2) = inface
179  case (6)
180  ! Subcell top (cell top)
181  inface = npolyverts + 3 ! want Domain(2) = -Domain(2); Boundary(2) = inface
182  end select
183 
184  if (inface .eq. -1) then
185  particle%iboundary(level_feature) = 0
186  else if (inface .eq. 0) then
187  particle%iboundary(level_feature) = 0
188  else
189  if ((inface .ge. 1) .and. (inface .le. 4)) then
190  ! Account for local cell rotation
191  inface = inface + cell%irvOrigin - 1
192  if (inface .gt. 4) inface = inface - 4
193  inface = cell%irectvert(inface) + infaceoff
194  if (inface .lt. 1) inface = inface + npolyverts
195  end if
196  particle%iboundary(level_feature) = inface
197  end if
198  end select