MODFLOW 6  version 6.7.0.dev3
USGS Modular Hydrologic Model
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 154 of file MethodCellTernary.f90.

155  use constantsmodule, only: dzero, done, dhalf
156  ! dummy
157  class(MethodCellTernaryType), intent(inout) :: this
158  type(ParticleType), pointer, intent(inout) :: particle
159  real(DP), intent(in) :: tmax
160  ! local
161  integer(I4B) :: i
162  real(DP) :: x, y, z, xO, yO
163  real(DP), allocatable :: xs(:), ys(:)
164 
165  ! (Re)allocate type-bound arrays
166  select type (cell => this%cell)
167  type is (cellpolytype)
168  call this%assess(particle, this%cell%defn, tmax)
169  if (.not. particle%advancing) return
170 
171  ! Number of vertices
172  this%nverts = cell%defn%npolyverts
173 
174  ! (Re)allocate type-bound arrays
175  if (allocated(this%xvert)) then
176  deallocate (this%xvert)
177  deallocate (this%yvert)
178  deallocate (this%vne)
179  deallocate (this%vv0x)
180  deallocate (this%vv0y)
181  deallocate (this%vv1x)
182  deallocate (this%vv1y)
183  deallocate (this%iprev)
184  deallocate (this%xvertnext)
185  deallocate (this%yvertnext)
186  end if
187 
188  allocate (this%xvert(this%nverts))
189  allocate (this%yvert(this%nverts))
190  allocate (this%vne(this%nverts))
191  allocate (this%vv0x(this%nverts))
192  allocate (this%vv0y(this%nverts))
193  allocate (this%vv1x(this%nverts))
194  allocate (this%vv1y(this%nverts))
195  allocate (this%iprev(this%nverts))
196  allocate (this%xvertnext(this%nverts))
197  allocate (this%yvertnext(this%nverts))
198 
199  ! New origin is the corner of the cell's
200  ! bounding box closest to the old origin
201  allocate (xs(this%nverts))
202  allocate (ys(this%nverts))
203  xs = cell%defn%polyvert(1, :)
204  ys = cell%defn%polyvert(2, :)
205  xo = xs(minloc(abs(xs), dim=1))
206  yo = ys(minloc(abs(ys), dim=1))
207  deallocate (xs)
208  deallocate (ys)
209 
210  ! Cell vertices
211  do i = 1, this%nverts
212  x = cell%defn%polyvert(1, i)
213  y = cell%defn%polyvert(2, i)
214  call transform(x, y, dzero, x, y, z, xo, yo)
215  this%xvert(i) = x
216  this%yvert(i) = y
217  end do
218 
219  ! Top, bottom, and thickness
220  this%ztop = cell%defn%top
221  this%zbot = cell%defn%bot
222  this%dz = this%ztop - this%zbot
223 
224  ! Shifted arrays
225  do i = 1, this%nverts
226  this%iprev(i) = i
227  end do
228  this%iprev = cshift(this%iprev, -1)
229  this%xvertnext = cshift(this%xvert, 1)
230  this%yvertnext = cshift(this%yvert, 1)
231 
232  ! Calculate vertex velocities.
233  call this%vertvelo()
234 
235  ! Transform particle coordinates
236  call particle%transform(xo, yo)
237 
238  ! Track the particle across the cell.
239  call this%track(particle, 2, tmax)
240 
241  ! Transform particle coordinates back
242  call particle%transform(xo, yo, invert=.true.)
243  call particle%reset_transform()
244  end select
This module contains simulation constants.
Definition: Constants.f90:9
real(dp), parameter dhalf
real constant 1/2
Definition: Constants.f90:68
real(dp), parameter dzero
real constant zero
Definition: Constants.f90:65
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 531 of file MethodCellTernary.f90.

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

◆ create_method_cell_ternary()

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

Definition at line 57 of file MethodCellTernary.f90.

58  ! dummy
59  type(MethodCellTernaryType), pointer :: method
60  ! local
61  type(CellPolyType), pointer :: cell
62  type(SubcellTriType), pointer :: subcell
63 
64  allocate (method)
65  call create_cell_poly(cell)
66  method%cell => cell
67  method%name => method%cell%type
68  method%delegates = .true.
69  call create_subcell_tri(subcell)
70  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 74 of file MethodCellTernary.f90.

75  class(MethodCellTernaryType), intent(inout) :: this
76  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 80 of file MethodCellTernary.f90.

81  ! dummy
82  class(MethodCellTernaryType), intent(inout) :: this
83  type(ParticleType), pointer, intent(inout) :: particle
84  integer(I4B), intent(in) :: next_level
85  class(MethodType), pointer, intent(inout) :: submethod
86 
87  select type (subcell => this%subcell)
88  type is (subcelltritype)
89  call this%load_subcell(particle, subcell)
90  end select
91 
92  call method_subcell_tern%init( &
93  fmi=this%fmi, &
94  cell=this%cell, &
95  subcell=this%subcell, &
96  events=this%events, &
97  tracktimes=this%tracktimes)
98  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 248 of file MethodCellTernary.f90.

249  ! dummy
250  class(MethodCellTernaryType), intent(inout) :: this
251  type(ParticleType), pointer, intent(inout) :: particle
252  class(SubcellTriType), intent(inout) :: subcell
253  ! local
254  integer(I4B) :: ic, isc, iv0
255  real(DP) :: x0, y0, x1, y1, x2, y2, xi, yi
256  real(DP) :: x1rel, y1rel, x2rel, y2rel
257  real(DP) :: di2, d02, d12, di1, d01
258  real(DP) :: alphai, betai
259 
260  select type (cell => this%cell)
261  type is (cellpolytype)
262  ic = cell%defn%icell
263  subcell%icell = ic
264  isc = particle%itrdomain(level_subfeature)
265  if (isc .le. 0) then
266  xi = particle%x
267  yi = particle%y
268  do iv0 = 1, this%nverts
269  x0 = this%xvert(iv0)
270  y0 = this%yvert(iv0)
271  x1 = this%xvertnext(iv0)
272  y1 = this%yvertnext(iv0)
273  x2 = this%xctr
274  y2 = this%yctr
275  x1rel = x1 - x0
276  y1rel = y1 - y0
277  x2rel = x2 - x0
278  y2rel = y2 - y0
279  di2 = xi * y2rel - yi * x2rel
280  d02 = x0 * y2rel - y0 * x2rel
281  d12 = x1rel * y2rel - y1rel * x2rel
282  di1 = xi * y1rel - yi * x1rel
283  d01 = x0 * y1rel - y0 * x1rel
284  alphai = (di2 - d02) / d12
285  betai = -(di1 - d01) / d12
286  ! assumes particle is in the cell, so no check needed for beta
287  if ((alphai .ge. dzero) .and. &
288  (alphai + betai .le. done)) then
289  isc = iv0
290  exit
291  end if
292  end do
293  if (isc .le. 0) then
294  print *, "error -- initial triangle not found in cell ", ic, &
295  " for particle at ", particle%x, particle%y, particle%z
296 
297  call pstop(1)
298  else
299  particle%itrdomain(level_subfeature) = isc
300  end if
301  end if
302  subcell%isubcell = isc
303 
304  ! Set coordinates and velocities at vertices of triangular subcell
305  iv0 = isc
306  subcell%x0 = this%xvert(iv0)
307  subcell%y0 = this%yvert(iv0)
308  subcell%x1 = this%xvertnext(iv0)
309  subcell%y1 = this%yvertnext(iv0)
310  subcell%x2 = this%xctr
311  subcell%y2 = this%yctr
312  subcell%v0x = this%vv0x(iv0)
313  subcell%v0y = this%vv0y(iv0)
314  subcell%v1x = this%vv1x(iv0) ! the indices here actually refer to subcells, not vertices
315  subcell%v1y = this%vv1y(iv0)
316  subcell%v2x = this%vctrx
317  subcell%v2y = this%vctry
318  subcell%ztop = this%ztop
319  subcell%zbot = this%zbot
320  subcell%dz = this%dz
321  subcell%vzbot = this%vzbot
322  subcell%vztop = this%vztop
323  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 102 of file MethodCellTernary.f90.

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

◆ vertvelo()

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

Definition at line 327 of file MethodCellTernary.f90.

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