350 class(MethodCellTernaryType),
intent(inout) :: this
355 real(DP),
allocatable,
dimension(:) :: xvals
356 real(DP),
allocatable,
dimension(:) :: yvals
363 real(DP),
allocatable,
dimension(:) :: wk1
364 real(DP),
allocatable,
dimension(:) :: wk2
365 real(DP),
allocatable,
dimension(:) :: unextxnext
366 real(DP),
allocatable,
dimension(:) :: unextynext
367 real(DP),
allocatable,
dimension(:) :: le
368 real(DP),
allocatable,
dimension(:) :: unextx
369 real(DP),
allocatable,
dimension(:) :: unexty
371 real(DP),
allocatable,
dimension(:) :: areasub
373 real(DP),
allocatable,
dimension(:) :: li
374 real(DP),
allocatable,
dimension(:) :: unintx
375 real(DP),
allocatable,
dimension(:) :: uninty
376 real(DP),
allocatable,
dimension(:) :: xmid
377 real(DP),
allocatable,
dimension(:) :: ymid
378 real(DP),
allocatable,
dimension(:) :: lm
379 real(DP),
allocatable,
dimension(:) :: umx
380 real(DP),
allocatable,
dimension(:) :: umy
381 real(DP),
allocatable,
dimension(:) :: kappax
382 real(DP),
allocatable,
dimension(:) :: kappay
383 real(DP),
allocatable,
dimension(:) :: vm0x
384 real(DP),
allocatable,
dimension(:) :: vm0y
385 real(DP),
allocatable,
dimension(:) :: vm1x
386 real(DP),
allocatable,
dimension(:) :: vm1y
388 integer(I4B) :: nvert
390 select type (cell => this%cell)
391 type is (cellpolytype)
394 allocate (le(this%nverts))
395 allocate (unextx(this%nverts))
396 allocate (unexty(this%nverts))
397 allocate (areasub(this%nverts))
398 allocate (li(this%nverts))
399 allocate (unintx(this%nverts))
400 allocate (uninty(this%nverts))
401 allocate (xmid(this%nverts))
402 allocate (ymid(this%nverts))
403 allocate (lm(this%nverts))
404 allocate (umx(this%nverts))
405 allocate (umy(this%nverts))
406 allocate (kappax(this%nverts))
407 allocate (kappay(this%nverts))
408 allocate (vm0x(this%nverts))
409 allocate (vm0y(this%nverts))
410 allocate (vm1x(this%nverts))
411 allocate (vm1y(this%nverts))
412 allocate (unextxnext(this%nverts))
413 allocate (unextynext(this%nverts))
414 allocate (wk1(this%nverts))
415 allocate (wk2(this%nverts))
420 wk1 = this%xvertnext - this%xvert
421 wk2 = this%yvertnext - this%yvert
422 le = dsqrt(wk1 * wk1 + wk2 * wk2)
427 areacell = area(this%xvert, this%yvert)
428 sixa = areacell * dsix
429 wk1 = -(this%xvert * this%yvertnext - this%xvertnext * this%yvert)
430 nvert =
size(this%xvert)
431 this%xctr = sum((this%xvert + this%xvertnext) * wk1) / sixa
432 this%yctr = sum((this%yvert + this%yvertnext) * wk1) / sixa
438 do i = 1, this%nverts
439 xvals(1) = this%xvert(i)
440 xvals(2) = this%xvertnext(i)
442 yvals(1) = this%yvert(i)
443 yvals(2) = this%yvertnext(i)
445 areasub(i) = area(xvals, yvals)
449 term =
done / (cell%defn%porosity * cell%defn%retfactor * this%dz)
450 do i = 1, this%nverts
451 this%vne(i) = -cell%defn%faceflow(i) * term / le(i)
455 divcell = sum(le * this%vne) / areacell
458 wk1 = this%xvert - this%xctr
459 wk2 = this%yvert - this%yctr
460 li = dsqrt(wk1 * wk1 + wk2 * wk2)
464 unextxnext = cshift(unintx, 1)
465 unextynext = cshift(uninty, 1)
468 xmid = 5.d-1 * (this%xvert + this%xctr)
469 ymid = 5.d-1 * (this%yvert + this%yctr)
472 wk1 = cshift(xmid, 1) - xmid
473 wk2 = cshift(ymid, 1) - ymid
474 lm = dsqrt(wk1 * wk1 + wk2 * wk2)
491 call this%calc_thru_hcsum(vm0i0, divcell, le, li, lm, areasub, areacell, &
492 unintx, uninty, unextx, unexty, &
493 unextxnext, unextynext, &
494 kappax, kappay, vm0x, vm0y, vm1x, vm1y, hcsum0)
496 vm0ival = vm0i0 + perturb
497 call this%calc_thru_hcsum(vm0ival, divcell, le, li, lm, areasub, areacell, &
498 unintx, uninty, unextx, unexty, &
499 unextxnext, unextynext, &
500 kappax, kappay, vm0x, vm0y, vm1x, vm1y, hcsum)
502 jac = (hcsum - hcsum0) / perturb
503 vm0ival = vm0i0 - hcsum0 / jac
504 call this%calc_thru_hcsum(vm0ival, divcell, le, li, lm, areasub, areacell, &
505 unintx, uninty, unextx, unexty, &
506 unextxnext, unextynext, &
507 kappax, kappay, vm0x, vm0y, vm1x, vm1y, hcsum)
512 this%vv0x = 2.d0 * vm0x - this%vctrx
513 this%vv0y = 2.d0 * vm0y - this%vctry
514 this%vv1x = 2.d0 * vm1x - this%vctrx
515 this%vv1y = 2.d0 * vm1y - this%vctry
518 term =
done / (cell%defn%retfactor * cell%defn%porosity * areacell)
519 this%vzbot = cell%defn%faceflow(this%nverts + 2) * term
520 this%vztop = -cell%defn%faceflow(this%nverts + 3) * term
541 deallocate (unextxnext)
542 deallocate (unextynext)