330 class(MethodCellTernaryType),
intent(inout) :: this
335 real(DP),
allocatable,
dimension(:) :: xvals
336 real(DP),
allocatable,
dimension(:) :: yvals
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
351 real(DP),
allocatable,
dimension(:) :: areasub
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
368 integer(I4B) :: nvert
370 select type (cell => this%cell)
371 type is (cellpolytype)
374 allocate (le(this%nverts))
375 allocate (unextx(this%nverts))
376 allocate (unexty(this%nverts))
377 allocate (areasub(this%nverts))
378 allocate (li(this%nverts))
379 allocate (unintx(this%nverts))
380 allocate (uninty(this%nverts))
381 allocate (xmid(this%nverts))
382 allocate (ymid(this%nverts))
383 allocate (lm(this%nverts))
384 allocate (umx(this%nverts))
385 allocate (umy(this%nverts))
386 allocate (kappax(this%nverts))
387 allocate (kappay(this%nverts))
388 allocate (vm0x(this%nverts))
389 allocate (vm0y(this%nverts))
390 allocate (vm1x(this%nverts))
391 allocate (vm1y(this%nverts))
392 allocate (unextxnext(this%nverts))
393 allocate (unextynext(this%nverts))
394 allocate (wk1(this%nverts))
395 allocate (wk2(this%nverts))
400 wk1 = this%xvertnext - this%xvert
401 wk2 = this%yvertnext - this%yvert
402 le = dsqrt(wk1 * wk1 + wk2 * wk2)
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
418 do i = 1, this%nverts
419 xvals(1) = this%xvert(i)
420 xvals(2) = this%xvertnext(i)
422 yvals(1) = this%yvert(i)
423 yvals(2) = this%yvertnext(i)
425 areasub(i) = area(xvals, yvals)
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)
435 divcell = sum(le * this%vne) / areacell
438 wk1 = this%xvert - this%xctr
439 wk2 = this%yvert - this%yctr
440 li = dsqrt(wk1 * wk1 + wk2 * wk2)
444 unextxnext = cshift(unintx, 1)
445 unextynext = cshift(uninty, 1)
448 xmid = 5.d-1 * (this%xvert + this%xctr)
449 ymid = 5.d-1 * (this%yvert + this%yctr)
452 wk1 = cshift(xmid, 1) - xmid
453 wk2 = cshift(ymid, 1) - ymid
454 lm = dsqrt(wk1 * wk1 + wk2 * wk2)
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)
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)
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)
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
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
521 deallocate (unextxnext)
522 deallocate (unextynext)