MODFLOW 6  version 6.7.0.dev1
USGS Modular Hydrologic Model
MethodCellPollockQuad.f90
Go to the documentation of this file.
2 
3  use kindmodule, only: dp, i4b
4  use errorutilmodule, only: pstop
5  use constantsmodule, only: done, dzero
6  use methodmodule, only: methodtype
10  use celldefnmodule, only: celldefntype
12  use particlemodule, only: particletype
13  implicit none
14 
15  private
17  public :: create_method_cell_quad
18 
20  contains
21  procedure, public :: apply => apply_mcpq
22  procedure, public :: deallocate
23  procedure, public :: load => load_mcpq
24  procedure, public :: load_subcell
25  procedure, public :: pass => pass_mcpq
27 
28 contains
29 
30  !> @brief Create a new Pollock quad-refined cell method
31  subroutine create_method_cell_quad(method)
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
45  end subroutine create_method_cell_quad
46 
47  !> @brief Deallocate the Pollock quad-refined cell method
48  subroutine deallocate (this)
49  class(methodcellpollockquadtype), intent(inout) :: this
50  deallocate (this%name)
51  end subroutine deallocate
52 
53  !> @brief Load subcell into tracking method
54  subroutine load_mcpq(this, particle, next_level, submethod)
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  trackctl=this%trackctl, &
69  tracktimes=this%tracktimes)
70  submethod => method_subcell_plck
71  end subroutine load_mcpq
72 
73  !> @brief Pass particle to next subcell if there is one, or to the cell face
74  subroutine pass_mcpq(this, particle)
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(3)
84  isc = particle%idomain(3)
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%idomain(3) = 4
97  particle%iboundary(3) = 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%idomain(3) = 3
102  particle%iboundary(3) = 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%idomain(3) = 2
126  particle%iboundary(3) = 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%idomain(3) = 1
131  particle%iboundary(3) = 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%idomain(3) = 2
139  particle%iboundary(3) = 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%idomain(3) = 3
152  particle%iboundary(3) = 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%idomain(3) = 1
164  particle%iboundary(3) = 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%idomain(3) = 4
169  particle%iboundary(3) = 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(2) = 0
186  else if (inface .eq. 0) then
187  particle%iboundary(2) = 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(2) = inface
197  end if
198  end select
199  end subroutine pass_mcpq
200 
201  !> @brief Solve the quad-rectangular cell via Pollock's method
202  subroutine apply_mcpq(this, particle, tmax)
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  ! Check termination/reporting conditions
213  call this%check(particle, cell%defn, tmax)
214  if (.not. particle%advancing) return
215 
216  ! Transform model coordinates to local cell coordinates
217  ! (translated/rotated but not scaled relative to model)
218  xorigin = cell%xOrigin
219  yorigin = cell%yOrigin
220  zorigin = cell%zOrigin
221  sinrot = cell%sinrot
222  cosrot = cell%cosrot
223  call particle%transform(xorigin, yorigin, zorigin, &
224  sinrot, cosrot)
225 
226  ! Track the particle over the cell
227  call this%track(particle, 2, tmax)
228 
229  ! Transform cell coordinates back to model coordinates
230  call particle%transform(xorigin, yorigin, zorigin, &
231  sinrot, cosrot, invert=.true.)
232  call particle%reset_transform()
233  end select
234  end subroutine apply_mcpq
235 
236  !> @brief Load the rectangular subcell from the rectangular cell
237  subroutine load_subcell(this, particle, subcell)
238  ! dummy
239  class(methodcellpollockquadtype), intent(inout) :: this
240  type(particletype), pointer, intent(inout) :: particle
241  class(subcellrecttype), intent(inout) :: subcell
242  ! local
243  real(DP) :: dx, dy, dz, areax, areay, areaz
244  real(DP) :: dxprel, dyprel
245  integer(I4B) :: isc, npolyverts, m1, m2
246  real(DP) :: qextl1, qextl2, qintl1, qintl2
247  real(DP) :: factor, term
248 
249  select type (cell => this%cell)
250  type is (cellrectquadtype)
251  factor = done / cell%defn%retfactor
252  factor = factor / cell%defn%porosity
253  npolyverts = cell%defn%npolyverts
254 
255  isc = particle%idomain(3)
256  ! Subcells 1, 2, 3, and 4 are Pollock's subcells A, B, C, and D,
257  ! respectively
258 
259  dx = cell%dx
260  dy = cell%dy
261  ! If not already known, determine subcell number
262  if (isc .le. 0) then
263  dxprel = particle%x / dx
264  dyprel = particle%y / dy
265 
266  if (dyprel .ge. 5d-1) then
267  if (dxprel .le. 5d-1) then
268  isc = 4
269  else
270  isc = 1
271  end if
272  else
273  if (dxprel .le. 5d-1) then
274  isc = 3
275  else
276  isc = 2
277  end if
278  end if
279 
280  subcell%isubcell = isc
281  particle%idomain(3) = isc
282  end if
283  dx = 5d-1 * dx
284  dy = 5d-1 * dy
285  dz = cell%defn%top - &
286  cell%defn%bot
287  areax = dy * dz
288  areay = dx * dz
289  areaz = dx * dy
290  qintl1 = cell%qintl(isc)
291  ! qintl list wraps around, so isc+1=5 is ok
292  qintl2 = cell%qintl(isc + 1)
293  qextl1 = cell%qextl1(isc)
294  qextl2 = cell%qextl2(isc)
295 
296  subcell%dx = dx
297  subcell%dy = dy
298  subcell%dz = dz
299  subcell%sinrot = dzero
300  subcell%cosrot = done
301  subcell%zOrigin = dzero
302  select case (isc)
303  case (1)
304  subcell%xOrigin = dx
305  subcell%yOrigin = dy
306  term = factor / areax
307  subcell%vx1 = qintl1 * term
308  subcell%vx2 = -qextl2 * term
309  term = factor / areay
310  subcell%vy1 = -qintl2 * term
311  subcell%vy2 = -qextl1 * term
312  case (2)
313  subcell%xOrigin = dx
314  subcell%yOrigin = dzero
315  term = factor / areax
316  subcell%vx1 = -qintl2 * term
317  subcell%vx2 = -qextl1 * term
318  term = factor / areay
319  subcell%vy1 = qextl2 * term
320  subcell%vy2 = -qintl1 * term
321  case (3)
322  subcell%xOrigin = dzero
323  subcell%yOrigin = dzero
324  term = factor / areax
325  subcell%vx1 = qextl2 * term
326  subcell%vx2 = -qintl1 * term
327  term = factor / areay
328  subcell%vy1 = qextl1 * term
329  subcell%vy2 = qintl2 * term
330  case (4)
331  subcell%xOrigin = dzero
332  subcell%yOrigin = dy
333  term = factor / areax
334  subcell%vx1 = qextl1 * term
335  subcell%vx2 = qintl2 * term
336  term = factor / areay
337  subcell%vy1 = qintl1 * term
338  subcell%vy2 = -qextl2 * term
339  end select
340  m1 = npolyverts + 2
341  m2 = m1 + 1
342  term = factor / areaz
343  subcell%vz1 = 2.5d-1 * cell%defn%faceflow(m1) * term
344  subcell%vz2 = -2.5d-1 * cell%defn%faceflow(m2) * term
345  end select
346  end subroutine load_subcell
347 
subroutine, public create_cell_rect_quad(cell)
Create a new rectangular-quad cell.
This module contains simulation constants.
Definition: Constants.f90:9
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
real(dp), parameter done
real constant 1
Definition: Constants.f90:76
subroutine pstop(status, message)
Stop the program, optionally specifying an error status code.
Definition: ErrorUtil.f90:24
This module defines variable data types.
Definition: kind.f90:8
procedure subroutine, public create_method_cell_quad(method)
Create a new Pollock quad-refined cell method.
subroutine load_mcpq(this, particle, next_level, submethod)
Load subcell into tracking method.
subroutine pass_mcpq(this, particle)
Pass particle to next subcell if there is one, or to the cell face.
subroutine load_subcell(this, particle, subcell)
Load the rectangular subcell from the rectangular cell.
subroutine apply_mcpq(this, particle, tmax)
Solve the quad-rectangular cell via Pollock's method.
Particle tracking strategies.
Definition: Method.f90:2
Subcell-level tracking methods.
type(methodsubcellpollocktype), pointer, public method_subcell_plck
subroutine, public create_subcell_rect(subcell)
Create a new rectangular subcell.
Definition: SubcellRect.f90:28
Base grid cell definition.
Definition: CellDefn.f90:10
Base type for particle tracking methods.
Definition: Method.f90:31
Particle tracked by the PRT model.
Definition: Particle.f90:78