MODFLOW 6  version 6.7.0.dev0
USGS Modular Hydrologic Model
Data Types | Functions/Subroutines
methodcellternarymodule Module Reference

Data Types

type  methodcellternarytype
 

Functions/Subroutines

subroutine, public create_method_cell_ternary (method)
 Create a tracking method. More...
 
subroutine destroy_mct (this)
 Destroy the tracking method. More...
 
subroutine load_mct (this, particle, next_level, submethod)
 Load subcell into tracking method. More...
 
subroutine pass_mct (this, particle)
 Pass particle to next subcell if there is one, or to the cell face. More...
 
subroutine apply_mct (this, particle, tmax)
 Apply the ternary method to a polygonal cell. More...
 
subroutine load_subcell (this, particle, subcell)
 Loads a triangular subcell from the polygonal cell. More...
 
subroutine vertvelo (this)
 Calculate vertex velocities. More...
 
subroutine calc_thru_hcsum (this, vm0ival, divcell, le, li, lm, areasub, areacell, unintx, uninty, unextx, unexty, unintxnext, unintynext, kappax, kappay, vm0x, vm0y, vm1x, vm1y, hcsum)
 

Function/Subroutine Documentation

◆ apply_mct()

subroutine methodcellternarymodule::apply_mct ( class(methodcellternarytype), intent(inout)  this,
type(particletype), intent(inout), pointer  particle,
real(dp), intent(in)  tmax 
)
private

Definition at line 150 of file MethodCellTernary.f90.

151  use constantsmodule, only: dzero, done, dhalf
152  ! dummy
153  class(MethodCellTernaryType), intent(inout) :: this
154  type(ParticleType), pointer, intent(inout) :: particle
155  real(DP), intent(in) :: tmax
156  ! local
157  integer(I4B) :: i
158 
159  ! (Re)allocate type-bound arrays
160  select type (cell => this%cell)
161  type is (cellpolytype)
162  ! Check termination/reporting conditions
163  call this%check(particle, this%cell%defn)
164  if (.not. particle%advancing) return
165 
166  ! Number of vertices
167  this%nverts = cell%defn%npolyverts
168 
169  ! (Re)allocate type-bound arrays
170  if (allocated(this%xvert)) then
171  deallocate (this%xvert)
172  deallocate (this%yvert)
173  deallocate (this%vne)
174  deallocate (this%vv0x)
175  deallocate (this%vv0y)
176  deallocate (this%vv1x)
177  deallocate (this%vv1y)
178  deallocate (this%iprev)
179  deallocate (this%xvertnext)
180  deallocate (this%yvertnext)
181  end if
182 
183  allocate (this%xvert(this%nverts))
184  allocate (this%yvert(this%nverts))
185  allocate (this%vne(this%nverts))
186  allocate (this%vv0x(this%nverts))
187  allocate (this%vv0y(this%nverts))
188  allocate (this%vv1x(this%nverts))
189  allocate (this%vv1y(this%nverts))
190  allocate (this%iprev(this%nverts))
191  allocate (this%xvertnext(this%nverts))
192  allocate (this%yvertnext(this%nverts))
193 
194  ! Cell vertices
195  do i = 1, this%nverts
196  this%xvert(i) = cell%defn%polyvert(1, i)
197  this%yvert(i) = cell%defn%polyvert(2, i)
198  end do
199 
200  ! Top, bottom, and thickness
201  this%ztop = cell%defn%top
202  this%zbot = cell%defn%bot
203  this%dz = this%ztop - this%zbot
204 
205  ! Shifted arrays
206  do i = 1, this%nverts
207  this%iprev(i) = i
208  end do
209  this%iprev = cshift(this%iprev, -1)
210  this%xvertnext = cshift(this%xvert, 1)
211  this%yvertnext = cshift(this%yvert, 1)
212  end select
213 
214  ! Calculate vertex velocities.
215  call this%vertvelo()
216 
217  ! Track the particle across the cell.
218  call this%track(particle, 2, tmax)
219 
constantsmodule
This module contains simulation constants.
Definition: Constants.f90:9
constantsmodule::dhalf
real(dp), parameter dhalf
real constant 1/2
Definition: Constants.f90:68
constantsmodule::dzero
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
constantsmodule::done
real(dp), parameter done
real constant 1
Definition: Constants.f90:76

◆ calc_thru_hcsum()

subroutine methodcellternarymodule::calc_thru_hcsum ( class(methodcellternarytype), intent(inout)  this,
real(dp)  vm0ival,
real(dp)  divcell,
real(dp), dimension(:)  le,
real(dp), dimension(:)  li,
real(dp), dimension(:)  lm,
real(dp), dimension(:)  areasub,
real(dp)  areacell,
real(dp), dimension(:)  unintx,
real(dp), dimension(:)  uninty,
real(dp), dimension(:)  unextx,
real(dp), dimension(:)  unexty,
real(dp), dimension(:)  unintxnext,
real(dp), dimension(:)  unintynext,
real(dp), dimension(:)  kappax,
real(dp), dimension(:)  kappay,
real(dp), dimension(:)  vm0x,
real(dp), dimension(:)  vm0y,
real(dp), dimension(:)  vm1x,
real(dp), dimension(:)  vm1y,
real(dp)  hcsum 
)

Definition at line 519 of file MethodCellTernary.f90.

524  ! dummy
525  class(MethodCellTernaryType), intent(inout) :: this
526  real(DP) :: vm0ival
527  real(DP) :: divcell
528  real(DP) :: hcsum
529  real(DP), dimension(:) :: le
530  real(DP), dimension(:) :: li
531  real(DP), dimension(:) :: lm
532  real(DP), dimension(:) :: areasub
533  real(DP) :: areacell
534  real(DP), dimension(:) :: unintx
535  real(DP), dimension(:) :: uninty
536  real(DP), dimension(:) :: unextx
537  real(DP), dimension(:) :: unexty
538  real(DP), dimension(:) :: unintxnext
539  real(DP), dimension(:) :: unintynext
540  real(DP), dimension(:) :: kappax
541  real(DP), dimension(:) :: kappay
542  real(DP), dimension(:) :: vm0x
543  real(DP), dimension(:) :: vm0y
544  real(DP), dimension(:) :: vm1x
545  real(DP), dimension(:) :: vm1y
546  ! local
547  real(DP), allocatable, dimension(:) :: vm0i
548  real(DP), allocatable, dimension(:) :: vm0e
549  real(DP), allocatable, dimension(:) :: vm1i
550  real(DP), allocatable, dimension(:) :: vm1e
551  real(DP), allocatable, dimension(:) :: uprod
552  real(DP), allocatable, dimension(:) :: det
553  real(DP), allocatable, dimension(:) :: wt
554  real(DP), allocatable, dimension(:) :: bi0x
555  real(DP), allocatable, dimension(:) :: be0x
556  real(DP), allocatable, dimension(:) :: bi0y
557  real(DP), allocatable, dimension(:) :: be0y
558  real(DP), allocatable, dimension(:) :: bi1x
559  real(DP), allocatable, dimension(:) :: be1x
560  real(DP), allocatable, dimension(:) :: bi1y
561  real(DP), allocatable, dimension(:) :: be1y
562  real(DP), allocatable, dimension(:) :: be01x
563  real(DP), allocatable, dimension(:) :: be01y
564  real(DP) :: emxx
565  real(DP) :: emxy
566  real(DP) :: emyx
567  real(DP) :: emyy
568  real(DP) :: rx
569  real(DP) :: ry
570  real(DP) :: emdet
571  integer(I4B) :: i
572  integer(I4B) :: ip
573 
574  ! Allocate local arrays
575  allocate (vm0i(this%nverts))
576  allocate (vm0e(this%nverts))
577  allocate (vm1i(this%nverts))
578  allocate (vm1e(this%nverts))
579  allocate (uprod(this%nverts))
580  allocate (det(this%nverts))
581  allocate (wt(this%nverts))
582  allocate (bi0x(this%nverts))
583  allocate (be0x(this%nverts))
584  allocate (bi0y(this%nverts))
585  allocate (be0y(this%nverts))
586  allocate (bi1x(this%nverts))
587  allocate (be1x(this%nverts))
588  allocate (bi1y(this%nverts))
589  allocate (be1y(this%nverts))
590  allocate (be01x(this%nverts))
591  allocate (be01y(this%nverts))
592 
593  ! Set vm0i(1)
594  vm0i(1) = vm0ival
595 
596  ! Get remaining vm0i's sequentially using divergence conditions
597  do i = 2, this%nverts
598  ip = this%iprev(i)
599  vm0i(i) = (li(ip) * vm0i(ip) - le(ip) * this%vne(ip) &
600  + areasub(ip) * divcell) / li(i)
601  end do
602 
603  ! Get vm1i's from vm0i's using continuity conditions
604  vm1i = cshift(vm0i, 1)
605 
606  ! Get centroid velocity by setting up and solving 2x2 linear system
607  uprod = unintx * unextx + uninty * unexty
608  det = done - uprod * uprod
609  bi0x = (unintx - unextx * uprod) / det
610  be0x = (unextx - unintx * uprod) / det
611  bi0y = (uninty - unexty * uprod) / det
612  be0y = (unexty - uninty * uprod) / det
613  uprod = unintxnext * unextx + unintynext * unexty
614  det = done - uprod * uprod
615  bi1x = (unintxnext - unextx * uprod) / det
616  be1x = (unextx - unintxnext * uprod) / det
617  bi1y = (unintynext - unexty * uprod) / det
618  be1y = (unexty - unintynext * uprod) / det
619  be01x = 5.d-1 * (be0x + be1x)
620  be01y = 5.d-1 * (be0y + be1y)
621  wt = areasub / areacell
622  emxx = dtwo - sum(wt * be01x * unextx)
623  emxy = -sum(wt * be01x * unexty)
624  emyx = -sum(wt * be01y * unextx)
625  emyy = dtwo - sum(wt * be01y * unexty)
626  rx = sum(wt * (bi0x * vm0i + bi1x * vm1i + be01x * this%vne))
627  ry = sum(wt * (bi0y * vm0i + bi1y * vm1i + be01y * this%vne))
628  emdet = emxx * emyy - emxy * emyx
629  this%vctrx = (emyy * rx - emxy * ry) / emdet
630  this%vctry = (emxx * ry - emyx * rx) / emdet
631 
632  ! Get vm0e's using "known" conditions
633  vm0e = 5.d-1 * (this%vne + unextx * this%vctrx + unexty * this%vctry)
634 
635  ! Get vm1e's from uniformity along exterior edges
636  vm1e = vm0e
637 
638  ! Transform vm0 and vm1 to (x, y) coordinates
639  vm0x = bi0x * vm0i + be0x * vm0e
640  vm0y = bi0y * vm0i + be0y * vm0e
641  vm1x = bi1x * vm1i + be1x * vm0e
642  vm1y = bi1y * vm1i + be1y * vm0e
643 
644  ! Calculate head-cycle summation (which is proportional to
645  ! the curl of the head gradient)
646  hcsum = sum(lm * (kappax * (vm0x + vm1x) + kappay * (vm0y + vm1y)))
647 
648  ! Deallocate local arrays
649  deallocate (vm0i)
650  deallocate (vm0e)
651  deallocate (vm1i)
652  deallocate (vm1e)
653  deallocate (uprod)
654  deallocate (det)
655  deallocate (wt)
656  deallocate (bi0x)
657  deallocate (be0x)
658  deallocate (bi0y)
659  deallocate (be0y)
660  deallocate (bi1x)
661  deallocate (be1x)
662  deallocate (bi1y)
663  deallocate (be1y)
664  deallocate (be01x)
665  deallocate (be01y)
666 

◆ create_method_cell_ternary()

subroutine, public methodcellternarymodule::create_method_cell_ternary ( type(methodcellternarytype), pointer  method)

Definition at line 55 of file MethodCellTernary.f90.

56  ! dummy
57  type(MethodCellTernaryType), pointer :: method
58  ! local
59  type(CellPolyType), pointer :: cell
60  type(SubcellTriType), pointer :: subcell
61 
62  allocate (method)
63  call create_cell_poly(cell)
64  method%cell => cell
65  method%name => method%cell%type
66  method%delegates = .true.
67  call create_subcell_tri(subcell)
68  method%subcell => subcell
Here is the call graph for this function:
Here is the caller graph for this function:

◆ destroy_mct()

subroutine methodcellternarymodule::destroy_mct ( class(methodcellternarytype), intent(inout)  this)
private

Definition at line 72 of file MethodCellTernary.f90.

73  class(MethodCellTernaryType), intent(inout) :: this
74  deallocate (this%name)

◆ load_mct()

subroutine methodcellternarymodule::load_mct ( class(methodcellternarytype), intent(inout)  this,
type(particletype), intent(inout), pointer  particle,
integer(i4b), intent(in)  next_level,
class(methodtype), intent(inout), pointer  submethod 
)
private

Definition at line 78 of file MethodCellTernary.f90.

79  ! dummy
80  class(MethodCellTernaryType), intent(inout) :: this
81  type(ParticleType), pointer, intent(inout) :: particle
82  integer(I4B), intent(in) :: next_level
83  class(MethodType), pointer, intent(inout) :: submethod
84 
85  select type (subcell => this%subcell)
86  type is (subcelltritype)
87  call this%load_subcell(particle, subcell)
88  end select
89  call method_subcell_tern%init( &
90  fmi=this%fmi, &
91  cell=this%cell, &
92  subcell=this%subcell, &
93  trackctl=this%trackctl, &
94  tracktimes=this%tracktimes)
95  submethod => method_subcell_tern

◆ load_subcell()

subroutine methodcellternarymodule::load_subcell ( class(methodcellternarytype), intent(inout)  this,
type(particletype), intent(inout), pointer  particle,
class(subcelltritype), intent(inout)  subcell 
)

Definition at line 223 of file MethodCellTernary.f90.

224  ! dummy
225  class(MethodCellTernaryType), intent(inout) :: this
226  type(ParticleType), pointer, intent(inout) :: particle
227  class(SubcellTriType), intent(inout) :: subcell
228  ! local
229  integer(I4B) :: ic
230  integer(I4B) :: isc
231  integer(I4B) :: iv0
232  real(DP) :: x0
233  real(DP) :: y0
234  real(DP) :: x1
235  real(DP) :: y1
236  real(DP) :: x2
237  real(DP) :: y2
238  real(DP) :: x1rel
239  real(DP) :: y1rel
240  real(DP) :: x2rel
241  real(DP) :: y2rel
242  real(DP) :: xi
243  real(DP) :: yi
244  real(DP) :: di2
245  real(DP) :: d02
246  real(DP) :: d12
247  real(DP) :: di1
248  real(DP) :: d01
249  real(DP) :: alphai
250  real(DP) :: betai
251 
252  select type (cell => this%cell)
253  type is (cellpolytype)
254  ic = cell%defn%icell
255  subcell%icell = ic
256  isc = particle%idomain(3)
257  if (isc .le. 0) then
258  xi = particle%x
259  yi = particle%y
260  do iv0 = 1, this%nverts
261  x0 = this%xvert(iv0)
262  y0 = this%yvert(iv0)
263  x1 = this%xvertnext(iv0)
264  y1 = this%yvertnext(iv0)
265  x2 = this%xctr
266  y2 = this%yctr
267  x1rel = x1 - x0
268  y1rel = y1 - y0
269  x2rel = x2 - x0
270  y2rel = y2 - y0
271  di2 = xi * y2rel - yi * x2rel
272  d02 = x0 * y2rel - y0 * x2rel
273  d12 = x1rel * y2rel - y1rel * x2rel
274  di1 = xi * y1rel - yi * x1rel
275  d01 = x0 * y1rel - y0 * x1rel
276  alphai = (di2 - d02) / d12
277  betai = -(di1 - d01) / d12
278  ! assumes particle is in the cell, so no check needed for beta
279  if ((alphai .ge. dzero) .and. &
280  (alphai + betai .le. done)) then
281  isc = iv0
282  exit
283  end if
284  end do
285  if (isc .le. 0) then
286  print *, "error -- initial triangle not found in cell ", ic, &
287  " for particle at ", particle%x, particle%y, particle%z
288 
289  call pstop(1)
290  else
291  particle%idomain(3) = isc
292  end if
293  end if
294  subcell%isubcell = isc
295 
296  ! Set coordinates and velocities at vertices of triangular subcell
297  iv0 = isc
298  subcell%x0 = this%xvert(iv0)
299  subcell%y0 = this%yvert(iv0)
300  subcell%x1 = this%xvertnext(iv0)
301  subcell%y1 = this%yvertnext(iv0)
302  subcell%x2 = this%xctr
303  subcell%y2 = this%yctr
304  subcell%v0x = this%vv0x(iv0)
305  subcell%v0y = this%vv0y(iv0)
306  subcell%v1x = this%vv1x(iv0) ! the indices here actually refer to subcells, not vertices
307  subcell%v1y = this%vv1y(iv0)
308  subcell%v2x = this%vctrx
309  subcell%v2y = this%vctry
310  subcell%ztop = this%ztop
311  subcell%zbot = this%zbot
312  subcell%dz = this%dz
313  subcell%vzbot = this%vzbot
314  subcell%vztop = this%vztop
315  end select
Here is the call graph for this function:

◆ pass_mct()

subroutine methodcellternarymodule::pass_mct ( class(methodcellternarytype), intent(inout)  this,
type(particletype), intent(inout), pointer  particle 
)
private

Definition at line 99 of file MethodCellTernary.f90.

100  ! dummy
101  class(MethodCellTernaryType), intent(inout) :: this
102  type(ParticleType), pointer, intent(inout) :: particle
103  ! local
104  integer(I4B) :: isc
105  integer(I4B) :: exitFace
106  integer(I4B) :: inface
107 
108  exitface = particle%iboundary(3)
109  isc = particle%idomain(3)
110 
111  select case (exitface)
112  case (0)
113  ! Subcell interior (cell interior)
114  inface = -1
115  case (1)
116  ! Subcell face 1 (cell face)
117  inface = isc
118  if (inface .eq. 0) inface = this%nverts
119  case (2)
120  ! Subcell face --> next subcell in "cycle" (cell interior)
121  isc = isc + 1
122  if (isc .gt. this%nverts) isc = 1
123  particle%idomain(3) = isc
124  particle%iboundary(3) = 3
125  inface = 0
126  case (3)
127  ! Subcell face --> preceding subcell in "cycle" (cell interior)
128  isc = isc - 1
129  if (isc .lt. 1) isc = this%nverts
130  particle%idomain(3) = isc
131  particle%iboundary(3) = 2
132  inface = 0
133  case (4)
134  ! Subcell bottom (cell bottom)
135  inface = this%nverts + 2
136  case (5)
137  ! Subcell top (cell top)
138  inface = this%nverts + 3
139  end select
140  if (inface .eq. -1) then
141  particle%iboundary(2) = 0
142  else if (inface .eq. 0) then
143  particle%iboundary(2) = 0
144  else
145  particle%iboundary(2) = inface
146  end if

◆ vertvelo()

subroutine methodcellternarymodule::vertvelo ( class(methodcellternarytype), intent(inout)  this)
private

Definition at line 319 of file MethodCellTernary.f90.

320  use constantsmodule, only: dzero, done, dhalf
321  ! dummy
322  class(MethodCellTernaryType), intent(inout) :: this
323  ! local
324  real(DP) :: term
325  integer(I4B) :: i
326  real(DP) :: perturb
327  real(DP), allocatable, dimension(:) :: xvals
328  real(DP), allocatable, dimension(:) :: yvals
329  real(DP) :: sixa
330  real(DP) :: vm0i0
331  real(DP) :: vm0ival
332  real(DP) :: hcsum0
333  real(DP) :: hcsum
334  real(DP) :: jac
335  real(DP), allocatable, dimension(:) :: wk1
336  real(DP), allocatable, dimension(:) :: wk2
337  real(DP), allocatable, dimension(:) :: unextxnext
338  real(DP), allocatable, dimension(:) :: unextynext
339  real(DP), allocatable, dimension(:) :: le
340  real(DP), allocatable, dimension(:) :: unextx
341  real(DP), allocatable, dimension(:) :: unexty
342  real(DP) :: areacell
343  real(DP), allocatable, dimension(:) :: areasub
344  real(DP) :: divcell
345  real(DP), allocatable, dimension(:) :: li
346  real(DP), allocatable, dimension(:) :: unintx
347  real(DP), allocatable, dimension(:) :: uninty
348  real(DP), allocatable, dimension(:) :: xmid
349  real(DP), allocatable, dimension(:) :: ymid
350  real(DP), allocatable, dimension(:) :: lm
351  real(DP), allocatable, dimension(:) :: umx
352  real(DP), allocatable, dimension(:) :: umy
353  real(DP), allocatable, dimension(:) :: kappax
354  real(DP), allocatable, dimension(:) :: kappay
355  real(DP), allocatable, dimension(:) :: vm0x
356  real(DP), allocatable, dimension(:) :: vm0y
357  real(DP), allocatable, dimension(:) :: vm1x
358  real(DP), allocatable, dimension(:) :: vm1y
359 
360  select type (cell => this%cell)
361  type is (cellpolytype)
362 
363  ! Allocate local arrays
364  allocate (le(this%nverts)) ! lengths of exterior (cell) edges
365  allocate (unextx(this%nverts)) ! x components of unit normals to exterior edges
366  allocate (unexty(this%nverts)) ! y components of unit normals to exterior edges
367  allocate (areasub(this%nverts)) ! subcell areas
368  allocate (li(this%nverts)) ! lengths of interior edges ("spokes")
369  allocate (unintx(this%nverts)) ! x components of unit normals to interior edges
370  allocate (uninty(this%nverts)) ! y components of unit normals to interior edges
371  allocate (xmid(this%nverts)) ! x coordinates of midpoints
372  allocate (ymid(this%nverts)) ! y coordinates of midpoints
373  allocate (lm(this%nverts)) ! lengths of midpoint connectors
374  allocate (umx(this%nverts)) ! x components of midpoint-connector (cw) unit vectors
375  allocate (umy(this%nverts)) ! y components of midpoint-connector (cw) unit vectors
376  allocate (kappax(this%nverts)) ! x components of kappa vectors
377  allocate (kappay(this%nverts)) ! y components of kappa vectors
378  allocate (vm0x(this%nverts)) ! x component of vm0
379  allocate (vm0y(this%nverts)) ! y component of vm0
380  allocate (vm1x(this%nverts)) ! x component of vm1
381  allocate (vm1y(this%nverts)) ! y component of vm1
382  allocate (unextxnext(this%nverts)) ! vector of "next" interior unit-normal x coordinates defined for convenience
383  allocate (unextynext(this%nverts)) ! vector of "next" interior unit-normal y coordinates defined for convenience
384  allocate (wk1(this%nverts))
385  allocate (wk2(this%nverts))
386  allocate (xvals(3))
387  allocate (yvals(3))
388 
389  ! Exterior edge unit normals (outward) and lengths
390  wk1 = this%xvertnext - this%xvert
391  wk2 = this%yvertnext - this%yvert
392  le = dsqrt(wk1 * wk1 + wk2 * wk2)
393  unextx = -wk2 / le
394  unexty = wk1 / le
395 
396  ! Cell area
397  areacell = area(this%xvert, this%yvert)
398 
399  ! Cell centroid (in general, this is NOT the average of the vertex coordinates)
400  sixa = areacell * 6.d0
401  wk1 = -(this%xvert * this%yvertnext - this%xvertnext * this%yvert)
402  this%xctr = sum((this%xvert + this%xvertnext) * wk1) / sixa
403  this%yctr = sum((this%yvert + this%yvertnext) * wk1) / sixa
404 
405  ! Subcell areas
406  do i = 1, this%nverts
407  xvals(1) = this%xvert(i)
408  xvals(2) = this%xvertnext(i)
409  xvals(3) = this%xctr
410  yvals(1) = this%yvert(i)
411  yvals(2) = this%yvertnext(i)
412  yvals(3) = this%yctr
413  areasub(i) = area(xvals, yvals)
414  end do
415 
416  ! Cell-edge normal velocities (outward)
417  term = done / (cell%defn%porosity * cell%defn%retfactor * this%dz)
418  do i = 1, this%nverts
419  this%vne(i) = -cell%defn%faceflow(i) * term / le(i)
420  end do
421 
422  ! Cell divergence (2D)
423  divcell = sum(le * this%vne) / areacell
424 
425  ! Interior edge (cw) unit normals and lengths
426  wk1 = this%xvert - this%xctr
427  wk2 = this%yvert - this%yctr
428  li = dsqrt(wk1 * wk1 + wk2 * wk2)
429  unintx = wk2 / li
430  uninty = -wk1 / li
431  ! Shifted arrays for convenience
432  unextxnext = cshift(unintx, 1)
433  unextynext = cshift(uninty, 1)
434 
435  ! Midpoints of interior edges
436  xmid = 5.d-1 * (this%xvert + this%xctr)
437  ymid = 5.d-1 * (this%yvert + this%yctr)
438 
439  ! Unit midpoint-connector (cw) vectors and lengths
440  wk1 = cshift(xmid, 1) - xmid
441  wk2 = cshift(ymid, 1) - ymid
442  lm = dsqrt(wk1 * wk1 + wk2 * wk2)
443  umx = wk1 / lm
444  umy = wk2 / lm
445 
446  ! Kappa vectors (K tensor times unit midpoint-connector vectors) do not
447  ! account for anisotropy, which is consistent with the way internal face
448  ! flow calculations are done in MP7. The isotropic value of K does not
449  ! matter in this case because it cancels out of the calculations, so
450  ! K = 1 is assumed for simplicity.
451  kappax = umx
452  kappay = umy
453 
454  ! Use linearity to find vm0i[0] such that curl of the head gradient
455  ! is zero
456  perturb = 1.d-2
457  ! Calculations at base value
458  vm0i0 = 0.d0
459  call this%calc_thru_hcsum(vm0i0, divcell, le, li, lm, areasub, areacell, &
460  unintx, uninty, unextx, unexty, &
461  unextxnext, unextynext, &
462  kappax, kappay, vm0x, vm0y, vm1x, vm1y, hcsum0)
463  ! Calculations at perturbed value
464  vm0ival = vm0i0 + perturb
465  call this%calc_thru_hcsum(vm0ival, divcell, le, li, lm, areasub, areacell, &
466  unintx, uninty, unextx, unexty, &
467  unextxnext, unextynext, &
468  kappax, kappay, vm0x, vm0y, vm1x, vm1y, hcsum)
469  ! Calculations at root value
470  jac = (hcsum - hcsum0) / perturb
471  vm0ival = vm0i0 - hcsum0 / jac
472  call this%calc_thru_hcsum(vm0ival, divcell, le, li, lm, areasub, areacell, &
473  unintx, uninty, unextx, unexty, &
474  unextxnext, unextynext, &
475  kappax, kappay, vm0x, vm0y, vm1x, vm1y, hcsum)
476 
477  ! Project linearly to get corner (vertex) velocities. Note that velocity
478  ! vv1 is at the next vertex cw from vv0, so vv0(i) and vv1(i) are the
479  ! two vertex velocities used by triangular subcell i.
480  this%vv0x = 2.d0 * vm0x - this%vctrx
481  this%vv0y = 2.d0 * vm0y - this%vctry
482  this%vv1x = 2.d0 * vm1x - this%vctrx
483  this%vv1y = 2.d0 * vm1y - this%vctry
484 
485  ! Set top and bottom velocities
486  term = done / (cell%defn%retfactor * cell%defn%porosity * areacell)
487  this%vzbot = cell%defn%faceflow(this%nverts + 2) * term
488  this%vztop = -cell%defn%faceflow(this%nverts + 3) * term
489 
490  ! Deallocate local arrays
491  deallocate (le)
492  deallocate (unextx)
493  deallocate (unexty)
494  deallocate (areasub)
495  deallocate (li)
496  deallocate (unintx)
497  deallocate (uninty)
498  deallocate (xmid)
499  deallocate (ymid)
500  deallocate (lm)
501  deallocate (umx)
502  deallocate (umy)
503  deallocate (kappax)
504  deallocate (kappay)
505  deallocate (vm0x)
506  deallocate (vm0y)
507  deallocate (vm1x)
508  deallocate (vm1y)
509  deallocate (unextxnext)
510  deallocate (unextynext)
511  deallocate (wk1)
512  deallocate (wk2)
513  deallocate (xvals)
514  deallocate (yvals)
515 
516  end select
Here is the call graph for this function: